# Poisson counts and exponential waiting times # in case of failures par(mfrow=c(1,1),mar = c(0.1,0.1,0.1,0.1)) CEX.g =0.7 LWD=8 n.trials = 25 W=600; DW=100 XLIM=c(-W/30,3*W); DY=0.65 DY.0=0.25 YLIM=c(-n.trials/10,n.trials*(1+DY)) Y.0 = (1+DY.0)*n.trials; mu=80 # mu of time r.v dt=1 DT=100 n.items=12 RAINBOW=rainbow(n.items) #hard RAINBOW=c( "#47CF69", "#E41C83", "#EF9949", "#1E5994", "#6DFDE4", "#2B5D25", "#83F805", "#FDCEBA", "#653011", "#FF8FEE", "#861379", "#BEF56E") #soft RAINBOW=c( "#6BCC4E", "#A85AC8", "#B84836", "#4E5B34", "#6998C1", "#7ECCA4", "#51354D", "#C6487B", "#C8B09B", "#C4893A", "#C2C84D", "#BF8FBF") n=n.items*n.trials x.origin = 1.23*W ; d.x= (1.5*W)/W y.origin = 0 ; d.y = n.trials set.seed(1234567) # 4 set.seed(2234567) # 1 set.seed(3234567) # 1 set.seed(4234567) # 2 set.seed(5234567) # 2 set.seed(6234567) # 2 set.seed(7234567) # 2 set.seed(8234567) # 1 set.seed(9234567) # 2 set.seed(9876543) # 1 set.seed(8876543) # 5 set.seed(7876543) # 0 set.seed(6876543) # 3 set.seed(5876543) # 1 set.seed(4876543) # 4 set.seed(3876543) # 1 set.seed(2876543) # 2 set.seed(1876543) # 2 set.seed(0876543) # 2 set.seed(9999999) # 0 set.seed(8888888) # 0 set.seed(7777777) # 3 set.seed(6666666) # 1 set.seed(5555555) # 0 set.seed(4444444) # 1 set.seed(3333333) # 3 plot( XLIM, YLIM, col="white",xlab="Km", ylab="Proportion in g 0, 1, 2, ... ") #y = (1:n.trials)/(n.trials+1) #for (k in 1:n.items) lines(qgamma(y,k,1/mu),y*n.trials) # width of time window (in Km rather than time) interval = matrix(rexp(n,rate=1/mu),n.trials,n.items) t.end = matrix(NA,n.trials,n.items) ; t.start=t.end n.FAIL = rep(0,n.trials) x.l=3*XLIM[1] text(x.l,n.trials/2,"A",adj=c(0,0.5)) for(r in 1:n.trials) { t.end[r,] = cumsum(interval[r,]) t.start[r,] = c(0,t.end[r,(1:(n.items-1))]) n.fail=0 y=r -0.5 for(i in 1:n.items) { if( t.start[r,i] < W ) { segments(t.start[r,i],y, min(t.end[r,i],W),y,col=RAINBOW[i], lwd=LWD, lend=1) if(t.end[r,i] < W) { points( min(t.end[r,i],W)-0.02,y,col="black",pch=19,cex=0.1) n.fail=n.fail+1 if(n.fail==n.items) points( min(t.end[r,i],W)+1, Y.0 - 4 - runif(1,0,0.5), col="black",pch=19,cex=0.2) } } } if(r==1) text(qgamma(0.001,n.items,1/mu),Y.0 - 4 - 1/2, "B:",adj=c(1,0.5)) if(n.fail < n.items) points( t.end[r,i], Y.0 - 4 - runif(1,0,0.5), col="black",pch=19,cex=0.2) if(n.fail < n.items) text(1.06*W,y,toString(n.fail+1),adj=c(1,0.5), col=RAINBOW[n.fail+1], cex=CEX.g) if(r==1) text(1.09*W, n.trials/2, "lequel article en cours d'utilisation à la fin", adj=c(0.5,1), cex=0.75, srt=90) n.FAIL[r] = n.fail } for(t in seq(0,3*W,DW)) { segments(t,Y.0,t,Y.0-1,lwd=0.5,col="lightblue") text(t, Y.0-1.5,toString(t), adj=c(0.5,1.25),col="lightblue",cex=0.5) if( t <= W) text(t, y.origin-0.025*d.y,toString(t), adj=c(0.5,1),col="lightblue",cex=0.5) } t.range=seq(0,(3*W),20) DY = n.trials*(DY-DY.0)*1.00 pdf = dgamma(t.range,shape=3,scale=mu) max.pdf = max(pdf) for ( k in c(12,9,6,3,18)) { pdf = DY*dgamma(t.range,shape=k,scale=mu)/max.pdf maX = max(pdf) at.t = t.range[which(pdf==maX)] if(k != n.items) lines(t.range,Y.0+ pdf, lty="dotted",lwd=0.75,col="grey20") text(at.t, Y.0+maX,toString(k), adj=c(0,-0.25), cex=0.5) if(k==n.items) { polygon(c(0,t.range[t.range<=W],W), c(Y.0,Y.0+ pdf[t.range<=W],Y.0), col="grey85", border=NA) polygon(c(W,t.range[t.range>=W],3*W), c(Y.0,Y.0+ pdf[t.range>=W],Y.0), col="grey65", border=NA) P = pgamma(W,n.items,1/mu,lower.tail=FALSE) w.mid = qgamma((1-P)+P/2, n.items,1/mu) text(w.mid,Y.0+1.75*maX, paste(toString( round(100*P,1) ),"%",sep=""), adj=c(0,0), col="grey65") segments(w.mid,Y.0+1.65*maX, 0.92*w.mid, Y.0+0.4*maX, lwd=1, col="grey50") w.mid = qgamma((1-P), n.items,1/mu) text(w.mid,Y.0+1.75*maX, paste(toString( round(100*(1-P),1) ),"%",sep=""), adj=c(0,0), col="grey85") segments(w.mid,Y.0+1.65*maX, 0.9*w.mid, Y.0+0.1*maX, lwd=1, col="grey50") text(x.l,Y.0+maX,"C",adj=c(0,0.5)) } } segments(0,Y.0,3*W,Y.0,col="grey80",lwd=1) X=2.9*W ; DX= 0.04*W r = 0 ; I=0 for(i in 0:(n.items-1)) { NF = sum(n.FAIL==i) if(NF>0) { for (j in 1:NF) { r=r+1 y = r-0.5 segments(X,y, X+DX, y, lwd=LWD, col=RAINBOW[i+1], lend=1) text(X+3*DX,y,toString(i+1),adj=c(1,0.5), col=RAINBOW[i+1], cex=CEX.g) } } } text(3*W,y.origin,"A'",adj=c(1,1.5)) ##### y=0:n.items t = seq(0,W,DW/5) nn=length(t) m=matrix(NA,nrow=length(y),ncol=length(t)) m[1,] = rep(0,length(t)) for(Y in y[-1]) m[Y+1,] = ppois(Y-1, t/mu) for(Y in y[-1]){ bottom = m[Y,] top = m[Y+1,length(t):1] TOP = -sort(-top) mid=((TOP+bottom)/2) mid.t=t[1:nn] target = mid[6+2*Y] e = mid-target at.which.t = mid.t[which(abs(e)==min(abs(e)))] xx= c(t,-sort(-t)) yy=c(bottom, top ) polygon(x.origin+xx*d.x, y.origin+yy*d.y , col=RAINBOW[Y],border=NA) xxx = x.origin+at.which.t*d.x yyy = y.origin+target*d.y rect(xxx-16, yyy-0.4,xxx+16,yyy+0.4,col="white",border=NA) text(xxx,yyy,toString(Y), cex=0.6, font=2,col=RAINBOW[Y]) if(Y==n.items) { segments( x.origin+W*d.x, y.origin+P*d.y, x.origin+1.06*W*d.x, y.origin+(P+0.015)*d.y) text( x.origin+1.06*W*d.x, y.origin+(P+0.015)*d.y, paste(toString( round(100*P,1) ),"%",sep=""), adj=c(-0.1,0.5), col="grey65") } } rect(x.origin,y.origin, x.origin+W*d.x, y.origin+1*d.y,lwd=0.5,border="grey30") text(x.origin+(W/2)*d.x,y.origin+1*d.y,"D",adj=c(0.5,-0.25)) text(x.origin, y.origin-0.025*d.y, "Jour:", adj=c(0,1), cex=0.5, col="lightblue") text(x.origin+ (W/2)*d.x, y.origin-0.110*d.y, "nombre d'événements attendu", adj=c(0.5,1), cex=0.75) for(T in seq(DT,W,DT)){ text(x.origin+T*d.x, y.origin-0.025*d.y,toString(T), adj=c(0.5,1), cex=0.5, col="lightblue") text(x.origin+T*d.x, y.origin-0.075*d.y,toString(T/mu), adj=c(0.5,1), cex=0.5) segments(x.origin+T*d.x, 0,x.origin+T*d.x,n.trials, cex=0.5,col="white",lwd=1) } for (p in seq(0,1,0.05)) segments( x.origin+W*d.x, p*n.trials, x.origin+W*d.x+5, p*n.trials, lwd=0.4, col="grey60" ) for (p in seq(0,1,0.1)) { segments( x.origin+W*d.x, p*n.trials, x.origin+W*d.x+10, p*n.trials, lwd=0.8, col="grey40" ) text(x.origin+W*d.x+13, p*n.trials, paste(toString(p*100),"%",sep=""), adj=c(0,0.5), cex=0.6, col="lightblue") segments( x.origin+0*d.x, p*n.trials, x.origin+W*d.x, p*n.trials, lwd=0.75, col="lightblue", lty="dotted" ) } #rect(0.98*x.origin, 10*YLIM[1],3.2*W, y.origin+1.035*d.y,lwd=3)