SciPy Bernoulli random number generation, different behavior within a for loop for different sample sizes

Question

The purpose of this code is to demonstrate CLT.

If I do the following:

num_samples = 10000

sample_means = np.empty(num_samples)
for i in range(num_samples):
  mean = np.mean(st.bernoulli.rvs(p=0.5, size=100))
  sample_means[i] = mean

sample_demeaned = np.subtract(sample_means, 0.5)
denominator = np.divide(0.5, np.sqrt(100))
z_ed = np.divide(sample_demeaned, denominator)

plt.hist(z_ed, bins=40, edgecolor='k', density=True)

x = np.linspace(st.norm.ppf(0.001), st.norm.ppf(0.999), 10000)
y = st.norm.pdf(x)
plt.plot(x, y, color='red')

I get:

However, if I try to do it with a for loop for different sample sizes:

num_samples = 10000
sample_sizes = np.array([5, 20, 75, 100])

sample_std_means = np.empty(shape=(num_samples, len(sample_sizes)))


for col, size in enumerate(sample_sizes):
  sample_means = np.empty(num_samples)
  for i in range(num_samples):
    mean = np.mean(st.bernoulli.rvs(p=0.5, size=size))
    sample_means[i] = mean
  sample_demeaned = np.subtract(sample_means, 0.5)
  denominator = np.divide(0.5, np.sqrt(size))
  z_ed = np.divide(sample_demeaned, denominator)
  sample_std_means[:, col] = sample_means

And then plot each of them in a 2x2 grid:

fig, axes = plt.subplots(nrows=2, ncols=2, figsize=(10, 7))

x = np.linspace(st.norm.ppf(0.001), st.norm.ppf(0.999), 10000)
y = st.norm.pdf(x)

for i, ax in enumerate(axes.flatten()):
  ax.hist(sample_std_means[i], bins=40, edgecolor='k', color='midnightblue')
  ax.set_ylabel('Density')
  ax.set_xlabel(f'n = {sample_sizes[i]}')
  ax.plot(x, y, color='red')
  ax.set_xlim((-3, 3))

plt.show()

I get the following image:

I cannot debug the discrepancy here. Any help is highly appreciated.

Please note that scipy.stats and numpy have been imported as st and np respectively in both code blocks.

user_12046409 · Answer

First, note that one numpy's strong points is that it allows operations which mix arrays and single numbers. This is called broadcasting. So, for example sample_demeaned = np.subtract(sample_means, 0.5) can be written more concise as sample_demeaned = sample_means - 0.5.

Several issues are going wrong:

sample_std_means[:, col] = sample_means should use the just calculated z_ed instead of sample_means.
ax.hist(sample_std_means[i], ...) uses the i'th row of the array. That row only contains 4 elements. You'd want sample_std_means[;,i] to take the i'th column.
The pdf is drawn in its normalized form (with an area below the curve equal to one). However, the histogram's height is proportional to the number of samples. Its total area is num_samples * bin_width, where the histogram's default bin width is the length from the first to the last element divided by the number of bins. To get both the pdf and histogram with similar sizes, either the histogram should be normalized (using density=True) or the pdf should be multiplied by the expected area of the histogram.

import numpy as np
import scipy.stats as st
import matplotlib.pyplot as plt

num_samples = 10000
sample_sizes = np.array([5, 20, 75, 100])

sample_std_means = np.empty(shape=(num_samples, len(sample_sizes)))

for col, size in enumerate(sample_sizes):
    sample_means = np.empty(num_samples)
    for i in range(num_samples):
        sample_means[i] = np.mean(st.bernoulli.rvs(p=0.5, size=size))
    sample_demeaned = sample_means - 0.5
    z_ed = sample_demeaned / (0.5 / np.sqrt(size))
    sample_std_means[:, col] = z_ed

fig, axes = plt.subplots(nrows=2, ncols=2, figsize=(10, 7))

x = np.linspace(st.norm.ppf(0.001), st.norm.ppf(0.999), 1000)
y = st.norm.pdf(x)

for i, ax in enumerate(axes.flatten()):
    ax.hist(sample_std_means[:, i], bins=40, edgecolor='k', color='midnightblue', density=True)
    ax.set_ylabel('Density')
    ax.set_xlabel(f'n = {sample_sizes[i]}')
    # bin_width = (sample_std_means[:, i].max() - sample_std_means[:, i].min()) / 40
    # ax.plot(x, y * num_samples * bin_width, color='red')
    ax.plot(x, y, color='red')
    ax.set_xlim((-3, 3))
plt.show()

Now note the weird empty bars in the histograms. A histogram works best for continuous distributions. But the mean of n Bernoulli trials can have at most n+1 different outcomes. When all trials would be True, the mean would be n/n = 1. When all would be False, the mean would be 0. Combined, the possible means are 0, 1/n, 2/n, ..., 1. The histogram of such a discrete distribution should take these values into account for the boundaries between the bins.

The following code creates a scatter plot, using the position of the means and a random y-value to visualize how many there are per x. Also, the position of the bin boundaries is calculated and visualized by dotted vertical lines.

fig, axes = plt.subplots(nrows=2, ncols=2, figsize=(10, 7))

for i, ax in enumerate(axes.flatten()):
    ax.scatter(sample_std_means[:, i], np.random.uniform(0, 1, num_samples), color='r', alpha=0.5, lw=0, s=1)
    # there are n+1 possible mean values for n bernoulli trials
    # n+2 boundaries will be needed to separate the bins
    bins = np.arange(-1, sample_sizes[i]+1) / sample_sizes[i]
    bins += (bins[1] - bins[0]) / 2 # shift half a bin
    bins -= 0.5 # subtract the mean
    bins /= (0.5 / np.sqrt(sample_sizes[i])) # correction factor
    for b in bins:
        ax.axvline(b, color='g', ls=':')
    ax.set_xlabel(f'n = {sample_sizes[i]}')
    ax.set_xlim((-3, 3))

And here are the histograms using these bins:

    ax.hist(sample_std_means[:, i], bins=bins, edgecolor='k', color='midnightblue', density=True)

SciPy Bernoulli random number generation, different behavior within a for loop for different sample sizes

1 Answers