The statistics of a Galton Board

V. Hunter Adams

In [244]:
from IPython.core.display import HTML

HTML("""
<style>
.output_png {
    display: table-cell;
    text-align: center;
    vertical-align: middle;
}
</style>
""")

HTML('''<script>
code_show=true; 
function code_toggle() {
 if (code_show){
 $('div.input').hide();
 } else {
 $('div.input').show();
 }
 code_show = !code_show
} 
$( document ).ready(code_toggle);
</script>
<form action="javascript:code_toggle()"><input type="submit" value="Click here to toggle on/off the raw code."></form>''')
Out[244]:
In [247]:
import numpy
import matplotlib.pyplot as plt

def plotBinomial(num_rows, overlay=False, showplot=True):
        
    # Compute binomial probabilities
    p = 0.5**(float(num_rows))
    probabilities = numpy.zeros(num_rows+1)
    for i in range(num_rows+1):
        binomial_coefficient = ((numpy.math.factorial(num_rows))/
                                (numpy.math.factorial(i)*numpy.math.factorial(num_rows-i)))
        probabilities[i] = p*binomial_coefficient
        
    # Compute Gaussian probabilities
    gaussian_xvals = numpy.linspace(-1, num_rows+1, 100)
    gaussian_sigma = numpy.sqrt(num_rows/4)
    gaussian_yvals = ((1./(gaussian_sigma*numpy.sqrt(2*numpy.pi)))*
                      numpy.exp(-0.5*(((gaussian_xvals-(num_rows/2.))**2.)/(gaussian_sigma**2.))))
    
    # Should we plot?
    if showplot:
        
        # Setup and beautify
        fig, ax = plt.subplots(figsize=(21, 5))
        ax.tick_params(axis='both', which='major', labelsize=15)
        ax.set_xlabel('Bin number', fontsize=18)
        ax.set_ylabel('Probability mass', fontsize=18)
        ax.set_title(str(num_rows)+'-row Galton Board bin probabilities', fontsize=24)
        
        # Include the Gaussian in the plot?
        if overlay:
            ax.plot(gaussian_xvals, gaussian_yvals, 'r-')
            
        # Plot the binomials
        ax.bar(numpy.arange(num_rows+1), probabilities,
               tick_label=numpy.arange(num_rows+1) if num_rows<25 else None,
               label='Binomial distribution')
        
        # If we've included the Gaussian, add a legend
        if overlay:
            ax.legend(['Gaussian', 'Binomial'], fontsize=18)
        plt.show()
    else:
        return [probabilities, gaussian_yvals]

Setting up the problem

In an idealized Galton Board, like the one shown below, a ball bounces against a single peg in each row of the board. When the ball strikes the peg, it has a 50 percent probability of bouncing to the left, and a 50 percent probability of bouncing to the right. In this idealization, the outcome of the bounce in each row is completely independent from the bounces in previous rows. An array of bins catch the balls as the fall through the bottom row of the board.

missing
Fig. 1: Galton Board overview

We can label each of these bins. Starting from the left, we label the bins with increasing integers starting at 0. The number of bins will be equal to one more than the number of rows in our Galton Board. The ball falls onto the top-center peg, at a position which is half the number of rows in the board. For the 5-row Galton Board shown above, for instance, the ball falls onto the board at position 2.5. Each peg strike has a 50 percent probability of increasing its position by $\frac{1}{2}$ (ball falls to the right) and a 50 percent probability of decreasing its position by $\frac{1}{2}$ (ball falls to the left). The bin in which it ultimately lands represents the sum of $N$ such random trials, where $N$ is the number of rows in the board.

Note that we've chosen the $\frac{1}{2}$ step size so that the ball always falls out of the board at an integer position, which corresponds to our labeled bins. We might have chosen different conventions for our bin labels and step sizes. The chosen one, I think, makes the final distributions the most intuitive.

Distribution for a one-row Galton Board

A one-row Galton Board isn't much of Galton Board, its just a peg against which a ball bounces. But it's the simplest case and offers a good starting point! After bouncing against the peg, the ball has a 50 percent probability of decreasing its position by $\frac{1}{2}$ (bouncing left) and a 50 percent probability of increasing its position by $\frac{1}{2}$ (bouncing right). This situation describes a Bernoulli trial - a random experiment with exactly two possible outcomes. To keep the math general for the moment, we'll say that the ball bounces to the left with a probability $p$. Thus, the position of the ball after the first bounce, a random variable that we'll call $X_1$, is given by:

\begin{align} X_1 \sim \text{Bernoilli}(p) \end{align}

If we set $p=0.5$, we can generate the (boring and predictable) distribution of balls in the two bins underneath our single peg. For our one-row Galton Board, we assume that the ball is dropped from the x-position 0.5 (half the number of rows). Thus, the ball will leave the Board at position 0 or 1, and has equal probability of landing in either.

In [248]:
plotBinomial(1)

Distribution for a two-row Galton Board

