Challenge: Resampling Approach to Compare Mean Values of the Datasets

We can also use the resampling approach to test the hypothesis with non-Gaussian datasets. Resampling is a technique for sampling from an available data set to generate additional samples, each of which is considered representative of the underlying population.

Approach description

Let's describe the most simple resampling method to check main hypothesis that two datasets X and Y have equal mean values:

Concatenate both arrays (X and Y) into one big array;
Shuffle that entire array so observations from each group are spread randomly throughout that array instead of being separated at the breaking point;
Arbitrarily split the array in the breaking point (X_length), assign observations below index len(X_length) to Group A and the rest to Group B;
Subtract the mean of this new Group A from the mean of the new Group B. This would give us one permutation test statistic;
Repeat those steps N times to simulate the main hypothesis distribution;
Calculate test statistics on initial sets X and Y;
Determine critical values of the main hypothesis distribution;
Check if the test statistic calculated on initial sets falls into a critical area of the main hypothesis distribution. If it falls then reject the main hypothesis.

Let's apply this approach in code:

Task

Swipe to start coding

Your task is to implement described above resampling algorithm and to check the corresponding hypothesis on two datasets:

Use the np.concatenate() method to merge X and Y arrays.
Use the .shuffle() method of the np.random module to shuffle data in the merged array.
Use np.quantile() method to calculate left critical value.
Use the created resampling_test() function to check the hypothesis on generated data.

Solution

import numpy as np

import matplotlib.pyplot as plt

from scipy.stats import chi2, expon

def resampling_test(X, Y, N, alpha):

# Concatenate both arrays into one big array

data = np.concatenate((X, Y))

X_length = len(X)

# Calculate the observed test statistic on initial sets X and Y

obs_stat = np.mean(X) - np.mean(Y)

# Create an array to store permutation test statistics

perm_stats = np.zeros(N)

# Perform N permutations

for i in range(N):

# Shuffle the entire array

np.random.shuffle(data)

# Split the shuffled array into two groups

perm_X = data[:X_length]

perm_Y = data[X_length:]

# Calculate the test statistic for this permutation

perm_stat = np.mean(perm_X) - np.mean(perm_Y)

perm_stats[i] = perm_stat

# Calculate the two-tailed critical values of the main hypothesis distribution

crit_l = np.quantile(perm_stats, alpha / 2)

crit_r = np.quantile(perm_stats, 1 - alpha / 2)

# Plot the distribution of permutation test statistics and critical areas

plt.hist(perm_stats, bins=30)

plt.axvline(x=crit_l, color='r', linestyle='--', label='Critical values')

plt.axvline(x=crit_r, color='r', linestyle='--')

plt.axvline(x=obs_stat, color='g', linestyle='-', label='Observed statistic')

Everything was clear?

Thanks for your feedback!

Section 4. Chapter 5