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}$.
import numpy as np
import pylab
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.random_integers(0, 1, 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 = pylab.hist(nheads, bin_edges, normed=True, align='left')
pylab.hist(binomial, bins, normed=True, color='y', align='left', alpha=0.5)
k = np.arange(n+1)
bdist = np.zeros(n+1)
for r in k:
nCr = np.math.factorial(n) / np.math.factorial(r) / np.math.factorial(n-r)
bdist[r] = nCr * p**r * (1. - p)**(n-r)
pylab.plot(k, bdist, 'o', color='r')
pylab.xlim(0,n)
pylab.show()
In the above code, the histogram figure is be plotted with pylab
, 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.