Exercise 4.1: Writing functions for bootstrap replicates

---

It will be useful to have some functions in your arsenal to do statistical inference with bootstrapping. In this exercise, you will write some handy functions. You should test these functions out on real data.

a) In the lessons, we wrote a function, draw_bs_rep() to draw a single bootstrap replicate out of a single set of repeated measurements. Update this function to have a size keyword argument so that you can draw many bootstrap replicates and return a Numpy array of the replicates. Here are step-by-step instructions.

1. Define a function with call signature draw_bs_reps(data, func, rg, size=1, args=()), where func is a function that takes in an array and returns a statistic; it has call signature func(data, *args). Examples that could be passed in as func are np.mean, np.std, np.median, or a user-defined function. rg is an instance of a Numpy random number generator. size is the number of replicates to generate.

2. Write a good doc string.

3. Define n to be the length of the input data array.

4. Use a list comprehension to compute a list of bootstrap replicates.

5. Return the replicates as a Numpy array.

b) Write a function analogous to the one in part (a) except for pairs bootstrap. The call signature should be draw_bs_pairs(data1, data2, func, rg, size=1, args=()), where func has call signature func(data1, data2, *args).

You will want to include these in a module so you can use it over and over again. I will not be providing this functionality in the bootcamp_utils module; I want you to write this yourself. (Or, you can install the dc_stat_think module that I wrote using pip, which has this and many other useful functions for bootstrapping.)

Solution

import numpy as np
import pandas as pd

a) I show the function below.

def draw_bs_reps(data, func, rg, size=1, args=()):
    """
    Generate bootstrap replicates out of `data` using `func`.

    Parameters
    ----------
    data : array_like
        One-dimensional array of data.
    func : function
        Function, with call signature `func(data, *args)` to compute
        replicate statistic from resampled `data`.
    rg : random number generator instance
        Either `np.random` or the result of `np.random.default_rng()`.
    size : int, default 1
        Number of bootstrap replicates to generate.
    args : tuple, default ()
        Arguments to be passed to `func`.

    Returns
    -------
    output : ndarray
        Bootstrap replicates computed from `data` using `func`.
    """
    return np.array(
        [
            func(rg.choice(data, replace=True, size=len(data)), *args)
            for _ in range(size)
        ]
    )

Let’s try this function out on the beak depth data from 1975 to get bootstrap replicates of the mean.

rg = np.random.default_rng()

df = pd.read_csv('data/grant_complete.csv', comment='#')
inds = (df['year'] == 1975) & (df['species'] == 'scandens')
bd_1975 = df.loc[inds, 'beak depth (mm)'].values

# Compute replicates
bs_reps = draw_bs_reps(bd_1975, np.mean, rg, size=2000)

# 95% confidence interval
print(np.percentile(bs_reps, [2.5, 97.5]))

[8.8478046  9.07747701]

Nice!

b) I show the function below.

def draw_bs_pairs(data1, data2, func, rg, size=1, args=()):
    """
    Perform pairs bootstrap for single statistic.

    Parameters
    ----------
    data1 : array_like
        First data set. Must be same length as `data2`.
    data2 : array_like
        Second data set. Must be same length as `data1`.
    func : function
        Function, with call signature `func(data1, data2, *args)` to
        compute replicate statistic from pairs bootstrap sample. It
        must return a single, scalar value.
    rg : random number generator instance
        Either `np.random` or the result of `np.random.default_rng()`.
    size : int, default 1
        Number of pairs bootstrap replicates to draw.
    args : tuple, default ()
        Arguments to be passed to `func`.

    Returns
    -------
    output : ndarray
        Bootstrap replicates.
    """
    n = len(data1)

    # Set up array of indices to sample from
    inds = np.arange(n)

    # Initialize replicates
    bs_replicates = np.empty(size)

    # Generate replicates
    for i in range(size):
        bs_inds = rg.choice(inds, n)
        bs_1, bs_2 = data1[bs_inds], data2[bs_inds]
        bs_replicates[i] = func(bs_1, bs_2, *args)

    return bs_replicates

Let’s take it for a spin with the correlation coefficient between the length and depth of beaks in 1975.

bl_1975 = df.loc[inds, 'beak length (mm)'].values

def corrcoef(data1, data2):
    cov = np.cov(data1, data2)
    return cov[0, 1] / np.sqrt(cov[0, 0] * cov[1, 1])


# Compute replicates
bs_reps = draw_bs_pairs(bd_1975, bl_1975, corrcoef, rg, size=2000)

# 95% confidence interval
print(np.percentile(bs_reps, [2.5, 97.5]))

[0.45686897 0.75179156]

Looks good!

Computing environment

%load_ext watermark
%watermark -v -p numpy,pandas,jupyterlab

Python implementation: CPython
Python version       : 3.9.12
IPython version      : 8.3.0

numpy     : 1.21.5
pandas    : 1.4.2
jupyterlab: 3.3.2