Sums from a uniform random distribution

Recently, I came across some pages about a problem that I simulated in Python years ago, but never solved analytically, although a reader left a hint in the Comments.

Here we start by taking a look at what this author calls "random number blackjack" or the problem of the sum of uniformly distributed random numbers.

For the rest of the post, all the random numbers that we will talk about are given by a random number generator which draws from uniformly distributed real numbers in the range [0,1), i.e. greater than or equal to zero, and less than one.  

We also consider the sum of one or more independent values from such a distribution, which forms a different distribution.

the game

On each turn of this unusual game, you draw a random number.  You receive $100 for each turn you play, including the last.  Suppose you drew 0.45, 0.35, 0.25, on the third turn you would go bust, since the sum of the numbers drawn would be 1.05 > 1.0.  The payout would be $300.

And the question is, what is your expected winnings if you play the game many times?  Is it worth it to play (accepting some risk of bad luck) if the initial buy-in requires $250?

initial thoughts

First of all, since the generator always gives a number less than 1.0, you always play at least two turns.  The expected value (mean) of the first number is 0.5.  For a continuous distribution, the mean has a technical definition, but it is always on the axis of symmetry of the distribution, which is obviously 0.5 here.

The expected value of the sum of several random numbers is the sum of their expectations.

So, for example, E[3] = 1.5.  Most of the time we will quit after two or three turns, but there will occasionally be an extended run of smaller numbers and a corresponding increase in the number of turns.

[spoiler alert, you may want to look at random number blackjack to try to work out the answer]

What famous number do we know that lies between two and three?  This is so trivial to simulate in Python I won't even post an example.

On one run I got 2.718531.

So it looks like the result is equal to e.  (The internet says that adding more rounds doesn't help the accuracy because of limitations in the random number generator).

I came across this problem in a slightly different form in the terrific introductory calculus book, Calculus Made Easy, which was originally written by Sylvanus P. Thompson (and available as a Project Gutenberg download).  

It was added by Martin Gardner when he edited the classic work (1998 ed., p. 153) and is simply a note about how e shows up everywhere.

But the problem is at least as old as Feller's classic text Probability (which I don't have, unfortunately).

related problem

Tim Black solves a related problem (here).  

Recall that for a standard (fair) die or dice, the expected value of a throw is the sum of each value times the probability that it will occur.  For a single die, the average is (1 + ... 6).1/6 = 21/6 or 3.5.  

For two dice, we have

1(2) + 2(3) + ... + 6(7) + ... + 2(11) + 1(12)

The distribution is no longer uniform, but it is still symmetric around the value of 7, which is the mean.  The expected values from adding two random draws from a uniform distribution are observed to add, but the resulting distribution is no longer uniform. 

Suppose we know a probability distribution for the sum of n random numbers, for some value of n, and then calculate the probability that the sum is greater than or equal to 1.  We can then obtain the expected value over a number of trials as that value n times the probability we calculated.

The probability distribution for the sum of two random numbers from the generator has a mean of 1.0, so the probability of exceeding 1.0 is 0.5.  That event has a weight of 2, so the contribution to the total expected value for a number of trials is 2(0.5) = 1.  So, in the same way as we did for the dice, we have that P(2) = 0.5.  

We also have that P(1) = 1.  That's another 1 to add to the expected value overall.

So now, what is the probability distribution for the sum of three random numbers?  That gets a little trickier.  The difficulty is that the probability distribution changes as n changes.  Eventually, it becomes normal, but how different is it for small n like 3, 4, or 5?

Here is where our analyst has a great idea.

Imagine that we change the game slightly.  We still have $100 as the payout at each stage.

From a stream of random numbers, as the numbers are drawn we write into another stream the sum at each stage.  So in the example above we would get 0, 0.45, 0.80,  and then at the third draw the sum is 1.05.  Rather than write the last value, subtract 1 first, then write that down and keep going.  

Notice that over this stretch we have a valid game, a sequence of increasing values followed by one that must be smaller than the last.

The values are

0 0.45 0.80 0.05

The values are in ascending order until the last one, which might be anything smaller than 0.80.  This must be true for a single round from the game according to the rules we have set up.

Since there are n! ways of arranging n values, and only one of those arrangements has the numbers in strictly ascending order, the probability of the event (for a random uniform distribution) is 1/n!.  In other words, starting at the beginning of a stream of random numbers

the probability of obtaining a result of 1 is 1.

the probability of obtaining a result of 2 is 1/2!.

the probability of obtaining a result of 3 is 1/3!.

