import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import scipy.stats as spicy
df = pd.read_csv("https://raw.githubusercontent.com/roualdes/data/refs/heads/master/penguins.csv")
rng = np.random.default_rng()
ndx = ~df["bill_length_mm"].isna()
N = np.sum(ndx)
x = df.loc[ndx, "bill_length_mm"]
R = 1_000
means = np.zeros(R)
for r in range(R):
    rx = rng.choice(x, size = N) # uniformly re-sampling our data N times
    means[r] = np.mean(rx)
np.quantile(means, [0.025, 0.975])
array([43.37855994, 44.51141813])
x
0      39.1
1      39.5
2      40.3
3      36.7
4      39.3
       ... 
337    55.8
338    43.5
339    49.6
340    50.8
341    50.2
Name: bill_length_mm, Length: 342, dtype: float64
def bootstrap(arr, T, R = 1_000):
    N = np.size(arr)
    Ts = np.zeros(R)
    for r in range(R):
        # rx = rng.choice(arr, size = N)
        idx = rng.integers(N, size = N)
        if type(arr) is np.ndarray:
            Ts[r] = T(arr[idx])
        else:
            Ts[r] = T(arr.iloc[idx])
    return Ts
mns = bootstrap(x, np.mean)
np.quantile(mns, [0.025, 0.975])
array([43.35813596, 44.48751462])
sds = bootstrap(x, np.std)
np.quantile(sds, [0.025, 0.975])
array([5.13680784, 5.73536427])
sds = bootstrap(x, np.median)
np.quantile(sds, [0.025, 0.975])
array([42.9, 45.4])
plt.hist(means, density = True, histtype = "step");

# bootstrap confidence interval
np.quantile(means, [0.025, 0.975])
array([43.36223684, 44.48954678])
# compared to the standard method
xbar = np.mean(df["bill_length_mm"])
s = np.std(df["bill_length_mm"])
N = np.sum(~df["bill_length_mm"].isna())
t = spicy.t(df = N - 1).ppf([0.025, 0.975]) # ppf(0.975)
xbar + t * s / np.sqrt(N)
array([43.34209692, 44.50176273])