Trying to understand Bootstrapping w/ Python

bootstrappython

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!

Best Answer

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.

Related Question