Chapter 2 Exercises
The "how many" trick
Solutions

1. Imagine simulating a series of random numbers from a uniform distribution on the interval (0,1) and adding them up. Adapt the code in Example 2.5.1 to calculate the expected number of terms required to make the sum exceed one. [Hint: this will be one more than the largest number for which the sum < 1.] Does the resulting expected value look familiar?

model {
for (i in 1:max) {
x[i] ~ dunif(0,1)
}   
cum[1] <- x[1]
for (i in 2:max) {
cum[i] <- cum[i-1] + x[i]
}
for (i in 1:max) {
cum.step[i] <- i*step(1 - cum[i])
}
num <- ranked(cum.step[], max) + 1
}

list(max = 20)

   node   mean   sd   MC error   2.5%   median   97.5%   start   sample
   num   2.715   0.871   0.003782   2.0   2.0   5.0   1   50000

The closed form solution is
e (Derbyshire, J. Prime Obsession: Bernhard Riemann and the Greatest Unsolved Problem in Mathematics. New York: Penguin, 2004, pp. 366--367). This remarkable number has a habit of turning up in unexpected places.

2. Add one line to your code to allow calculation of the probability that n terms will sum to less than one. What is the probability for n=5?

We can add the line in red below and the mean values of lessthan1[i] will give the required probabilities:

model {
for (i in 1:max) {
x[i] ~ dunif(0,1)
}   
cum[1] <- x[1]
for (i in 2:max) {
cum[i] <- cum[i-1] + x[i]
}
for (i in 1:max) {
cum.step[i] <- i*step(1 - cum[i])
lessthan1[i] <- step(1 - cum[i])
}
num <- ranked(cum.step[], max) + 1
}

list(max = 20)

   node   mean   sd   MC error   2.5%   median   97.5%   start   sample
   lessthan1[1]   1.0   0.0   4.472E-13   1.0   1.0   1.0   1   50000
   lessthan1[2]   0.4993   0.5   0.002173   0.0   0.0   1.0   1   50000
   lessthan1[3]   0.1653   0.3715   0.001716   0.0   0.0   1.0   1   50000
   lessthan1[4]   0.04088   0.198   8.464E-4   0.0   0.0   1.0   1   50000
   lessthan1[5]   0.0079   0.08853   3.752E-4   0.0   0.0   0.0   1   50000
   lessthan1[6]   0.0013   0.03603   1.63E-4   0.0   0.0   0.0   1   50000
   lessthan1[7]   2.2E-4   0.01483   6.478E-5   0.0   0.0   0.0   1   50000
   lessthan1[8]   6.0E-5   0.007746   3.44E-5   0.0   0.0   0.0   1   50000
   lessthan1[9]   0.0   0.0   4.472E-13   0.0   0.0   0.0   1   50000
   lessthan1[10]   0.0   0.0   4.472E-13   0.0   0.0   0.0   1   50000
   lessthan1[11]   0.0   0.0   4.472E-13   0.0   0.0   0.0   1   50000
   lessthan1[12]   0.0   0.0   4.472E-13   0.0   0.0   0.0   1   50000
   lessthan1[13]   0.0   0.0   4.472E-13   0.0   0.0   0.0   1   50000
   lessthan1[14]   0.0   0.0   4.472E-13   0.0   0.0   0.0   1   50000
   lessthan1[15]   0.0   0.0   4.472E-13   0.0   0.0   0.0   1   50000
   lessthan1[16]   0.0   0.0   4.472E-13   0.0   0.0   0.0   1   50000
   lessthan1[17]   0.0   0.0   4.472E-13   0.0   0.0   0.0   1   50000
   lessthan1[18]   0.0   0.0   4.472E-13   0.0   0.0   0.0   1   50000
   lessthan1[19]   0.0   0.0   4.472E-13   0.0   0.0   0.0   1   50000
   lessthan1[20]   0.0   0.0   4.472E-13   0.0   0.0   0.0   1   50000
   
   So the probability that 5 terms will sum to less than one is under 1%