Multiplying the value of each result by its probability we end up with 1 + 1 + 1/2!.  The expected value of a series of games is the sum of all the possibilities.  

E = sum 1/n!

This is just the infinite series for e.

simulations

I wrote two simulations to show results relevant to this problem.  The first one shows the distribution of sums of n = 1, 2, 3 or 4 random numbers.  As you can see from the figure

even 3 at a time, the sums look pretty close to a normal distribution.  The Central Limit Theorem says that they will tend to normal, and there is a bunch of theory that I don't understand that says if the draws are from a uniform distribution then the convergence is very rapid.

I got curious about this alternate game, so I wrote a simulation which shows that the sum of random numbers, when computed as above, by discarding the whole numbers from the result, appears to be still random uniform.  (gist here).  The original data and the summed series are plotted in the same hisotgram with transparency 0.5.  The new data is random uniform or close to it.

I don't know what the theoretical explanation for this is.  However, if it's true, then rather than do the sums, we can just draw from the random uniform distribution, and tally up the runs where all the values are increasing, until the last one  If we do the bookkeeping correctly, we get e as the result.

That means the original problem has the same as the alternative one.

serious analysis

I have reworked what is on the Mathworld page as follows:

That's where I am so far.  There's plenty more to investigate.  

The sum of two random numbers from a uniform distribution has a distribution that is given by convolution of the individual distributions.  But then each distribution for n > 2 is formed by another convolution.  Ultimately, the distributions tend to the normal.  

I don't see how you get to something as simple as 1 - 1/n! from that, although Tim Black gave us a different path above, which is why I wrote this post.

[Update:  I've been dense.  The "different path" is in fact the means by which the integral is evaluated.  It is not done by writing some complex expression and then seeking the antiderivative and evaluating it.  Instead, we know that the value for the cumulative distribution function at the upper bound must be 1, and at the lower bound it must be 1/n!. ]

There is a suggestion that this sort of thing is done more easily with generating or characteristic functions.  

Probably the first, simple thing would be to run the simulation using random numbers and not bother with the sum part, as we also looked at here.  [Update:  the result is as I suspected.  See gist.  If we simply find the length of runs in a stream of random numbers from a uniform distribution, where they are in increasing order, and then find the mean of those lengths, the result is e.]

Poisson and airline reservations

In the problems for Chapter 5.1, Grinstead and Snell have the following:

An airline finds that 4 percent of the passengers that make reservations on a particular flight will not show up. Consequently, their policy is to sell 100 reserved seats on a plane that has only 98 seats. Find the probability that every person who shows up for the flight will find a seat available.

---

Attempt #1:
Each seat that is sold is a Bernoulli trial with Prob(no-show) = 0.04.

N = 100, so the expected number of no-shows for a flight = N p = 4.
The variance is N p (1-p) = 3.84, sd = 1.95

If the number of no-shows exceeds ≈ 1 sd, there won't be a seat for someone.
I expect this to happen about 0.5*e^-1 of the time = 0.18.
The desired probability is 1 minus this.

The problem with this is that no-shows are not normally distributed!

---

Attempt #2:

If we use the Poisson approximation we have

P = λk / k! e-λ lambda = 4 P(0 no-show) = e-4 = 0.018 P(1 no-show) = 4 * 0.018 = 0.073 P(nobody gets bumped) = 1 - 0.073 - 0.018 = 1 - 0.091 = 0.919

---

Attempt #3 is a simulation (Python code below.
Note the number of seats simulated is 10⁸ !).
The program prints the probability that no one gets bumped, and the st dev for the simulated values.

$ python problem28.py 0.087023 0.00877316767194

import numpy as np import random def overbooked(): rL = list() for i in range(100): if random.random() > 0.96: pass else: rL.append(1) return sum(rL) > 98 N = 1000 L = list() for i in range(1000): S = sum([overbooked() for i in range(N)]) L.append(S*1.0/N) A = np.array(L) print np.mean(A), np.std(A)

Plotting the normal distribution with Python

It is nice to be able to add a plot of the normal distribution on top of another plot, say a histogram of your data. I've done it before from R (here) using code like this (which assumes we have some data in an array M):

plot(function(x) dnorm(x,0,sd(M)), min(M),max(M),lwd=3, add=T,col='red')

I wanted to find out how to do this using numpy and matplotlib. It turns out that while R has these functions built-in, numpy doesn't seem to have them. On this page showing an example exactly like what I want to do, the SciPy docs implement the pdf for the normal distribution as if it had come straight out of wikipedia. So, we'll do the same!

