Fermat and his library

This morning I woke up to a delightful tweet from fermatslibrary about sample random uniform numbers and how many it takes, on average, to sum to 1.

If you look at the embedded picture, you can see the math sketched out but of course it’s alwasy more fun to simulate.

Some R Code

How can we generate a random uniform [0,1] number in R? There are many distributions built in and this one is included. Note, I’m setting a seed to ensure if you’re following along at home you will get the same answer. runif is the function which generates random values, by default min is 0 and max is 1 (just what we want).

set.seed(1001)
runif(n = 1,min = 0,max = 1)
## [1] 0.9856888

We can visually ‘verify’ the distribution we’re drawing from with a histogram. While it’s not totally flat, it’s fairly obvious that it’s uniformly sampling values between 0 and 1.

hist(runif(10000), main = "Histogram of 10,000 Random Uniform [0,1] Values",xlab="")

How many are needed?

To determine how many values are needed to sum to 1, we’ll use a while loop. This is a loop that will cycle until a logical condition is satisified. It’s like a for loop except it has the potential to be an infinite loop.

summedValue = 0
steps = 0
while(summedValue < 1){
  summedValue = summedValue + runif(1)
  steps = steps+1
}
cat("Total Summed Value:",summedValue,"\n")
## Total Summed Value: 1.753755
cat("It took ",steps," step",ifelse(steps==1,"","s")," to get there\n",sep="")
## It took 2 steps to get there

This is great! We can see to total value, and how many values were drawn to get to 1. But now we need to do this lots of times.

Functional programming

By wrapping our above code in a function we’ll be able to take advantage of functional programming. Note how the function is mostly a wrapper and all that’s changed is we replaced the cat steps with a return argument (what the function will return as a value).

countUnif = function(){
  summedValue = 0
  steps = 0
  while(summedValue < 1){
    summedValue = summedValue+runif(1)
    steps = steps+1
  }
  return(steps)
}

countUnif()
## [1] 3

Now we can simulate this many, many times. The replicate function makes it very easy to generate 1,000,000 experiments.

results = replicate(1000000,countUnif())
cat("Average values required from simulation:",mean(results),"\n")
## Average values required from simulation: 2.717411
cat("Average values required from math      :",exp(1),"\n")
## Average values required from math      : 2.718282

Wow! That’s really close. But I’d really like to see how the value stabilizes over time. So, I’d like a cumulative average.

cumulativeAverage = cumsum(results)/(1:1000000)

Visualize!

plot(cumulativeAverage,type='l',main='Average Random Uniform [0,1] Values Needed to Sum to 1',ylab='Cumulative Average')
abline(h=exp(1),col='red')

Zoom and enhance! Here we can see that we’re not necessarily stabilizing,but the movement is all in a very small range.

plot(cumulativeAverage,type='l',main='Average Random Uniform [0,1] Values Needed to Sum to 1',ylab='Cumulative Average',ylim=c(2.715,2.72))
abline(h=exp(1),col='red')

Conclusion

A fun and easy simluation to do in R, allowing functional programming practice, to demonstrate a neat math fundamental.