Calculate a binomial in Python to determine the probability

Question

Calculate a binomial in Python to determine the probability of getting: 7, 8, 9, 10, 11, 12, or 13 low‐birthweight babies in 100 deliveries, if the probability of this outcome is 0.1. Arrange the values in a table. Plot these probabilities (vertical axis) against number of low birthweight babies. Comment on the shape of this graph.

What I did was:

runs = 50
n= 100
p = .1
binomial = np.random.binomial(n, p, runs)
prob_7 = sum(binomial==7)/runs
print('The probability of 7 premature babies is: ' + str(prob_7))

Output: The probability of 7 premature babies is: 0.06

But I don't know if this is correct, because I am not sure about the quantity of runs (experiments) I set of 50.

Answer 1

Plotting a seaborn distplot needs an adjustment, as it is primarily meant for continuous distributions. The distplot will put the data in 16 equally size bins, that don't align with the integer numbers. For discrete distributions, distplot would need explicit bins, e.g. range(30). However, with that many bins, the default calculated kde will not be as desired.

You can also use numpy to count how many values there are of each frequency and then plot a bar chart. As the theoretical values can be easily calculated (comb(n, k)*(p**k)*((1-p)**(n-k))), we can plot these on top as a comparison. Notice that even for 10,000 runs there still is a visible difference between the theoretical and the experimental value.

Here is some sample code that can be used as a base for experimentation. It is using 10000 runs.

from matplotlib import pyplot as plt
import numpy as np
import seaborn as sns
from scipy.special import comb

runs = 10000
n = 100
p = 0.1
binomial = np.random.binomial(n, p, runs)

fig, (ax1, ax2) = plt.subplots(ncols=2, figsize=(14, 6))
sns.distplot(binomial, kde=True, norm_hist=True, color='black', bins=range(26),
             hist_kws={"linewidth": 15, 'alpha': 1, 'color': 'g'}, ax=ax1)

values, value_count = np.unique(binomial, return_counts=True)
frequencies = value_count / runs
filter = (values >= 7) & (values <= 13)
ax2.bar(values[filter], frequencies[filter], color='crimson')
ax2.bar(values[~filter], frequencies[~filter], color='dodgerblue')

k = np.arange(101)
theoretical = comb(n, k) * (p ** k) * ((1-p) ** (n-k))
ax2.plot(k[:26], theoretical[:26], color='limegreen', markersize=4, marker='o', alpha=0.6, zorder=3)

for ax in (ax1, ax2):
    ax.set_ylabel('Frequency')
    ax.set_xlabel('Number of births')
    ax.set_xticks(range(26))
    ax.set_xlim(0, 26)
    ax.set_ylim(0, frequencies.max() * 1.1)

print('theoretical sum for 7 ≤ value ≤ 13 :', sum(theoretical[7:14]))
print('experimental sum for 7 ≤ value ≤ 13 :', sum([freq for freq, i in zip(frequencies, values) if 7 <= i <= 13]))

plt.show()

Output:

theoretical sum for 7 ≤ value ≤ 13 :  0.7589675919647877
experimental sum for 7 ≤ value ≤ 13 : 0.7580999999999999

The experiment can be taken one step further. Repeat doing the experiment over 50 runs. And study the distribution of the outcome. Note that this again only has a discrete result: averaging 50 integers. So, again we need to provide some bins adjusted to distribution: numbers with exactly 2 decimals, the last decimal being even. The code adds a green line at the mean of the distribution, and a yellow line at the theoretical solution.

from matplotlib import pyplot as plt
import seaborn as sns
import numpy as np
from scipy.special import comb

runs = 50
repeats = 2000
n = 100
p = 0.1
binomial = np.random.binomial(n, p, (repeats, runs))
between7and13 = np.sum((binomial >= 7) & (binomial <= 13), axis=1) / runs
print("Mean :", between7and13.mean())
print("Median :", np.median(between7and13))
k = np.arange(101)
theoretical = comb(n, k) * (p ** k) * ((1 - p) ** (n - k))
print('theoretical sum for 7 ≤ value ≤ 13 :', sum(theoretical[7:14]))

sns.distplot(between7and13, kde=True, norm_hist=True, color='black', bins=np.arange(0.55, 0.95, 0.02),
             hist_kws={'alpha': 1, 'color': 'tomato'})
plt.axvline(np.median(between7and13), color='gold', lw=1)
plt.axvline(sum(theoretical[7:14]), color='lime', lw=1)
plt.title(f'Outcome of {repeats} experiments')
plt.xlabel(f'Mean of "7 ≤ Number of births ≤ 13" over {runs} runs')
plt.show()

Note that the mean (and the median) is very close to the desired theoretical result. The distribution itself is rather wide though, and the mode can be higher than desired.

Answer 2

You're asked to calculate $P(7\leq X \leq 13)$, but you only find $P(X=7)$, where $X\sim \text{Bin(n,p)}$ denotes the number of low-birth weight babies.

So, your script should be:

sum((binomial>=7) & (binomial <= 13))/runs

And, it's better to increase your number of runs to decrease your estimate's variance. You can also calculate it analytically as @whuber commented.

Edit: If you want to calculate $P(X=x)$, then your code is correct other than the plotting. Just increase your number of runs, and add some plotting code.