Project
A Stat667
Computer solution
Coding the coin toss experiment
The core function is to simulate either the N geometric random variates or the coin tossing game as stated and returns a value of either 1 (last coin toss is of a single coin) or 0 (last coin toss is not of a single coin). Four methods are to be compared for their efciencies, namely toss, toss2, toss3, and toss4 in R functions: "toss"<function(n, prob = 0.5) {if(n <= 1) return(n) repeat { n <- n - rbinom(1, n, prob) if(n <= 1) break } n } "toss2"<function(n, prob = 0.5) {if(n <= 1) return(n) Recall(n - rbinom(1, n, prob), prob) } "toss3"<function(n, prob = 0.5) {x <- rgeom(n, prob) if (sum(x==max(x)) == 1) 1 else 0 } "toss4"<function(n, prob = 0.5) {MAX <- -1 nmax <- 1 for(i in 1:n) { g <- rgeom(1, prob) if(g > MAX) { MAX <- g nmax <- 1 } else if(g == MAX) nmax <- nmax + 1 } if(nmax == 1) 1 else 0 }
1
The last two are using geometric distributions and the rst two are simulations of coin tosses. The comparison of the CPU times (run on a SUN Linux box) are done by the following R session:
# There will be n.trials trials. Mean CPU times of these n.trials # trials for each of the four methods and each of the various values of # N are to be compared and tabulated. cpu.time <- function(expr, gcFirst=TRUE) sum(system.time(expr, gcFirst)[-3]) n.trials <- 100 cpu <- matrix(0,n.trials,4) # the 4 columns store CPU times, in n.trials # trials, for the four methods colnames(cpu) <- paste("method",1:4) for(N in c(5000,10000,20000,50000)){ # for each value of N for(i in 1:n.trials){ # for each trial cpu[i,1] <- cpu.time(toss(N)) cpu[i,2] <- cpu.time(toss2(N)) cpu[i,3] <- cpu.time(toss3(N)) cpu[i,4] <- cpu.time(toss4(N)) } cat("N =",N,"\n\n") print(apply(cpu,2,mean)) # print the mean CPU times for the 4 # methods for current N }
The comparisons are summarized in Table 1. Note that methods 1 and 2 (toss) are the most efcient Table 1: Average CPU times (in hundredth seconds) for the four methods with various values of N. method N 1 2 3 4 5000 0.0001 0.0001 0.0020 0.0517 10000 0.0004 0.0004 0.0035 0.1021 20000 0.0002 0.0002 0.0072 0.2032 50000 0.0002 0.0001 0.0197 0.5099 (in CPU time) and the latter consumes more memory due to its recursive nature. Hence method 1 will be adopted in simulation. Method 4 gives old fashion element-wise computation. Method 3 use vectorized computation. However, the overhead of function evaluations (rgeom, max, and sum) may cause its slowdown.
The simulation
The coin toss trial will be repeated n.trials = 1000 times and thus obtain an estimate pN . This complete a simulation cycle. For each of the various values of N , n.sims = 100 simulations are performed. The n.sims = 100 values of pN for each N will be plotted against time. The R program named TOSS is given in Appendix C. The calculation of a pN can be vectorized to improve run time: "pN.hat"<function(N, n.trials, prob = 0.5, ptail = 1 - prob) 2
{N <- numeric(n.trials) bigN <- N > 1 repeat { if ((n <- sum(bigN)) == 0) break N[bigN] <- rbinom(n, N[bigN], ptail) bigN[bigN] <- N[bigN] > 1 } mean(N) } The comparison of CPU times between vectorized and non-vectorized versions is given in the R session named PNHAT in Appendix C. The vectorized version is approximately 1000 times faster. The R program named IMPROVE in Appendix C makes use of the vectorization computation. The plots for pN s (program in PLOT in Appendix C) are given in Figure 1 in Appendix A.
Calculation of the exact value of pN
We will use the coin tossing model. Let = P (head). The state space is {0, 1, . . . , N }. The state variable Xn is the number of coins landed with tails at nth toss with X0 denoting the initial number of coins (i.e., N ). Note that 0 is the absorbing state and that Xn+1 |Xn b(Xn , 1 ). Let pi be the probability that the last coin toss is of a single coin with i coins initially. Then p1 = 1 and that
i i1
pi =
j=1
Pij pj =
j=1
fi (j)pj + (1 )i pi ,
where fi () is the p.d.f. b(i, of 1 ). It follows then that
i1
pi =
j=1
fi (j)pj / 1 (1 )i .
Hence, if we know pj for j = 1, . . . , i 1, then we are able to calculate pi . Program this in the following R function: "pN"<function(n, prob = 1, phead = 0.5, ptail = 1 - phead) {ans <- numeric(n); k <- length(prob) ans[1:k] <- prob for(i in (k + 1):n) { j <- 1:(i-1) ans[i] <- sum(dbinom(j, i, ptail) * ans[j]) / (1 - ptaili) } ans }
3
Thus, the invocation probs <- pN(100) in R stores pN for N = 1, . . . , 100 in probs. Subsequently, the invocation probs <- pN(1000,probs) uses existing pN s, calculates pN for N = 101, . . . , 1000, and stores all pN s back to probs. The pN s for N = 2, . . . , 1000 are plotted in Figure 2 in Appendix B. The rectangles on the right of the panels indicate the same range in the y-axis. Thus, each latter panel gives a blowup view (in y-axis of the previous panel. The program used to generate this gure is given in Appendix D.
4
Appendix A
Figure 1: Simulation results for various values of N.
N = 10
0.76 0.76
N = 1000
0.72
0.68
0
20
40
60
80
100
0.68 0
0.72
pN
pN
20
40
60
80
100
simulation #
simulation #
N = 50
0.76 0.76
N = 2000
0.72
0.68
0
20
40
60
80
100
0.68 0
0.72
pN
pN
20
40
60
80
100
simulation #
simulation #
N = 100
0.76 0.76
N = 5000
0.72
0.68
0
20
40
60
80
100
0.68 0
0.72
pN
pN
20
40
60
80
100
simulation #
simulation #
N = 500
0.76 0.68 0 0.72 pN
20
40
60
80
100
simulation #
5
Appendix B
Figure 2: Plot of theoretical values of pN s
0.71 pN 0.67 2 0.69
200
400 N
600
800
1000
45 pN = 0.721?? 40 35 30 11 200 400 N 600 800 1000
60 55 50 45 40 35 30 41 200 400 N 600 800 1000
pN = 0.7213??
6
Appendix C
TOSS
# There will be n.trials trials in each of n.sims simulations for each value # of N. For each simulation, proportion of 1s over n.trials trials is # calculated for the estimation of pN. n.trials <- 1000 n.sims <- 100 N <- c(10,50,100,500,1000,2000,5000) nN <- length(N) pNs <- matrix(0,n.sims,nN) for(j in 1:nN) for(i in 1:n.sims){ for(m in 1:n.trials) pNs[i,j] <- pNs[i,j]+toss(N[j]) } pNs <- pNs/n.trials
PNHAT
> cpu <- numeric(100) > for(i in 1:100)cpu[i] <- cpu.time(pN.hat(10000,10000)) > mean(cpu) # mean CPU for using pN.hat (i.e., vectorization) [1] 0.001683168 > expr <- expression({x <- 0;for(i in 1:10000)x <- x+toss(10000);x/01000}) > for(j in 1:100)cpu[j] <- cpu.time(eval(expr)) > mean(cpu) # mean CPU for using explicit loop (i.e., non-vectorization) [1] 1.384059
IMPROVE
n.trials <- 1000 n.sims <- 100 N <- c(10,50,100,500,1000,2000,5000) nN <- length(N) pNs <- matrix(0,n.sims,nN) for(j in 1:nN) for(i in 1:n.sims) pNs[i,j] <- pN.hat(N[j],n.trials)
PLOT
pdf(file=toss.pdf, width=6, height=8.5) par(mfcol=c(4,2)) # column-wise 4 by 2 multiple figures ylim <- range(pNs) # range of the simulated pNs for(i in 1:7) plot(pNs[,i],xlab=simulation #,ylab=bquote(p[N]), ylim=ylim,type=l, main=paste("N =",N[i])) dev.off()
7
Appendix D
N <- 1000; fwid <- c(.02,.04); xlim <- c(0,N) probs <- pN(N) pdf(file="pNGraph.pdf",width=4.5,height=6) par(mfrow=c(3,1),mar=c(5,4,1,3)+.1,xpd=T) y.range <- range(my <- probs[mx <- 41:N]) plot(2:N,probs[-1],xlab=N,ylab=bquote(p[N]), xlim=xlim,xaxt=n,type=l) axis(1,at=c(2,seq(200,N,200))) usr <- par(usr) x <- usr[2] + (usr[2]-usr[1])*fwid rect(x[1],y.range[1],x[2],y.range[2]) axis(4,label=F,tcl=-0.3) plot(11:N,probs[-(1:10)],xlab=N,ylab=bquote(p[N]==0.721*plain("??")), xlim=xlim,axes=F,type=l) axis(1,at=c(11,seq(200,N,200))) axis(2,at=(y<-seq(0.72130,0.72145,0.00005)), label=as.character(round((y-0.721)*105)),las=1) box() rect(x[1],y.range[1],x[2],y.range[2]) axis(4,label=F,tcl=-0.3) plot(mx,my,xlab=N,ylab=bquote(p[N]==0.7213*plain("??")), xlim=xlim,axes=F,type=l) axis(1,at=c(mx[1],seq(200,N,200))) axis(2,at=(y<-seq(0.721330,0.721360,0.000005)), label=as.character(round((y-0.7213)*106)),las=1) box() rect(x[1],y.range[1],x[2],y.range[2]) axis(4,label=F,tcl=-0.3) dev.off()
8