Trying to understand Bootstrapping w/ Python I am trying to understand when (and how) to use Bootstrapping. I read on some other questions that you shouldn't use Bootstrapping for small confidence intervals, and I wanted to try it by myself.
My idea was to:

*

*take multiple samples from a normal population (with mean 100 and std 5)

*bootstrap each of those samples and the mean of those samples in an array

*Based on those means, build a confidence interval for the mean of my population

But when I run the following code:
import numpy as np
import pandas as pd
import scipy.stats as stats

d = {'Sample #': [], 'Mean': [], 'Std':[], 'Lower Bound':[], 'Higher Bound': [],'Test Result': []}
df = pd.DataFrame(data=d)

sample_size = 100
population_mean = 100
population_std = 5

for a in range(100):
  new_sample = np.random.normal(population_mean, population_std, sample_size)
  bootstrapped_sample_means = []
  for z in range(1000):
      x = np.random.choice(new_sample, size=sample_size, replace=True)
      bootstrapped_sample_means .append(x.mean())
  mean, std = np.mean(bootstrapped_sample_means), np.std(bootstrapped_sample_means)
  lower_bound = mean - 1.96*std/np.sqrt(1000)
  higher_bound = mean + 1.96*std/np.sqrt(1000)
  row = {'Sample #': a, 'Mean': mean, 'Std':std, 'Lower Bound':lower_bound, 'Higher Bound': higher_bound,'Test Result': (population_mean > higher_bound) | (population_mean < lower_bound)}
  df = df.append(row, ignore_index = True)

I end up finding that even for larger sample size - most of the times the population mean is not inside the confidence intervals.
I guess something is wrong in my code (or in my understanding of how things work) so any help is greatly appreciated!
 A: This is a comment to complement @Lynn's answer (+1) which explains the big mistake in your implementation of the bootstrap.
There is another error. It's more subtle as it doesn't have an obvious effect on the coverage of the confidence interval: You should use the sample statistic, not the average of the bootstrapped statistics as the center of the bootstrap confidence interval.
Here is the updated python code:
import numpy as np

np.random.seed(1234)

sample_size = 100
bootstrap_reps = 1000
population_mean = 100
population_std = 5

sample = np.random.normal(population_mean, population_std, sample_size)

# Function to generate bootstrap samples
bootstrap_sample = lambda: np.random.choice(sample, size=sample_size, replace=True)

# The sample mean is an estimator of the population mean
estimator = np.mean

sample_statistic = estimator(sample)
bootstrapped_statistics = [
    estimator(bootstrap_sample())
    for _ in range(bootstrap_reps)
]

# Use the sample statistic [here the sample mean],
# not the average of the bootstrapped statistics,
# to construct the confidence interval
sample_statistic, np.mean(bootstrapped_statistics)
# (100.17556141562721, 100.18181492004352)

# Lower and upper limit of the 95% confidence interval for the population parameter
z_alpha = 1.96
ci_lower = sample_statistic - z_alpha * np.std(bootstrapped_statistics)
ci_upper = sample_statistic + z_alpha * np.std(bootstrapped_statistics)

A: When you calculate your lower and upper bounds, you are dividing your standard deviations of your means by $\sqrt{1000}$, which I assume you are doing because you want to convert your standard deviation to a standard error. But the standard deviation you are using is the standard deviation (or standard error) of the mean, so you don't need to do this. If you replace those lines with:
lower_bound = mean - 1.96*std
higher_bound = mean + 1.96*std

Your results should look better.
