Resampling Techniques

Resampling involves techniques that repeatedly draw new samples from an original dataset. These are useful techniques to be aware of as they can be useful when we have small sample sizes but, the underlying distribution is not normal. These types of techniques are best understood by looking at examples.

Bootstrapping

Let’s start with an example of a resampling technique called “bootstrap” resampling. Bootstrap resampling or bootstrapping is one of the most common resampling techniques and involves random resampling with replacement (this means that you can resample the same data more than once). There are many different ways to use bootstrapping, here are a couple:

• Bootstrapping for hypothesis testing involves constructing a large number of resamples of the original dataset. These resamples should be of equal size to a specific sample of interest (which is itself drawn from the original dataset) and be drawn by random sampling with replacement from the original dataset. In this way, you can construct a sampling distribution and you never need to assume anything about the underlying distribution of the data as it is already built-in to the original dataset. This approach can be used for small sample sizes and is often done when you have a long observational climatology or climate model control integration to draw samples from.

• Bootstrapping for confidence intervals involves resampling your specific sample to create a sampling distribution based on your specific sample data. From this sampling distribution you can compute confidence intervals.

Let’s look at an example. For this example, we are going to examine a large-scale mode of atmospheric variability in the Southern Hemisphere, the Southern Annular Mode (SAM). You can download your own copy of the SAM data file here.

# Load packages
import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl
mpl.rc('font',size=16) #set default font size and weight for plots
import scipy.stats as st
import scipy.io as sio
import time
from IPython import display

# Load data (we are loading a summer season time series of the SAM for 1957-2025)
with open('SAM_summer.csv', 'r', encoding='utf-8-sig') as f:
    SAM = np.genfromtxt(f, delimiter=',')

X = SAM[:,1]
TIME_SAM = SAM[:,0]

A good habit to get into is the quickly check the size of your dataset before you begin.

print(X.shape)
print(TIME_SAM)

(69,)
[1957. 1958. 1959. 1960. 1961. 1962. 1963. 1964. 1965. 1966. 1967. 1968.
 1969. 1970. 1971. 1972. 1973. 1974. 1975. 1976. 1977. 1978. 1979. 1980.
 1981. 1982. 1983. 1984. 1985. 1986. 1987. 1988. 1989. 1990. 1991. 1992.
 1993. 1994. 1995. 1996. 1997. 1998. 1999. 2000. 2001. 2002. 2003. 2004.
 2005. 2006. 2007. 2008. 2009. 2010. 2011. 2012. 2013. 2014. 2015. 2016.
 2017. 2018. 2019. 2020. 2021. 2022. 2023. 2024. 2025.]

So, we have 69 summer’s from 1957-2025. Let’s plot the data to see what it looks like.

# Plot time series
plt.figure(figsize=(10,6))
plt.plot(TIME_SAM,X,color = 'black', linewidth = 1.5)
plt.xlabel('Year')
plt.ylabel('SAM Index')
plt.title('Time Series of SAM Index')
plt.ylim(-5,5)
plt.xlim(min(TIME_SAM), max(TIME_SAM))
plt.axhline(0,color='gray')

<matplotlib.lines.Line2D at 0x177c0ff10>

You should be able to clearly see a upward trend in the data over this time period, which is related to the depletion of stratospheric ozone. Let’s compute the sample means for the two halves of the time series, 1957-1990 and 1992-2025.

# mean 1954-1971
SAM1 = np.mean(X[0:35])

# mean 1989-2006
SAM2 = np.mean(X[36:])

print(np.round(SAM1,3),np.round(SAM2,3))

SAM_diff = SAM2 - SAM1

print("The difference is", np.round(SAM_diff,3))

-0.663 1.387
The difference is 2.051

We could use a \(t\)-test to test whether there is a significant difference in the SAM index between these two time periods, but instead let’s use bootstrapping.

So, we are going to randomly resample our original SAM time series and grab a pair of 33-year samples (with replacement) and compute the difference.

N=33
SAM_diff_bootstrap=[]
for i in np.arange(0,100000):
    SAM1_tmp = np.mean(np.random.choice(X,N))
    SAM2_tmp = np.mean(np.random.choice(X,N))
    SAM_diff_bootstrap.append(SAM2_tmp - SAM1_tmp)

Now, we can plot the distribution of differences to examine the probability of obtaining our observed difference. First, we generate the histogram…

# create bins 
bins = np.linspace(np.min(SAM_diff_bootstrap),np.max(SAM_diff_bootstrap),int((np.max(SAM_diff_bootstrap)-np.min(SAM_diff_bootstrap))/0.01))

# calculate the histogram
histSAM_diff,bins = np.histogram(SAM_diff_bootstrap,bins)

# convert counts to frequency
freqSAM_diff = histSAM_diff/len(SAM_diff_bootstrap)

… next, we plot it along with our observed SAM difference.

plt.figure(figsize=(10,6))

# xbins for plotting
xbins = np.linspace(np.min(SAM_diff_bootstrap),np.max(SAM_diff_bootstrap),len(histSAM_diff))

# plot the distribution
plt.plot(xbins,freqSAM_diff,'deepskyblue',linewidth=2)
plt.axvline(np.mean(SAM_diff),color='gray',label="Observed SAM Difference")
plt.ylabel('Frequency')
plt.xlabel('SAM Difference')
plt.legend(loc='upper left')
plt.ylim(0,0.015)
plt.title('PDF of Bootstrapped Sample SAM Differences')
plt.tight_layout()

Now, we can clearly visualize that our observed difference is in the very tail of our bootstrapped sampling distribution, indicating that the probability of obtaining such a large difference is very low. How low? Let’s see.

