Chapter 2 Exercises
The “how many” trick
<span style="color:blue"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%