Simulate an experiment carried out ntrials times in which, for each experiment, n coins are tossed and the total number of heads each time is recorded.
Plot the results of the simulation on a suitable histogram, and compare with the expected binomial distribution of heads.
Solution P6.6.1
The simulation can be carried out in two ways. Without using np.random.binomial() we can simply sample an integer randomly from $[0,1]$ 10 times and take the sum: if we associate 1 with the outcome heads and 0 with the outcome tails we then have the number of heads out of 10 coin-flips.
Alternatively, since the number of heads will be binomially-distributed, we can sample ntrials times from the binomial distribution described by $n=10$, $p=\frac{1}{2}$.
from math import factorial
import numpy as np
import matplotlib.pyplot as plt
p, n = 0.5, 10
ntrials = 1000
# Simulate the experiment by randomly selecting an integer out of [0, 1]
# n times for each of the ntrials
nheads = np.zeros(ntrials)
for i in range(ntrials):
nheads[i] = np.sum(np.random.randint(0, 2, n))
# Using NumPy to sample from the binomial distribution
binomial = np.random.binomial(n, p, ntrials)
bin_edges = np.arange(n + 2)
hist, bins, patches = plt.hist(
nheads, bin_edges, density=True, align="left", label="Simulated"
)
plt.hist(
binomial,
bins,
density=True,
align="left",
alpha=0.8,
label="Binomial draws",
)
k = np.arange(n + 1)
bdist = np.zeros(n + 1)
for r in k:
nCr = factorial(n) / factorial(r) / factorial(n - r)
bdist[r] = nCr * p**r * (1.0 - p) ** (n - r)
plt.plot(k, bdist, "o", color="r", label="Binomial theorem")
plt.xlim(0, n)
plt.xlabel("Number of heads")
plt.ylabel("Probability")
plt.legend()
plt.show()
In the above code, the histogram figure is be plotted with pyplot, taking care to set up the bin edges properly, normalize the data and to align the histogram bars correctly. The "exact" probabilities, bdist, are calculated from the binomial distribution formula given above.
Coin flipping simulation and comparison with the Binomial theorem.
