3. Resampling methods
3.1 Introduction to resampling methods
Probability example
In this exercise, we will review the difference between sampling with and without replacement. We will calculate the probability of an event using simulation, but vary our sampling method to see how it impacts probability.
Consider a bowl filled with colored candies – three blue, two green, and five yellow. Draw three candies at random, with replacement and without replacement. You want to know the probability of drawing a yellow candy on the third draw given that the first candy was blue and the second candy was green.
# Set up the bowl
success_rep, success_no_rep, sims = 0, 0, 10000
bowl = ['b'] * 3 + ['g'] * 2 + ['y'] * 5
for i in range(sims):
# Sample with and without replacement & increment success counters
sample_rep = np.random.choice(bowl, replace=True, size=3)
sample_no_rep = np.random.choice(bowl, replace=False, size=3)
if (sample_rep[0] == 'b') & (sample_rep[1] == 'g') & (sample_rep[2] == 'y'):
success_rep += 1
if (sample_no_rep[0] == 'b') & (sample_no_rep[1] == 'g') & (sample_no_rep[2] == 'y'):
success_no_rep += 1
# Calculate probabilities
prob_with_replacement = success_rep / sims
prob_without_replacement = success_no_rep / sims
print("Probability with replacement = {}, without replacement = {}".format(prob_with_replacement, prob_without_replacement))
# Probability with replacement = 0.0266, without replacement = 0.0415
Does the difference between sampling with and without replacement make sense now?
3.2 Bootstrapping
Running a simple bootstrap
Welcome to the first exercise in the bootstrapping section. We will work through an example where we learn to run a simple bootstrap. As we saw in the video, the main idea behind bootstrapping is sampling with replacement.
Suppose you own a factory that produces wrenches. You want to be able to characterize the average length of the wrenches and ensure that they meet some specifications. Your factory produces thousands of wrenches every day, but it’s infeasible to measure the length of each wrench. However, you have access to a representative sample of 100 wrenches. Let’s use bootstrapping to get the 95% confidence interval (CI) for the average lengths.
Examine the list wrench_lengths, which has 100 observed lengths of wrenches, in the shell.
# Draw some random sample with replacement and append mean to mean_lengths.
mean_lengths, sims = [], 1000
for i in range(sims):
temp_sample = np.random.choice(wrench_lengths, replace=True, size=100)
sample_mean = np.mean(temp_sample)
mean_lengths.append(sample_mean)
# Calculate bootstrapped mean and 95% confidence interval.
boot_mean = np.mean(mean_lengths)
boot_95_ci = np.percentile(mean_lengths, [2.5, 97.5])
print("Bootstrapped Mean Length = {}, 95% CI = {}".format(boot_mean, boot_95_ci))
# Bootstrapped Mean Length = 10.027059690070363, 95% CI = [ 9.78662216 10.24854356]
Non-standard estimators
In the last exercise, you ran a simple bootstrap that we will now modify for more complicated estimators.
Suppose you are studying the health of students. You are given the height and weight of 1000 students and are interested in the median height as well as the correlation between height and weight and the associated 95% CI for these quantities. Let’s use bootstrapping.
Examine the pandas DataFrame df with the heights and weights of 1000 students. Using this, calculate the 95% CI for both the median height as well as the correlation between height and weight.
# Sample with replacement and calculate quantities of interest
sims, data_size, height_medians, hw_corr = 1000, df.shape[0], [], []
for i in range(sims):
tmp_df = df.sample(n=data_size, replace=True)
height_medians.append(np.median(tmp_df.heights))
hw_corr.append(tmp_df.weights.corr(tmp_df.heights))
# Calculate confidence intervals
height_median_ci = np.percentile(height_medians, [2.5, 97.5])
height_weight_corr_ci = np.percentile(hw_corr, [2.5, 97.5])
print("Height Median CI = {} \nHeight Weight Correlation CI = {}".format( height_median_ci, height_weight_corr_ci))
# Height Median CI = [5.25262253 5.55928686]
# Height Weight Correlation CI = [0.93892136 0.95103152]
Bootstrapping regression
Now let’s see how bootstrapping works with regression. Bootstrapping helps estimate the uncertainty of non-standard estimators. Consider the R2 statistic associated with a regression. When you run a simple least squares regression, you get a value for R2. But let’s see how can we get a 95% CI for R2.
Examine the DataFrame df with a dependent variable yy and two independent variables X1 and X2 using df.head(). We’ve already fit this regression with statsmodels (sm) using:
reg_fit = sm.OLS(df['y'], df.iloc[:,1:]).fit()
Examine the result using reg_fit.summary() to find that R2=0.3504. Use bootstrapping to calculate the 95% CI.
df.head()
y Intercept X1 X2
0 1.217851 1.0 0.696469 0.286139
1 1.555250 1.0 0.226851 0.551315
2 0.888520 1.0 0.719469 0.423106
3 1.736052 1.0 0.980764 0.684830
4 1.632073 1.0 0.480932 0.392118
rsquared_boot, coefs_boot, sims = [], [], 1000
reg_fit = sm.OLS(df['y'], df.iloc[:,1:]).fit()
# Run 1K iterations
for i in range(sims):
# First create a bootstrap sample with replacement with n=df.shape[0]
bootstrap = df.sample(n=df.shape[0], replace=True)
# Fit the regression and append the r square to rsquared_boot
rsquared_boot.append(sm.OLS(bootstrap['y'],bootstrap.iloc[:,1:]).fit().rsquared)
# Calculate 95% CI on rsquared_boot
r_sq_95_ci = np.percentile(rsquared_boot, [2.5, 97.5])
print("R Squared 95% CI = {}".format(r_sq_95_ci))
# R Squared 95% CI = [0.31089312 0.40543591]
3.3 Jackknife resampling
Basic jackknife estimation – mean
Jackknife resampling is an older procedure, which isn’t used as often compared as bootstrapping. However, it’s still useful to know how to run a basic jackknife estimation procedure. In this first exercise, we will calculate the jackknife estimate for the mean. Let’s return to the wrench factory.
You own a wrench factory and want to measure the average length of the wrenches to ensure that they meet some specifications. Your factory produces thousands of wrenches every day, but it’s infeasible to measure the length of each wrench. However, you have access to a representative sample of 100 wrenches. Let’s use jackknife estimation to get the average lengths.
Examine the variable wrench_lengths in the shell.
# Leave one observation out from wrench_lengths to get the jackknife sample and store the mean length
mean_lengths, n = [], len(wrench_lengths)
index = np.arange(n)
for i in range(n):
jk_sample = wrench_lengths[index != i]
mean_lengths.append(np.mean(jk_sample))
# The jackknife estimate is the mean of the mean lengths from each sample
mean_lengths_jk = np.mean(np.array(mean_lengths))
print("Jackknife estimate of the mean = {}".format(mean_lengths_jk))
Jackknife confidence interval for the median
In this exercise, we will calculate the jackknife 95% CI for a non-standard estimator. Here, we will look at the median. Keep in mind that the variance of a jackknife estimator is n-1 times the variance of the individual jackknife sample estimates where n is the number of observations in the original sample.
Returning to the wrench factory, you are now interested in estimating the median length of the wrenches along with a 95% CI to ensure that the wrenches are within tolerance.
Let’s revisit the code from the previous exercise, but this time in the context of median lengths. By the end of this exercise, you will have a much better idea of how to use jackknife resampling to calculate confidence intervals for non-standard estimators.
n = 100
# Leave one observation out to get the jackknife sample and store the median length
median_lengths = []
for i in range(n):
jk_sample = wrench_lengths[index != i]
median_lengths.append(np.median(jk_sample))
median_lengths = np.array(median_lengths)
# Calculate jackknife estimate and it's variance
jk_median_length = np.mean(median_lengths)
jk_var = (n-1)*np.var(median_lengths)
# Assuming normality, calculate lower and upper 95% confidence intervals
jk_lower_ci = jk_median_length - 1.96*np.sqrt(jk_var)
jk_upper_ci = jk_median_length + 1.96*np.sqrt(jk_var)
print("Jackknife 95% CI lower = {}, upper = {}".format(jk_lower_ci, jk_upper_ci))
# Jackknife 95% CI lower = 9.138592467547202, upper = 10.754868069037098
3.4 Permutation testing
Generating a single permutation
In the next few exercises, we will run a significance test using permutation testing. As discussed in the video, we want to see if there’s any difference in the donations generated by the two designs – A and B. Suppose that you have been running both the versions for a few days and have generated 500 donations on A and 700 donations on B, stored in the variables donations_A and donations_B.
We first need to generate a null distribution for the difference in means. We will achieve this by generating multiple permutations of the dataset and calculating the difference in means for each case.
First, let’s generate one permutation and calculate the difference in means for the permuted dataset.
# Concatenate the two arrays donations_A and donations_B into data
len_A, len_B = len(donations_A), len(donations_B)
data = np.concatenate([donations_A, donations_B])
# Get a single permutation of the concatenated length
perm = np.random.permutation(len(donations_A) + len(donations_B))
# Calculate the permutated datasets and difference in means
permuted_A = data[perm[:len(donations_A)]]
permuted_B = data[perm[len(donations_A):]]
diff_in_means = np.mean(permuted_A) - np.mean(permuted_B)
print("Difference in the permuted mean values = {}.".format(diff_in_means))
# Difference in the permuted mean values = -0.13886241452516757.
Hypothesis testing – Difference of means
We want to test the hypothesis that there is a difference in the average donations received from A and B. Previously, you learned how to generate one permutation of the data. Now, we will generate a null distribution of the difference in means and then calculate the p-value.
For the null distribution, we first generate multiple permuted datasets and store the difference in means for each case. We then calculate the test statistic as the difference in means with the original dataset. Finally, we calculate the p-value as twice the fraction of cases where the difference is greater than or equal to the absolute value of the test statistic (2-sided hypothesis). A p-value of less than say 0.05 could then determine statistical significance.
perm
array([[ 850, 473, 1067, ..., 592, 644, 944],
[ 594, 535, 125, ..., 1198, 1052, 233],
[ 494, 246, 809, ..., 179, 448, 953],
...,
[ 923, 888, 393, ..., 314, 40, 258],
[ 915, 1140, 953, ..., 20, 526, 272],
[ 925, 480, 1102, ..., 371, 85, 379]])
permuted_A_datasets
array([[2.49959605e+00, 1.94959571e+00, 4.10970395e+00, ...,
1.28042313e+00, 3.09147765e+00, 3.00804073e-01],
...,
[6.33883845e-01, 5.64755610e-01, 1.75817616e+00, ...,
5.46401636e-01, 1.97500056e+00, 5.43297027e+00]])
reps
1000
# Generate permutations equal to the number of repetitions
perm = np.array([np.random.permutation(len(donations_A) + len(donations_B)) for i in range(reps)])
permuted_A_datasets = data[perm[:, :len(donations_A)]]
permuted_B_datasets = data[perm[:, len(donations_A):]]
# Calculate the difference in means for each of the datasets
samples = np.mean(permuted_A_datasets, axis=1) - np.mean(permuted_B_datasets, axis=1)
# Calculate the test statistic and p-value
test_stat = np.mean(donations_A) - np.mean(donations_B)
p_val = 2*np.sum(samples >= np.abs(test_stat))/reps
print("p-value = {}".format(p_val))
# p-value = 0.002
Hypothesis testing – Non-standard statistics
In the previous two exercises, we ran a permutation test for the difference in mean values. Now let’s look at non-standard statistics.
Suppose that you’re interested in understanding the distribution of the donations received from websites A and B. For this, you want to see if there’s a statistically significant difference in the median and the 80th percentile of the donations. Permutation testing gives you a wonderfully flexible framework for attacking such problems.
Let’s go through running a test to see if there’s a difference in the median and the 80th percentile of the distribution of donations. As before, you’re given the donations from the websites A and B in the variables donations_A and donations_B respectively.
# Calculate the difference in 80th percentile and median for each of the permuted datasets (A and B)
samples_percentile = np.percentile(permuted_A_datasets, 80, axis=1) - np.percentile(permuted_B_datasets, 80, axis=1)
samples_median = np.median(permuted_A_datasets, axis=1) - np.median(permuted_B_datasets, axis=1)
# Calculate the test statistic from the original dataset and corresponding p-values
test_stat_percentile = np.percentile(donations_A, 80) - np.percentile(donations_B, 80)
test_stat_median = np.median(donations_A) - np.median(donations_B)
p_val_percentile = 2*np.sum(samples_percentile >= np.abs(test_stat_percentile))/reps
p_val_median = 2*np.sum(samples_median >= np.abs(test_stat_median))/reps
print("80th Percentile: test statistic = {}, p-value = {}".format(test_stat_percentile, p_val_percentile))
print("Median: test statistic = {}, p-value = {}".format(test_stat_median, p_val_median))
# 80th Percentile: test statistic = 1.6951624543447839, p-value = 0.026
# Median: test statistic = 0.6434965714975927, p-value = 0.014
way to be a data scientist 