Here is a first shot at it, just the normal pdf, without any sample data. If you want "curves", (or what look like curves but are actually very closely spaced "points"), just make the variable dx smaller. The plot code is similar to what we did the other day, with the Poisson distribution.

import math import numpy as np import matplotlib.pyplot as plt def normal(mu,sigma): def f(x): z = 1.0*(x-mu)/sigma e = math.e**(-0.5*z**2) C = math.sqrt(2*math.pi)*sigma return 1.0*e/C return f X = 2 dx = 0.1 R = np.arange(-X,X+dx,dx) L = list() sdL = (0.5,1,2,3) for sd in sdL: f = normal(mu=0,sigma=sd) L.append([f(x) for x in R]) colors = ['r','b','purple'] for c,P in zip(colors,L): plt.plot(R,P,zorder=1,color='0.2',lw=1.5) plt.scatter(R,P,zorder=2,s=50,color=c) ax = plt.axes() ax.set_xlim(-2.1,2.1) #ax.set_ylim(-0.01,0.5) plt.savefig('example.png')

Checking an exponential approximation

In the previous post about the Poisson approximation to the binomial distribution, we said that "since n is large and p is small":

(1-p)n ≈ e-np

I wanted to take a look at the accuracy of the approximation. In the plot at the top, p varies from 0.005 to 0.03 as shown on the x-axis, for values of n in the series 10,30,100,300 (red, blue, purple, gray). The first thing to observe is that, if p is small, the approximation is very good. The error is < 1% for all values of p < 0.25 if n is equal to 10 or 30. However, it is not necessary that n be large. In fact, the error is much worse for n = 100 or 300. Still, for 100 trials with p = 0.01, the error is < 1%.

R code:

f <- function(p,n) { left = (1-p)**n right = exp(-p*n) (right - left)/left } plot(1,type='n', xlim=c(0,0.03),ylim=c(0,0.10), xlab='p',ylab='error') p = seq(0.005,0.04,by=0.005) points(p,f(p,n=10),pch=16,cex=3,col='red') points(p,f(p,n=30),pch=16,cex=3,col='dodgerblue') points(p,f(p,n=100),pch=16,cex=3,col='purple') points(p,f(p,n=300),pch=16,cex=3,col='gray60')

Poisson approximation to the binomial

I came across a nice introduction to probability that starts from sets (here). It contains a small section extending the binomial distribution to the Poisson, which explains the simplifying assumptions in an intuitive way. I'm going to do something similar here, except that I'm following the derivation in my favorite bioinformatics textbook, Higgs & Attwood.

This will complement a previous post, which explored the approximation of the binomial as the normal distribution.

Say we're carrying out a series of Bernoulli trials, flipping a (possibly unfair) coin, or just throwing a loaded die where the possible results are divided into two classes: success and failure. If the probability of success is the same on any trial, denoted p (and that of failure is q = 1-p), then the probability of a particular sequence like SFSSSFSFSS containing k successes and n-k failures in a total of n trials is

p*q*p*p*p*q*p*q*p*p

We can gather all the like terms together to give

pk qn-k = pk (1-p)n-k

We can obtain the total probability for all trials having the same number of successes. Multiply the previous result by the number of combinations, the number of ways of picking k successes out of n trials, which I will symbolize as _nC_k for convenience. This is called "n choose k."

nCk = n! / (n-k)! k!

The complete expression for the binomial distribution is:

P(k) = [n! / (n-k)! k!] pk (1-p)n-k

The Poisson approximation applies to cases where p is small and n is large enough that p/n, the mean number of successes, stays reasonable (like 1 or 2), even as we make n large and p shrinks closer and closer to zero. We will (eventually, below) use the symbol λ for that ratio.

But start by looking at the term due to combinations. By our assumptions n is large and p is small, so that k is reasonable for the cases we are interested in (near the mean and smaller). Suppose k = 2. Then:

n!/(n-k)! = n * n-1 * n-2 * n-3... ≈ n2 ------------ n-2 * n-3...

So, generalizing to k, and also remembering the factor of 1/k! from the original expression, the complete term for the combinations is approximately

nCk ≈ nk / k!

The second simplification comes from the probability term for failures. Separating the exponents:

(1-p)(n-k) = (1-p)n * (1-p)-k

Since n is large and p is small:

(1-p)n ≈ e-np (1-p)-k ≈ 1

(Note: the first approximation will need some justification, that I'll defer to a later post. Also, it's not clear to me why we can ignore p in the second case but not in the first).

From these two approximations we obtain:

P(k) = nk / k! * pk * e-np

Substitute λ = np to reveal the familiar Poisson distribution:

