"age" <- function (U) { if(missing(U)) U<-sample(deck) ld<-length(U) opars<-par(bg='lightgray') plot(c(0,ld+1),c(0,11),type='n',xlab='cards',ylab='age') xy<-par('usr') rect(xy[1],xy[3],xy[2],xy[4],col='black') cols<-rainbow(10) for(start in 1:10) { this<-start while(this<=ld) { card<-U[this] this<-this+card cat(this,' ') lines(c(this-card,this,this),c(0,card,0),col=cols[start],lwd=2) } cat('\n') cat(start,card,ld-(this-card),'\n') } invisible(par(opars)) } "bd" <- function (lam=c(rep(8,20),0),mu=0:20,T=1) { S<-length(lam)-1 nu<-c(1,cumprod(lam[-(S+1)]/mu[-1])) nu<-nu/sum(nu) sc<-0.1*T/max(nu) seq<-sample(0:S,1,prob=nu) opars<-par(bg='lightgray') plot(T*c(-0.1,1),c(0,S),type='n',xlab='time',ylab='state') xy<-par('usr'); rect(xy[1],xy[3],xy[2],xy[4],col='black') lines(c(0,0),c(0,S),col='white') points(-sc*nu,0:S,col='white') z1<-runbd(0,floor(S/4),T,lam,mu) lines(z1$pathx,z1$pathy,col='red',lwd=2) v1<-z1$visit z2<-runbd(0,seq,T,lam,mu) lines(z2$pathx,z2$pathy,col='green',lwd=2) v2<-z2$visit Thit<-2*T for(s in 0:(min(length(v1),length(v2))-1)) { v1s<-v1[[s+1]]; v2s<-v2[[s+1]] if(!is.null(v1s)&!is.null(v2s)) { v1s<-matrix(v1s,ncol=2,by=TRUE); v2s<-matrix(v2s,ncol=2,by=TRUE) for(i in 1:nrow(v1s)) for(j in 1:nrow(v2s)) if(!((v1s[i,2]v2s[j,2]))) { t<-max(v1s[i,1],v2s[j,1]) if(tThit lines(c(Thit,z2$pathx[i]),c(sthit,z2$pathy[i]),col='blue',lwd=2,lty=2) } invisible(par(opars)) } "bds" <- function (lam=c(rep(8,20),0),mu=0:20,T=1) { S<-length(lam)-1 nu<-c(1,cumprod(lam[-(S+1)]/mu[-1])) nu<-nu/sum(nu) sc<-0.1*T/max(nu) seq<-sample(0:S,1,prob=nu) opars<-par(bg='lightgray') plot(T*c(-0.1,1),c(0,S),type='n',xlab='time',ylab='state') xy<-par('usr'); rect(xy[1],xy[3],xy[2],xy[4],col='black') lines(c(0,0),c(0,S),col='white') points(-sc*nu,0:S,col='white') for(j in 1:3) { z1<-runbd(0,floor(j*S/4),T,lam,mu) lines(z1$pathx,z1$pathy,col=c('red','yellow','blue')[j],lwd=2) } invisible(par(opars)) } "deck" <- c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 5, 5, 5, 5, 6, 6, 6, 6, 7, 7, 7, 7, 8, 8, 8, 8, 9, 9, 9, 9, 10, 10, 10, 10) "deck2" <- c(1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 5, 5, 5, 5, 6, 6, 6, 6, 7, 7, 7, 7, 8, 8, 8, 8, 9, 9, 9, 9, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10) "gcbp" <- function(p=.9,x=0:5) { y<-dbinom(x,1,p) z<-dpois(x,p) r<-pmin(y,z) opars<-par(bg='lightgray') plot(c(min(x)-1,max(x)+1),range(y,z),type='n',xlab='value',ylab='probability') xy<-par('usr'); rect(xy[1],xy[3],xy[2],xy[4],col='black') for(i in 1:length(x)) { rect(x[i]-.3,0,x[i],y[i],col='green') rect(x[i],0,x[i]+.3,z[i],col='blue') rect(x[i]-.3,0,x[i]+.3,r[i],den=7,col='red',lwd=3) } title('Gamma coupling of Bernoulli and Poisson') repeat { dd<-readline() if(dd!='') break j<-sample(length(x),1,prob=y); u<-runif(1) if(u*y[j]lambda] tv<-n for(i in 1:length(n)) tv[i]<-sum(abs(dbinom(0:200,n[i],lambda/n[i])-dpois(0:200,lambda))) bound<-2*lambda^2/n lecam<-2*lambda/n plot(n,bound,log='xy',type='n',ylim=range(tv,bound),ylab='distance') lines(n,bound,lwd=2) if(lambda>1) lines(n,lecam,lwd=2,col='blue',lty=2) points(n,tv,col='red',cex=2,pch=16) title(main = substitute(paste('||Binomial(n,',lambda, '/n)-Poisson(',lambda,')|| for ',lambda,'=',lam),list(lam=lambda))) } "rebd" <- function (z1,z2,lam=c(rep(8,20),0),mu=0:20,T=1) { S<-length(lam)-1 nu<-c(1,cumprod(lam[-(S+1)]/mu[-1])) nu<-nu/sum(nu) sc<-0.1*T/max(nu) opars<-par(bg='lightgray') plot(T*c(-0.1,1),c(0,S),type='n',xlab='time',ylab='state') xy<-par('usr'); rect(xy[1],xy[3],xy[2],xy[4],col='black') lines(c(0,0),c(0,S),col='white') points(-sc*nu,0:S,col='white') lines(z1$pathx,z1$pathy,col='red',lwd=2) v1<-z1$visit lines(z2$pathx,z2$pathy,col='green',lwd=2) v2<-z2$visit Thit<-2*T for(s in 0:(min(length(v1),length(v2))-1)) { v1s<-v1[[s+1]]; v2s<-v2[[s+1]] if(!is.null(v1s)&!is.null(v2s)) { v1s<-matrix(v1s,ncol=2,by=TRUE); v2s<-matrix(v2s,ncol=2,by=TRUE) print(v1s) print(v2s) for(i in 1:nrow(v1s)) for(j in 1:nrow(v2s)) { cat(v1s[i,1],v1s[i,2],v2s[j,1],v2s[j,2]) if(!((v1s[i,2]v2s[j,2]))) { cat(' hit') t<-max(v1s[i,1],v2s[j,1]) if(tThit lines(c(Thit,z2$pathx[i]),c(sthit,z2$pathy[i]),col='blue',lwd=2,lty=2) } invisible(par(opars)) } "runbd" <- function (t,s,T,lam,mu) { visit<-list() pathx<-t; pathy<-s while(ts) visit[[s+1]]<-c(visit[[s+1]],t,tn) else visit[[s+1]]<-c(t,tn) s<-sn t<-tn } list(visit=visit,pathx=pathx,pathy=pathy) } "xs" <- function (U) { if(missing(U)) U<-sample(deck) ld<-length(U) cmx<-max(U) opars<-par(bg='lightgray') plot(c(0,ld+1),c(0,1+cmx),type='n',xlab='cards',ylab='excess life') xy<-par('usr') rect(xy[1],xy[3],xy[2],xy[4],col='black') cols<-rainbow(cmx) for(start in 1:cmx) { this<-start while(this<=ld) { card<-U[this] this<-this+card cat(this,' ') lines(c(this-card,this-card,min(this,ld)), c(0,card,max(0,this-ld)),col=cols[start],lwd=2) } cat('\n') cat(start,card,ld-(this-card),'\n') } invisible(par(opars)) }