# use st.percentileofscore() to find the probability of obtaining a value of NAO_diff or higher
1-st.percentileofscore(SAM_diff_bootstrap,np.mean(SAM_diff))/100.0

2.999999999997449e-05

Really low! This means that the difference between these two 33-year segments is quite anomalous. Note that here we are assuming that each year in our 33-year segments is independent, but in reality this is often not the case. We will discuss how to deal with this issue later in the course.

Jackknifing

Jackknife resampling is very similar to bootstrapping except that you systematically remove one value from your sample, and calculate the statistic, then put the value back into the sample and remove the next value, calculate the statistics…and on and on.

Let’s take a look at an example of this using the SAM index data from above.

Suppose we are interested in examining the trend in this data (we will talk about trends and regression analysis in more detail next week).

We can use a simple linear fit to estimate the trend.

# Compute Linear fit
trend = np.polyfit(TIME_SAM,X,1) #trend has two components: the slope and the intercept 

# Plot data with linear fit
plt.figure(figsize=(10,6))
plt.plot(TIME_SAM,X,color = 'black', linewidth = 1.5)
plt.plot(TIME_SAM,TIME_SAM*trend[0]+trend[1],'--',color = 'black', linewidth = 1.5)
plt.xlabel('Year')
plt.ylabel('SAM Index')
plt.title('Estimation of SAM best-fit with regression')
plt.ylim(-5,5);
plt.xlim(min(TIME_SAM), max(TIME_SAM));
plt.axhline(0,color='gray')

<matplotlib.lines.Line2D at 0x178253ca0>

How do we determine the 95% confidence levels on this trend line?

Next week, we will talk about how to do this using \(z\)- or \(t\)-statistics, but here we will use the jackknife resampling method to get at an estimate.

Following the jackknife method, we will remove one data point from our time series at a time, recalculate the linear fit and repeat. We will get a distribution of possible trends from which we can find the 95% confidence bounds.

# Initialize an array of size (69,2). The dimension of size 69 reflects the number of iterations
# and the dimension of size 2 reflects that we will store our slopes and intercepts.
M = np.zeros((len(X),2))

# Initialize plot
plt.figure(figsize=(10,6))
plt.xlabel('Year');
plt.ylabel('SAM Index')
plt.title('Estimation of SAM best-fit with regression')

# Loop over data points
for j, val in enumerate(X):
    
    # Remove one data point
    Xj = np.delete(X,j)
    Tj = np.delete(TIME_SAM,j)
    
    # Calculate the linear fit using np.polyfit()
    trendj = np.polyfit(Tj,Xj,1) 
    
    # Save the slope and the intercept in an array
    M[j,0] = trendj[0] #slope
    M[j,1] = trendj[1] #intercept
    
    # Plot data
    if j < len(X)-1:
        
        # Plot the data point we are removing
        plt.plot(TIME_SAM[j],val,'.',color = 'mediumorchid', markersize = 15)
        
        # Plot the trend line after the data point is removed
        plt.plot(Tj,Tj*trendj[0] + trendj[1],'--', color = np.random.random_sample(size = 3))
        plt.ylim(-5,5);
        display.clear_output(wait=True)
        display.display(plt.gcf())
        time.sleep(0.05)
    else:
        plt.plot(TIME_SAM[j],val,'.',color = 'mediumorchid', markersize = 15)
        plt.plot(Tj,Tj*trendj[0] + trendj[1],'--', color = np.random.random_sample(size = 3))
        display.clear_output(wait=True)

As you can see we get a slightly different slope and y-intercept each time we remove a data point.

We can now use stats.scoreofpercentile() to find the confidence bounds on the slope and y-intercept.

slopes = np.squeeze(M[:,0])
yints = np.squeeze(M[:,1])

CI_upper_slope = st.scoreatpercentile(slopes,97.5)
CI_lower_slope = st.scoreatpercentile(slopes,2.5)
print("CI's for slope:",CI_upper_slope, CI_lower_slope)

CI_upper_yint = st.scoreatpercentile(yints,97.5)
CI_lower_yint = st.scoreatpercentile(yints,2.5)
print("CI's for y-intercept:",CI_upper_yint, CI_lower_yint)

CI's for slope: 0.05010386823819755 0.04446816245566715
CI's for y-intercept: -88.18072497699862 -99.41553803176868

To get a sense of the entire ditribution of slopes and y-intercepts, let’s plot the histograms and add our confidence intervals as vertical lines.

# Plot histogram of slopes
plt.figure(figsize=(18,6))
plt.subplot(1,2,1)
xint = np.arange(.040,.055,.00025)
b, bin_edges = np.histogram(slopes,xint)
plt.plot(bin_edges[:-1],b/M.shape[0],color='deepskyblue')

# Plot 95% confidence intervals on slope
plt.axvline(CI_upper_slope,color='mediumorchid')
plt.axvline(CI_lower_slope,color='mediumorchid')

plt.xlabel('Slope')
plt.ylabel('Frequency')
plt.title('Distribution of SAM slopes from jackknife')

# Plot histogram of intercepts
plt.subplot(1,2,2)
xint = np.arange(-102.5,-85.,.5)
y, bin_edges = np.histogram(yints,xint)
plt.plot(bin_edges[:-1],y/M.shape[0],color='deepskyblue')

# Plot 95% confidence intervals on slope
plt.axvline(CI_upper_yint,color='mediumorchid')
plt.axvline(CI_lower_yint,color='mediumorchid')

plt.xlabel('y-intercept')
plt.ylabel('Frequency')
plt.title('Distribution of SAM y-intercept from jackknife')
plt.tight_layout()