The distribution for the sum of two random variables is given by the convolution of the probability density function for each random variable (or, in the case of our discrete random variables, the convolution of their probability mass functions). Let us consider the sum of two Bernoulli random trials.

\begin{align} X_1 &\sim \text{Bernoilli}(p)\\ X_2 &\sim \text{Bernoilli}(p)\\ Y &= X_1 + X_2 + \mu \end{align}

For our idealized Galton board, $X_1$ and $X_2$ each have a 50 percent probability of assuming $\frac{1}{2}$ or -$\frac{1}{2}$, but we'll keep this general and assume an outcome of -$\frac{1}{2}$ with probability $p$. $\mu$, in the above expression, is the position from which the ball is dropped (half the number of rows) and is not a random number. This means that the possible values which $Y$ could assume include 0, 1, or 2. We would like to find the probability of each of these outcomes, and to make a nice histogram of those probabilities. This histogram will be the probability mass function for the random variable $Y$.

Let's start by computing each of these probabilities separately, then we'll generalize. For a two-row Galton Board, $\mu=1$. To compute the probability that $Y=0$, we can sum over all possible values for $X_1$ to get:

\begin{align} P(Y=0) &= P(X_1+X_2=-1)\\ &= P(X_1=-\frac{1}{2},X_2=-\frac{1}{2}) + P(X_1=\frac{1}{2},X_2=-\frac{3}{2})\\ &= P(X_1=-\frac{1}{2})P(X_2=-\frac{1}{2}) + P(X_1=\frac{1}{2})P(X_2=-\frac{3}{2})\\ &= (1-p) \cdot (1-p) + p \cdot 0\\ &= (1-p)^2 \end{align}

(Note that if the required value for $X_2$ in order to achieve $Y=-1$ is impossible, we assign it a probability 0). If we wanted to compute the probability that $Y=0$, we can similarly sum over all possible values for $X_1$:

\begin{align} P(Y=1) &= P(X_1 + X_2 = 0)\\ &= P(X_1=-\frac{1}{2}, X_2=\frac{1}{2}) + P(X_1=\frac{1}{2},X_2=-\frac{1}{2})\\ &= P(X_1=-\frac{1}{2})P(X_2=\frac{1}{2}) + P(X_1=\frac{1}{2})P(X_2=-\frac{1}{2})\\ &= (1-p) \cdot p + p \cdot (1-p)\\ &= 2p(1-p) \end{align}

And finally, if we wanted the probability that $Y=1$:

\begin{align} P(Y=2) &= P(X_1+X_2=1)\\ &= P(X_1=-\frac{1}{2},X_2=\frac{3}{2}) + P(X_1=\frac{1}{2},X_2=\frac{1}{2})\\ &= P(X_1=-\frac{1}{2})P(X_2=\frac{3}{2}) + P(X_1=\frac{1}{2})P(X_2=\frac{1}{2})\\ &= (1-p) \cdot 0 + p \cdot p\\ &= p^2 \end{align}

Which looks like the following:

In [249]:
plotBinomial(2)

Distribution for a three-row Galton Board

Let's add one more random variable, so that we can see a pattern emerge. Suppose now that $Y = X_1 + X_2 + X_3 + \mu$. For a 3-row Galton Board, $\mu=1.5$ and $Y$ can assume values 0, 1, 2, and 3. Let us compute the probabilities of each of these values:

Probability that $Y=0$:

\begin{align} P(Y=0) &= P(X_1+X_2+X_3 = -1.5)\\ &= P(X_3=-\frac{1}{2}, X_1+X_2=-1) + P(X_3=\frac{1}{2}, X_1+X_2=-2)\\ &= P(X_3=-\frac{1}{2})\cdot P(X_1+X_2=-1) + P(X_3=\frac{1}{2})\cdot P(X_1+X_2=-2)\\ &= (1-p) \cdot (1-p)^2 + p \cdot 0\\ &= (1-p)^3 \end{align}

Probability that $Y=1$: \begin{align} P(Y=1) &= P(X_1+X_2+X_3=-0.5)\\ &= P(X_3=-\frac{1}{2}, X_1+X_2=0) + P(X_3=\frac{1}{2}, X_1+X_2=-1)\\ &= P(X_3=-\frac{1}{2})\cdot P(X_1+X_2=0) + P(X_3=\frac{1}{2})\cdot P(X_1+X_2)=-1\\ &= (1-p) \cdot 2p(1-p) + p(1-p)^2\\ &= 2p(1-p)^2 + p(1-p)^2\\ &= 3p(1-p)^2 \end{align}

Probability that $Y=2$: \begin{align} P(Y=2) &= P(X_1+X_2+X_3=0.5)\\ &= P(X_3=-\frac{1}{2}, X_1+X_2=1) + P(X_3=\frac{1}{2}, X_1+X_2=0)\\ &= P(X_3=-\frac{1}{2})\cdot P(X_1+X_2=1) + P(X_3=\frac{1}{2})\cdot P(X_1+X_2)=0\\ &= (1-p) \cdot p^2 + p\cdot 2p(1-p)\\ &= p^2(1-p) + 2p^2(1-p)\\ &= 3p^2(1-p) \end{align}

