# How many random numbers does it take?

## 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.

Pick a uniformly random number in [0,1] and repeat until the sum of the numbers picked is >1. You'll on average pick e≈2.718… numbers! pic.twitter.com/8ak1hYENCi

— Fermat's Library (@fermatslibrary) October 28, 2017

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.