P(k) = λk / k! * e-λ

One interesting property of the Poisson distribution is that the variance is equal to the mean, so that the curves spread out as λ increases.

R code:

poisson <- function(lambda) { p <- function(k) { e = exp(-lambda) f = factorial(k) g = lambda**k g*e/f } } x = 0:20 plot(x,type='n',xlim=c(0,20),ylim=c(0,0.4), ylab='P(x)',xlab='x') y = poisson(1)(x) lines(x,y,col='black',lwd=1) points(x,y,pch=16,cex=2,col='red') y = poisson(4)(x) lines(x,y,col='black',lwd=1) points(x,y,pch=16,cex=2,col='dodgerblue') y = poisson(10)(x) lines(x,y,col='black',lwd=1) points(x,y,pch=16,cex=2,col='purple')

Exponential density 3

In example 4.20 of Grinstead and Snell there is a nice conjunction of Bayes theorem and use of the exponential density. Recall that the exponential pdf is:



To find the probability that X happens (hard-drive failure, radioactive decay) within a certain time period, we integrate the pdf over the interval. For example, the probability that the failure happens after a particular time t is:



(Naturally, since the cdf(t) is 1 minus this value).

Here is one form of Bayes theorem:



Now, consider two events E and F defined as follows:

E is the event that failure happens after time r F is the event that failure happens after time r + s

Note that P(F and E) = P(F) because F is totally contained within E.

Then:



The probability of failure after time r + s, when we know that failure occurs after time r, does not depend on r at all but is only a function of s. This is the memoryless property of the exponential function, alluded to in a previous post.

Exponential density 2

Think of e as a function (rather than the irrational number 2.71828...).

In R:

e <- function(x) { 2.71828**(x) } plot(e,0,5,lwd=10,col='gray70') plot(exp,0,5,lwd=2,col='red',add=T)

Then, e can be defined as the function whose derivative is itself. To see this:

plot(exp,0,5,lwd=5,xlim=c(0,5),ylim=c(0,75)) par(lwd=2) colors=c('blue','darkgreen','red','magenta') f <- function(x) { y=exp(x) points(x,y,col=colors[x],pch=16,cex=2) lines(c(x-1,x+1),c(y-y,y+y), col=colors[x]) } for (i in 1:4) f(i)

What did we do? First, we plotted the exponential function exp between 0 and 5 (heavy black line). Then, for each x in the series 1:4, we calculated y=exp(x) and plotted the point x,y in color. At the same time, we constructed a line with slope (y+y)/2 = y, passing through x,y. It is clear that the slope of the curve is equal to the value of the function at that point.

Another fun way to see this is to look at the infinite series for e^x:



Can you see that the derivative with respect to x of this series is identical to the series itself?

Using simple calculus (the chain rule), we can also show that



Looking at the probability density function (pdf) for the exponential distribution,



and cumulative distribution function, we confirm that the pdf is the derivative of the cdf, as it should be. The cdf is:

Exponential density

I've been reading (and re-reading) An Introduction to Probability by Grinstead and Snell. It is a wonderful book, available from here as a pdf. It goes slowly, has lots of explanation and many problems, as well as interesting historical perspective. I like it so much that I bought a hard copy.

I'm trying to understand the exponential density better. In example 2.17 of the book they pose the problem of modeling the time-to-breakdown of a hard drive by the exponential density:



If the average time-to-breakdown is 30 months, and we have already run the computer for 15 months with no breakdown, what is the current expected time-to-breakdown?

We use R to explore the question. The R function rexp gives random samples from the exponential density with a rate parameter r (the inverse of lambda above). Think of each of these as a possible lifetime for our drive. Since we know that the lifetime exceeds 15 months, filter the vector x for values > 15 and save in y.

r = 1/30 x = rexp(100000,rate=r) sel = x > 15 y = x[sel]

Plot histograms of the density (not counts, which R calls freq).

hist(y,breaks=100,xlim=c(0,150),freq=F) hist(x,breaks=100,col='gray70',freq=F,add=T)

It is clear that the y distribution is the same as x, just shifted over by 15.



We confirm this by looking at summary statistics (adjusting y first by subtracting 15). We see they are essentially identical:

> summary(x) Min. 1st Qu. Median 4.915e-05 8.580e+00 2.075e+01 Mean 3rd Qu. Max. 3.001e+01 4.170e+01 3.780e+02 > summary(y-15) Min. 1st Qu. Median 1.309e-05 8.627e+00 2.087e+01 Mean 3rd Qu. Max. 3.004e+01 4.185e+01 3.630e+02 freq=F,add=T)