# # 12_09_06 12_09_07 12_09_10 12_09_11 12_09_12 # 12_09_13 12_10_30 13_03_19 13_03_21 13_03_25 # 13_03_26 13_03_27 13_06_11 15_04_10 15_07_08 # # confection des figures pour le chapitre 3 # source("betbin.code.r"); source("nornor.code.r"); source("rbsb1.code.r"); source("rbsb2.code.r"); source("rbsb3.code.r"); source("rbsb4.code.r"); rbsb3k("RESET"); # chemin <- "./" nofi <- function(x) { paste(chemin,x,".pdf",sep="");} voir <- TRUE; voir <- FALSE; cgre <- "white"; cora <- "darkgrey"; cred <- "black"; cyel <- "white"; cpin <- "grey"; # values N <- 24; r <- 10; # vraisemblance 1 if (!voir) { pdf(nofi("vrais1"),width=7,height=10);} par(mfrow=c(1,1)); pp <- seq(0.1,0.9,0.1); draw8binos(n=N,p=pp, sub="",xlab="",ylab="", #main=paste(paste("Pour \\pi =",paste(pp,collapse=" / ")), # "\n","lorsque N = ",N," et r = ",r, # sep="") ccol=c(cred,cora,cgre), main ="" ); abline(v=r,lwd=3,col=cred); text(22,0.20,labels=expression(" "*pi*" = 0.1"),cex=1.5); text(22,0.50,labels=expression(" "*pi*" = 0.2"),cex=1.5); text(22,0.75,labels=expression(" "*pi*" = 0.3"),cex=1.5); # text(2,1.00,labels=expression(" "*pi*" = 0.4"),cex=1.5); text(2,1.25,labels=expression(" "*pi*" = 0.5"),cex=1.5); text(2,1.50,labels=expression(" "*pi*" = 0.6"),cex=1.5); text(2,1.75,labels=expression(" "*pi*" = 0.7"),cex=1.5); text(2,2.00,labels=expression(" "*pi*" = 0.8"),cex=1.5); text(2,2.27,labels=expression(" "*pi*" = 0.9"),cex=1.5); if (!voir) { dev.off(); } else { pause("vrais1");} # test d'hypothèse if (!voir) { pdf(nofi("test"),width=7,height=6);} draw8binos(n=24,pp=0.50, sub="",xlab="",ylab="", main=expression("Distribution de R lorsque "*pi*" = 0.5"), ccol=c(cred,cora,cgre), ci=0.95,xlimi=c(3,21)); if (!voir) { dev.off(); } else { pause("test");} # intervalle de confiance if (!voir) { pdf(nofi("intconf"),width=7,height=10);} pp <- seq(0.1,0.9,by=0.1); part2 <- paste(" = ", paste(pp,collapse=" / "),"\")\n lorsque N = ",N, " et r = ",r, sep=""); afaire <- paste("mainti <- expression(\"Pour \"*pi*\"",part2, sep=""); # perte de temps ! #eval(parse(text=afaire)); draw8binos(n=24,pp=pp, sub="",xlab="",ylab="", ccol=c(cred,cora,cgre), main="Intervalle de confiance pour une variable discrète", ci=0.95,xlimi=c(0,24)); abline(v=10,lwd=3,col=cred); text(22,0.20,labels=expression(" "*pi*" = 0.1"),cex=1.5); text(22,0.50,labels=expression(" "*pi*" = 0.2"),cex=1.5); text(22,0.75,labels=expression(" "*pi*" = 0.3"),cex=1.5); # text(2,1.00,labels=expression(" "*pi*" = 0.4"),cex=1.5); text(2,1.25,labels=expression(" "*pi*" = 0.5"),cex=1.5); text(2,1.50,labels=expression(" "*pi*" = 0.6"),cex=1.5); text(2,1.75,labels=expression(" "*pi*" = 0.7"),cex=1.5); text(2,2.00,labels=expression(" "*pi*" = 0.8"),cex=1.5); text(2,2.27,labels=expression(" "*pi*" = 0.9"),cex=1.5); if (!voir) { dev.off(); } else { pause("intconf");} #### # new figures for the normal example #### # distribution jointe # yobs <- 60; sigy <- 5; muth <- 55; sdth <- 8; ymuth <- function(theta) {theta;}; # inverting it thmuy <- function(y) {(64*y+25*55)/(64+25);}; muy <- muth; sdy <- sqrt(sdth^2+sigy^2); rho <- sqrt(sdth^2 / (sdth^2 + sigy^2)); sigth <- sqrt(sdth^2*(1-rho^2)); titre <- ""; if (!voir) { pdf(nofi("jointe1-norm"),width=7,height=6);} thlab= expression(" "*mu*" "); etendue <- c(20,80); draw8nornor(mean.x=muy,sd.x=sdy, mean.y=thmuy,sd.y=sigth, ci=-1, colo=c(cgre,cred), xran=etendue,yran=etendue, xlab="Y",ylab=thlab ); abline(v=muy,h=thmuy(muy),lwd=2); abline(v=muy+1.96*sdy*c(-1,1),lty=2); abline(h=thmuy(muy)+1.96*sdth*c(-1,1),lty=2); abline(v=60+0.4*c(-1,1),col="grey",lwd=3); abline(h=thmuy(60),lwd=1,col=cred); if (!voir) { dev.off(); } else { pause("jointe1-norm");} # vraisemblance 1 par(mfrow=c(1,1)); if (!voir) { pdf(nofi("vrais1-norm"),width=7,height=10);} mu <- seq(30,70,5); sigma <- 5; yobs <- 60; para <- data.frame(mean=mu,sd=sigma); papa <- paste("mu=",mu,sep=""); papa[1] <- expression(" "*mu*"=30"); for (ii in bf(mu)) { afaire <- paste("papa[ii] <- expression(\" \"*mu*\"=",mu[ii],"\")",sep=""); eval(parse(text=afaire)); } draw8cprobas(para,c(10,90),dnorm,qnorm, yval=mu,palab=papa,sipalab=0.1, ci=-1,sti="",xlab="",ylab="", mti="",kechelle=10,grafo=c("1","2","white")); abline(v=yobs,lwd=3,col=cred); if (!voir) { dev.off(); } else { pause("vrais1-norm");} # test d'hypothèse if (!voir) { pdf(nofi("test-norm"),width=7,height=6);} xli <- c(30,80); yli <- c(0,1); ou <- cbind(xli,yli); plot(0:1,0:1,type="n",xlim=xli,ylim=yli,xlab="",ylab=""); draw8cproba(as.vector(ou),c(mean=55,sd=5),varange=xli, dnorm,qnorm,vax=NULL,ci=0.95, echelle=1,orientation="N", grafo=c("1","2",cgre,cred) ); abline(v=55); if (!voir) { dev.off(); } else { pause("test-norm");} # Intervalle de confiance if (!voir) { pdf(nofi("intconf-norm"),width=7,height=10);} mu <- seq(30,70,5); sigma <- 5; yobs <- 60; para <- data.frame(mean=mu,sd=sigma); papa <- paste("mu=",mu,sep=""); papa[1] <- expression(" "*mu*"=30") ; for (ii in bf(mu)) { afaire <- paste("papa[ii] <- expression(\" \"*mu*\"=",mu[ii],"\")",sep=""); eval(parse(text=afaire)); } draw8cprobas(para,c(10,90),dnorm,qnorm, yval=mu,palab=papa,sipalab=0.1, ci=0.95,sti="",xlab="",ylab="", mti="Intervalle de confiance pour une variable continue", kechelle=10,grafo=c("1","2",cgre,cred)); abline(v=yobs,lwd=3,col=cred); if (!voir) { dev.off(); } else { pause("intconf-norm");} ################################################## ################################################## # # Approximation linéaire # if (!voir) { pdf(nofi("approline"),width=7,height=5);} Yval <- function(x,co) { res <- co[1]; if (length(co) > 1) { for (dd in 2:length(co)) { res <- res + x^(dd-1)*co[dd]; } } res; } Ytit <- function(co) { res <- paste("Y =",co[1]); if (length(co) > 1) { for (dd in 2:length(co)) { if (dd == 2) { xx <- "*x"; dd1 <- ""; } else { xx <- "*x^"; dd1 <- as.character(dd-1); } if (co[dd] > 0) { res <- paste(res,"+",paste(co[dd],xx,dd1,"",sep="")); } if (co[dd] < 0) { res <- paste(res,"-",paste(-co[dd],xx,dd1,"",sep="")); } } } res; } deriv <- function(co) { res <- rep(NA,length(co)-1); for (ii in 1:length(res)) { res[ii] <- co[ii+1]*ii; } res; } a0 <- c(385,2,-10,1); a1 <- deriv(a0); nbp <- 100; xx <- seq(-5,10,length=100); yy <- Yval(xx,a0); plot(xx,yy,type="l",col=cred, lwd=6,xlab="x", #ylab=Ytit(a0), ylab="f(x)", cex.lab=1.4,main=""); ou <- c(-3,0,5,9); for (oo in ou) { y0 <- Yval(oo,a0); z0 <- Yval(oo,a1); cc <- c(y0-oo*z0,z0); zz <- Yval(xx,cc); lines(xx,zz); lines(c(oo,oo),c(-100,Yval(oo,a0)),lty=2); } if (!voir) { dev.off(); } else { pause("approline");} # # modèle 1 # if (!voir) { pdf(nofi("modele1"),width=7,height=5);} nbp <- 100; tt <- seq(15,25,length=nbp); pro <- c(0.005,0.025,0.25,0.75,0.975,0.995); colo <- c("grey95","grey90","grey60","grey90","grey95"); nbcol <- length(colo); y1 <- function(x,co) { co[1] + (x-co[2])^2*co[3];} co <- c(10,18,-0.1,0.7); yy <- y1(tt,co); cil <- matrix(NA,nbp,nbcol+1); for (ii in 1:nbp) { mu <- yy[ii]; sig <- co[4]; cil[ii,] <- qnorm(pro,mu,sig); } plot(tt,yy,type="l", ylim =range(c(cil,0)), lwd=2.5,xlab="Température",ylab="Croissance", cex.lab=1.3,main=""); for (ii in 1:nbcol) { polygon(c(tt,rev(tt)),c(cil[,ii],rev(cil[,ii+1])), col=colo[ii],border=NA); } lines(tt,yy,lwd=5,lty=3); if (!voir) { dev.off(); } else { pause("modele1");} # # les modèles 1 # if (!voir) { pdf(nofi("modeles1"),width=7,height=10);} par(mfrow=c(2,2)); coc <- matrix(c(10,18,-0.1,0.4, 12,18,-0.1,1.5, 10,21,-0.03,0.7, 8,22,0.1,0.9),4,4); for (mm in 1:4) { co <- coc[,mm]; yy <- y1(tt,co); cil <- matrix(NA,nbp,nbcol+1); for (ii in 1:nbp) { mu <- yy[ii]; sig <- co[4]; cil[ii,] <- qnorm(pro,mu,sig); } plot(tt,yy,type="l", ylim = c(0,15), lwd=2.5,xlab="Température",ylab="Croissance", cex.lab=1.3,main=paste("(",letters[mm],")")); for (ii in 1:nbcol) { polygon(c(tt,rev(tt)),c(cil[,ii],rev(cil[,ii+1])), col=colo[ii],border=NA); } lines(tt,yy,lwd=2,lty=3); } if (!voir) { dev.off(); } else { pause("modeles1");} # # Les figures pour la vraisemblance # # distributions Bêtas if (!voir) { pdf(nofi("betas"),width=7,height=10);} vv <- c(0.5,1,2,4,8,16); xx <- seq(0,1,length=200)[-c(1,200)]; par(mfrow=c(length(vv),length(vv)),mai=rep(0.18,4)); for (aa in vv) { for (bb in vv) { mainti <- paste("a =",aa,"; b =",bb); plot(xx,dbeta(xx,aa,bb),type="l", lwd=2,xlab="",ylab="",sub="", xaxt="n",yaxt="n", main=""); }} if (!voir) { dev.off(); } else { pause("betas");} # distribution jointe # # par rapport à p aa <- 4; bb <- 2; nn <- 24; rr <- 10; pp <- seq(0.1,0.9,0.1); titre <- paste("a =",aa,"; b =",bb,"; N =",nn); if (!voir) { pdf(nofi("jointe1"),width=7,height=6);} par(mfrow=c(1,1)); draw8betabino(n=nn,a=aa,b=bb, pp=pp, rr=NULL, mti="", xlab="",ylab="", sti="", type="joint", format="complete"); if (!voir) { dev.off(); } else { pause("jointe1");} # # par rapport à R if (!voir) { pdf(nofi("jointe2"),width=7,height=6); } draw8betabino(n=nn,a=aa,b=bb, pp=NULL, rr=0:nn, #mti=titre, mti="", xlab="",ylab="", sti="", type="joint", format="complete"); if (!voir) { dev.off(); } else { pause("jointe2");} # Distribution postériore titre <- paste(titre,"; r =",rr); if (!voir) { pdf(nofi("posteriore1"),width=7,height=6);} draw8betabino(n=nn,a=aa,b=bb, pp=NULL, rr=rr, #mti=titre, mti="", xlab="",ylab="", sti="",ci=0.95, type="joint", ccol=c(cred,cora,cgre), format="complete"); abline(v=rr-0.45,lwd=5,col=cred); if (!voir) { dev.off(); } else { pause("posteriore1");} # # comparaison de la priore et de la postériore aaa <- aa+r; bbb <- bb+N-r; xx <- seq(0,1,le=100); y1 <- dbeta(xx, aa, bb); y2 <- dbeta(xx,aaa,bbb); if (!voir) { pdf(nofi("posteriore2"),width=7,height=6);} par(cex.axis=1.2,cex.lab=2); plot(xx,y2,type="l", main="", xlab=expression(""*pi*""),ylab="" ); abline(h=0); lines(xx,y1,lty=2); if (!voir) { dev.off(); } else { pause("posteriore2");} # # faisons varier r if (!voir) { pdf(nofi("posteriore3"),width=7,height=6);} par(cex.axis=1.2,cex.lab=2); plot(xx,y2,type="l", ylim=c(0,8), main="", xlab=expression(""*pi*""),ylab="" ); abline(h=0); for (rrr in 0:N) { aaa <- aa+rrr; bbb <- bb+N-rrr; yy <- dbeta(xx,aaa,bbb); lines(xx,yy,lty=2); } if (!voir) {dev.off(); } else { pause("posteriore3");} # # faisons varier N à r/N fixé if (!voir) { pdf(nofi("posteriore4"),width=7,height=6);} par(cex.axis=1.2,cex.lab=2); aaa <- aa+rr; bbb <- bb+N-rr; y2 <- dbeta(xx,aaa,bbb); plot(xx,y2,type="l", ylim=c(0,10), main="", xlab=expression(""*pi*""),ylab="" ); abline(h=0); for (n in 1:10) { aaa <- aa+n*5; bbb <- bb+n*7; yy <- dbeta(xx,aaa,bbb); lines(xx,yy,lty=2); } if (!voir) { dev.off(); } else { pause("posteriore4");} # # analyse séquentielle library(igraph); dag1 <- graph.formula(P1-+D1, D1-+P2, P2-+D2, D2-+P3, P3-+D3, D3-+P4, PP-+DD, DD-+RR,E1,E2); lay <- matrix(c(1:7,1,4,7,4,4, 7:1,rep(10,3),8.5,2),12); col <- c(cyel,cgre,cora,cgre,cora,cgre,cgre, cyel,cgre,cgre,NA,NA); sha <- c("circle","square","circle","square","circle","square","circle", "circle","rectangle","circle"); V(dag1)$label <- c("Pr","D1","PP1","D2","PP2","D3","Pos", "Pr","D1+D2+D3","Pos","(global)","(séquentiel)"); if (!voir) { pdf(nofi("bayeseq"),width=7,height=7);} par(cex=1.2); plot(dag1,layout=lay,frame=FALSE, main="",edge.color="black", edge.color="black", vertex.shape=sha, vertex.color=col, vertex.size=c(rep(20,8),60,20,0,0) ); text(8,4,labels="(global)"); if (!voir) { dev.off(); } else { pause("bayeseq");} # # Prédictive dag2b <- graph.formula(Y+-Theta, Theta-+Z); lay <- matrix(c(1,2,3,1,2,1),ncol=2); lay[,1] <- lay[,1]*0.5; g1 <- cgre; g2 <- cyel; g3 <- cpin; col <- c(g1,g2,g3); sha <- c("square",rep("circle",2)); V(dag2b)$label <- V(dag2b)$name; if (!voir) { pdf(nofi("predictive"),width=7,height=7);} par(cex=1.2); plot(dag2b,layout=lay,frame=FALSE, main="",edge.color="black", edge.color="black", vertex.shape=sha, vertex.color=col, vertex.size=c(36.5,40,40) ); if (!voir) { dev.off(); } else { pause("predictive");} # # # Méta-analyse dag2 <- graph.formula(Theta-+D1, Theta-+D2, Theta-+D3, Theta-+D4, Theta-+D5); lay <- matrix(c(0,cos(2*pi*(1:5)/5), 0,sin(2*pi*(1:5)/5)),6); g1 <- "darkgrey"; g2 <- "lightgrey"; g3 <- "grey"; col <- c(cyel,g1,g2,g1,g3,g2); sha <- c("circle",rep("square",5)); V(dag2)$label <- V(dag2)$name; if (!voir) { pdf(nofi("bayesme"),width=7,height=7);} par(cex=1.2); plot(dag2,layout=lay,frame=FALSE, main="",edge.color="black", edge.color="black", vertex.shape=sha, vertex.color=col, vertex.size=c(60,20,40,25,30,20) ); text(8,4,labels="(global)"); if (!voir) { dev.off(); } else { pause("bayesme");} # # library(binom); print(binom.confint(10,24)); # # stop("All Figures were made");