Probability that $Y=3$: \begin{align} P(Y=3) &= P(X_1+X_2+X_3=1.5)\\ &= P(X_3=-\frac{1}{2}, X_1+X_2=2) + P(X_3=\frac{1}{2}, X_1+X_2=1)\\ &= P(X_3=-\frac{1}{2})\cdot P(X_1+X_2=2) + P(X_3=\frac{1}{2})\cdot P(X_1+X_2=1)\\ &= (1-p) \cdot 0 + p \cdot p^2\\ &= p^3 \end{align}

Which looks like the following:

In [250]:
plotBinomial(3)

Distribution for an n-row Galton Board

At this point, we can see a pattern emerging. In the expression below, $k$ is the position of the bin as labeled in the figure at the top of the webpage, and $n$ is the total number of rows in the Galton board.

\begin{align} P(Y=k) &= \begin{pmatrix}n\\k\end{pmatrix}p^{k}(1-p)^{n-k} \end{align}

Indeed, we can prove by induction that this is the case. As our base case, we've shown that the probability of landing in bin $k$ for a one-row Galton Board is distributed as $B(1,p)$. (We've actually also shown that the probability of landing in bin $k$ for a two-row Galton Board is distributed as $B(2,p)$, and that the probability of landing in bin $k$ for a three-row Galton Board is distributed as $B(3,p)$). With our base case, we make the induction hypothesis shown above.

For our induction step, we consider adding one more row to our Galton Board. In effect, this adds one more independent Bernoilli trial with probability of bouncing left $p$, and probability of bouncing right $(1-p)$. Let's call this additional Bernoilli trial $Z$, and define a new random variable $W = Y+Z$. By our hypothesis, $Y \sim B(n,p)$.

\begin{align} P(W=k) &= P\left(Y=k+\frac{1}{2}\right)p + P\left(Y=k-\frac{1}{2}\right)(1-p) \end{align}

We've shown that the distribution for the $n+1$-row Galton Board depends on that of the $n$-row Galton Board in a way that mirrors the Binomial distribution. Thus, each additional row that we add to the board preserves this distribution. This proves that the Binomial distribution describes the sum of Bernoilli trials. And in the case of our idealized Galton board $p=0.5$, and so this simplifies even further.

\begin{align} P(Y=k) &= \begin{pmatrix}n\\k\end{pmatrix}\frac{1}{2^n} \end{align}

As we increase the number of rows in our Galton Board, here's what happens to our distribution. Do you see an old friend emerging?

In [255]:
import matplotlib.pyplot as plt
from matplotlib import animation, rc

def barlist(n): 
    probs = plotBinomial(n, False, False)[0]
    probabilities = numpy.zeros(51)
    for i in range(len(probs)):
        probabilities[i] = probs[i]
    return probabilities

# fig=plt.figure(figsize=(20, 5))
fig, ax = plt.subplots(figsize=(21, 5))
ax.tick_params(axis='both', which='major', labelsize=15)
ax.set_xlabel('Bin number', fontsize=18)
ax.set_ylabel('Probability mass', fontsize=18)
ax.set_title(str(1)+'-row Galton Board bin probabilities', fontsize=24)

# Compute Gaussian probabilities
# gaussian_xvals = numpy.linspace(-1, 1+2, 100)
# gaussian_sigma = numpy.sqrt(1/4.)
# gaussian_yvals = ((1./(gaussian_sigma*numpy.sqrt(2*numpy.pi)))*
#                   numpy.exp(-0.5*(((gaussian_xvals-(num_rows/2.))**2.)/(gaussian_sigma**2.))))

# line, = ax.plot(gaussian_xvals, gaussian_yvals, 'r-')

n =50 #Number of frames
x = numpy.arange(51)
barcollection = plt.bar(x,barlist(1))

def animate(i):
    y=barlist(i+1)
    ax.set_title(str(i+1)+'-row Galton Board bin probabilities', fontsize=24)
    gaussian_xvals = numpy.linspace(-1, i+2, 100)
    gaussian_sigma = numpy.sqrt((i+1)/4.)
    gaussian_yvals = ((1./(gaussian_sigma*numpy.sqrt(2*numpy.pi)))*
                      numpy.exp(-0.5*(((gaussian_xvals-((i+1)/2.))**2.)/(gaussian_sigma**2.))))
#     line.set_data([gaussian_xvals, gaussian_yvals])
    for i, b in enumerate(barcollection):
        b.set_height(y[i])

anim=animation.FuncAnimation(fig,animate,repeat=False,blit=False,frames=n,
                             interval=100)
plt.close()
rc('animation', html='jshtml')
In [256]:
anim
Out[256]: