# Calculate a binomial in Python to determine the probability

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.

• Because the question asks you to calculate, running a short simulation is probably not what it intended you to do. You can compute the answer directly (and easily) using a formula for the Binomial distribution.
– whuber
Jan 31, 2020 at 16:54
• Thank you for your response whuber. I thought they were interchangeable. Jan 31, 2020 at 18:07
• Whuber, I am trying to plot only the birthweights I want but I managed to plot all the binomial distribution which is incorrect. ax = sns.distplot(binomial, kde=True, norm_hist=True, color='black', hist_kws={"linewidth": 15,'alpha':1, 'color': 'g'}) ax.set(xlabel='Number of births', ylabel='Frequency') Jan 31, 2020 at 18:48

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.

• Thank you JohanC I upvoted your answer Feb 1, 2020 at 17:43
• JohanC what is the difference between repeats and runs? Feb 1, 2020 at 17:49
• The runs means: the simple experiment does 100 random coin flips with probability p, and runs this 50 times. It then counts how many of those flips are between 7 and 13. The large experiment repeats the small experiment 2000 times and looks at those results (each individual result is a number close to the desired result of 0.7589...). Feb 1, 2020 at 18:04
• So the shape of the distribution is symmetrical I think the mean=median=mode Feb 1, 2020 at 19:21
• No, it can't be symmetrical. The highest possible value is 1. The lowest 0. 1 means all runs give a result between 7 and 13. 0 means no run at all ended between those limits. So the distribution is more stretched out at the left. The mode is more likely to be at the right. Feb 1, 2020 at 19:35

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.

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.

• thank you for your answer! Could please tell me more why when you increase your runs it also decreases variability? Jan 31, 2020 at 18:06
• I understand this interpretation, but a key word in the question suggests it's not the intended one: "Plot these probabilities..." To me, that says the question seeks the individual probabilities for $7, 8, \ldots, 13$ low-birthweight babies.
– whuber
Jan 31, 2020 at 18:06
• What about in the simulation having 0 low-birthweight babies? I get 0 always as an output. Jan 31, 2020 at 20:32
• If you want to calculate $P(X=0)$, @JoãoVitorGomes it's just sum(binomial==0)/nruns Feb 1, 2020 at 12:00
• @gunes, thank you for your response. I had done this, but output was always 0. I then calculated with binomial formula and is such a small probability that was not showing enough decimal places. Feb 1, 2020 at 13:36