Coin flip simulator

The following code will allow you to simulate coin flipping and visualize the effect of the law of large numbers (LLN). To aid in visualization, we decide that for every heads we record a one, otherwise a zero. Thus the sum of n coins is between n and zero. To visualize the LLN, we take the arithmetic average of all flips so far with every new flip. The law of large numbers states that this average will converge to the probabiliy that we flip a heads.

\( P[X = 1] = 1 - P[X = 0] = p \; \Rightarrow \; n^{-1}\sum_{i=1}^n X_i \rightarrow_{n \rightarrow \infty} E[X] = p \;\; (LLN) \)

To actually simulate the arithmetic average, we will first simulate all coin flips and then calculate partial sums. For a list of summable values \( A = (a_i)_{i=1}^n \), its partial sum list \( B=(b_i)_{i=1}^n \) is the a list where the \( i^{th} \) entry is the sum of all entries \( a_j \in A : j \leq i \). Succinctly,

\( B = (b_i)_{i=1}^n = \left( \sum_{j=1}^i a_j \right)_{i=1}^n \)

The code we'll use to implement our partial sums follows:

partSums = function(x) {
    n = length(x)
    m = matrix(rep(1, n * n), nrow = n, ncol = n)
    return((m - upper.tri(m)) %*% x)
}

partSums works by utilizing the fact that a lower triangular matrix of ones multiplied by a vector results in its partial sums. This is an \( O(n^2) \) algorithm, whereas an obvious \( O(n) \) algorithm exists. So why do I use it? Because it's an interesting reformulation.

\[ B = (b_i)_{i=1}^n = \left( \begin{array}{cccc} 1 & 0 & \cdots & 0 \\ 1 & 1 & \cdots & 0 \\ \vdots & \vdots & \ddots & 0 \\ 1 & 1 & \cdots & 1 \end{array}\right)_{n \times n} \left( \begin{array}{c} a_1 \\ a_2 \\ \vdots \\ a_n \end{array} \right)_{n \times 1} \]

The following code simulates our coin flips' running proportions. As said above, when \( n \) is large, these proportions will come ever nearer to \( p \). Notice that 'rbinom' simulates random binomial random variables. In the case that the max value for the binomial is 1, it may only take on 1 or 0. In this case, the probability of 1 is \( p \). I've asked for \( n \) of them. Notice that in the return line, I take the part sums and then divide by a vector, \( \left(1,2, \cdots , n\right)^T \).

flips = function(n, p) {
    x = rbinom(n, 1, p)
    return(partSums(x)/1:n)
}

In R, dividing a vector by a vector results in element-wise division, as follows.

\( (b_1,b_2,\cdots,b_n)^T/(1,2,\cdots,n)^T = \left( \frac{b_1}{1} , \frac{b_2}{2} , \cdots , \frac{b_n}{n} \right)^T \)

Now we can plot our coin flips. We'll do so in a way that demonstrates the LLN. The following code will draw \( p \) as a thick grey line, and plot 6 alternative paths.

plotFlips = function(n, p) {
    plot(1:n, rep(p, n), ylim = c(0, 1), type = "l", xlab = "flips", ylab = "proportion Heads", 
        main = paste("Evan's coin flip simulator, (n,p) = (", n, ",", p, ")", 
            sep = ""), col = "grey", lwd = 10)
    lines(flips(n, p), col = "black", lwd = 2)
    lines(flips(n, p), col = "magenta", lwd = 2)
    lines(flips(n, p), col = "cyan", lwd = 2)
    lines(flips(n, p), col = "red", lwd = 2)
    lines(flips(n, p), col = "blue", lwd = 2)
    lines(flips(n, p), col = "green", lwd = 2)
}

Now coin flip simulation is very easy!

plotFlips(100, 0.71)