<jupyter_start><jupyter_text>A Bootstrap Method for One Sample Hypothesis Testing A one-sample test determines whether or not a sample statistic is statistically different from a population value. We will now use a bootstrap method with the same purpose.A bootstrap hypothesis test starts with a test statistic `t(x)` (not necessarily an estimate of a parameter).<jupyter_code>import numpy as np
import pandas as pd
from scipy import stats as st
import matplotlib.pyplot as plt
import seaborn as sns;
sns.set_style("whitegrid")  
from statsmodels.distributions.empirical_distribution import ECDF<jupyter_output><empty_output><jupyter_text>Generating data Let's start generating ficticious data.<jupyter_code># Generating data
np.random.seed(1234)
data = np.random.randint(158, 175, 50)
sample_mean = np.mean(data)
print(sample_mean)
sns.histplot(x=data)
plt.axvline(x=np.mean(sample_mean), color='black');<jupyter_output><empty_output><jupyter_text>Our data does not follow a normal distribution.  Calculating the sample distribution using bootstrap Let us define `generateSamples` function to generate samples WITH REPLACEMENT using `data`.<jupyter_code>def generateSamples(sample_data, num_samples=10000, sample_size=100):
    '''
    It returns a DataFrame where each column is a sample.
    '''
    col = ['S'+str(k) for k in range(num_samples)]
    samp = np.random.choice(sample_data, replace=True, size=sample_size*num_samples)
    df = pd.DataFrame(samp.reshape(sample_size, num_samples), columns=col)
    return(df)<jupyter_output><empty_output><jupyter_text>Let us generate 10,000 samples, and save them into a DataFrame.<jupyter_code>dfSamples = generateSamples(data, num_samples=10000, sample_size=40)
print(dfSamples.shape)
dfSamples.head()<jupyter_output>(40, 10000)<jupyter_text>We are going to use the mean as our statistic test.`sample_distribution_mean` is a list with the mean of all samples.<jupyter_code>sample_distribution_mean = dfSamples.mean()<jupyter_output><empty_output><jupyter_text>This is the histogram of `sample_distribution_mean`:<jupyter_code># Plotting a histogram of sample_distribution_mean
sns.histplot(x=sample_distribution_mean, color='mediumseagreen')
plt.axvline(x=np.mean(sample_distribution_mean), color='black')
plt.title('Sample Distribution Mean = %.2f' %(np.mean(sample_distribution_mean)));
# kde plot of sample_distribution_mean
sns.kdeplot(x=sample_distribution_mean, color='mediumseagreen', shade=True)
plt.axvline(x=np.mean(sample_distribution_mean), color='black')
plt.title('Sample Distribution Mean = %.2f' %(np.mean(sample_distribution_mean)));<jupyter_output><empty_output><jupyter_text>The hypothesis Suppose we want to test the hypothesis:- `Ho: mean   = 170`- `Ha: mean != 170`The first step is to construct the sample distribution under Ho; we need to build a sample distribution with a mean of `170`.<jupyter_code># sample distribution under Ho
sample_distribution_Ho = sample_distribution_mean - np.mean(sample_distribution_mean) + 170
sns.kdeplot(x=sample_distribution_Ho, color='dodgerblue', shade=True)
plt.axvline(x=np.mean(sample_distribution_Ho), color='black')
plt.title('Sample Distribution Mean = %.2f' %(np.mean(sample_distribution_Ho)));<jupyter_output><empty_output><jupyter_text>Critical Regions We have three critical regions, according to the hypothesis we want to test:- left- right- two-sidedTo calculate the p-values, we will use two variables:- `sample_distribution_Ho`: sample distribution under Ho- `test_value`: hypothesized valueLet's define alpha = 0.05<jupyter_code>alpha = 0.05
plt.figure(figsize=(14,3))

plt.subplot(131)
ax = sns.kdeplot(x=sample_distribution_Ho, color='lightskyblue', shade=True, alpha=0.4)
cv1 = np.percentile(sample_distribution_Ho, alpha*100)              # left critical value   
plt.axvline(x = cv1, ymin=0, ymax=0.5, color='orangered', linewidth=2)
plt.text(cv1-2, 0.05, 'alpha=0.05', color='red', size='13')
plt.title('Critical Region: left')

plt.subplot(132)
ax = sns.kdeplot(x=sample_distribution_Ho, color='lightskyblue', shade=True, alpha=0.4)
cv2 = np.percentile(sample_distribution_Ho, 100-alpha*100)          # right critical value 
plt.axvline(x = cv2, ymin=0, ymax=0.5, color='orangered', linewidth=2)
plt.text(cv2+0.1, 0.05, 'alpha=0.05', color='red', size='13')
plt.title('Critical Region: right')

plt.subplot(133)
ax = sns.kdeplot(x=sample_distribution_Ho, color='lightskyblue', shade=True, alpha=0.4)
cv1 = np.percentile(sample_distribution_Ho, (alpha/2)*100)          # left critical value 
cv2 = np.percentile(sample_distribution_Ho, 100-alpha*100)          # right critical value      
plt.axvline(x = cv1, ymin=0, ymax=0.5, color='orangered', linewidth=2)
plt.axvline(x = cv2, ymin=0, ymax=0.5, color='orangered', linewidth=2)
plt.text(cv1-2.2, 0.05, 'alpha=0.025', color='red', size='13')
plt.text(cv2 + 0.1, 0.05, 'alpha=0.025', color='red', size='13')
plt.title('Critical Region: two-sided')

plt.show();<jupyter_output><empty_output><jupyter_text>Calculating the p-value The p-value is the probability of getting a test statistic more extreme than the one we have calculated from our sample, given Ho is true.We reject Ho if our p-value is small enough.For calculating the p-value we will:- Calculate the number of simulated test stats that are more extreme than our hypothesized `test_value`- Calculate it using the empirical distribution function (ECDF)Both options are equivalent. They get the same result.Let us start defining the `test_value`:<jupyter_code>test_value = 170<jupyter_output><empty_output><jupyter_text>p-value for alternative hypothesis: `smaller` (critical region: left) - `Ho: measure = test_value`- `Ha: measure <  test_value`In our case, we will be testing: - `Ho: mean = 170`- `Ha: mean < 170`<jupyter_code>plt.figure(figsize=(5,3))
ax = sns.kdeplot(x=sample_distribution_Ho, color='lightskyblue', shade=True, alpha=0.4)
plt.plot(sample_mean, 0.01, 'o', color='black')
cv1 = np.percentile(sample_distribution_Ho, alpha*100)              # left critical value 
plt.axvline(x = cv1, ymin=0, ymax=0.5, color='orangered', linewidth=2)
plt.text(cv1-2, 0.05, 'alpha=0.05', color='red', size='13')
plt.title('Critical Region: left');<jupyter_output><empty_output><jupyter_text>Let us calculate the number of simulated test stats that are less than or equal to our `sample_mean`.p-value is the area under the curve from -infinity to the `sample_mean`<jupyter_code>lowtail = np.sum(np.array(sample_distribution_Ho)<=sample_mean)
pval_s = lowtail/(len(sample_distribution_Ho))
print('Smaller alternative hypothesis: p-value = %.4f' %pval_s)<jupyter_output>Smaller alternative hypothesis: p-value = 0.0006<jupyter_text>Let us calculate the p-value using the empirical distribution.<jupyter_code># Creating the empirical distribution
ecdf = ECDF(sample_distribution_Ho)
pval_s = ecdf(sample_mean)
print('Smaller alternative hypothesis: p-value = %.4f' %pval_s)<jupyter_output>Smaller alternative hypothesis: p-value = 0.0006<jupyter_text>As you can see, both values are the same.  p-value for alternative hypothesis: `larger` (critical region: right) - `Ho: measure = test_value`- `Ha: measure > test_value`In our case, we will be testing: - `Ho: mean = 170`- `Ha: mean > 170`<jupyter_code>plt.figure(figsize=(5,3))
ax = sns.kdeplot(x=sample_distribution_Ho, color='lightskyblue', shade=True, alpha=0.4)
plt.plot(sample_mean, 0.01, 'o', color='black')
cv2 = np.percentile(sample_distribution_Ho, 100-alpha*100)      # right critical value
plt.axvline(x = cv2, ymin=0, ymax=0.5, color='orangered', linewidth=2)
plt.text(cv2+0.1, 0.05, 'alpha=0.05', color='red', size='13')
plt.title('Critical Region: right');<jupyter_output><empty_output><jupyter_text>p-value here is 1 - the area under the curve (pfd) from the `sample_mean` to infinity because it is a right critical region.Let us calculate the number of simulated test stats that are greater than to our original `sample_mean`.<jupyter_code>upptail = np.sum(np.array(sample_distribution_Ho)>sample_mean)
pval_l = upptail/(len(sample_distribution_Ho))
print('Larger alternative hypothesis: p-value = %.4f' %pval_l)<jupyter_output>Larger alternative hypothesis: p-value = 0.9994<jupyter_text>Let us calculate the p-value using the empirical distribution.<jupyter_code>pval_l = 1 - ecdf(sample_mean)
print('Larger alternative hypothesis: p-value = %.4f' %pval_l)<jupyter_output>Larger alternative hypothesis: p-value = 0.9994<jupyter_text>Both ways produce the same p-value.Notice that we can add both p-values (smaller and larger) and get 1. That is because they are complementaries.<jupyter_code>print('Smaller alternative hypothesis: p-value = %.4f' %pval_s)
print('Larger alternative hypothesis: p-value  = %.4f' %pval_l)
print('Adding both values: %.1f' %(pval_s + pval_l))<jupyter_output>Smaller alternative hypothesis: p-value = 0.0006
Larger alternative hypothesis: p-value  = 0.9994
Adding both values: 1.0<jupyter_text>p-value for alternative hypothesis: `two-sided` (critical region: two-sided) - `Ho: measure  = test_value`- `Ha: measure != test_value`In our case, we will be testing: - `Ho: mean  = 170`- `Ha: mean != 170`<jupyter_code>plt.figure(figsize=(5,3))
ax = sns.kdeplot(x=sample_distribution_Ho, color='lightskyblue', shade=True, alpha=0.4)
plt.plot(sample_mean, 0.01, 'o', color='black')
cv1 = np.percentile(sample_distribution_Ho, (alpha/2)*100)     # left critical value 
cv2 = np.percentile(sample_distribution_Ho, 100-alpha*100)     # right critical value
plt.axvline(x=cv1, ymin=0, ymax=0.5, color='orangered', linewidth=2)
plt.axvline(x=cv2, ymin=0, ymax=0.5, color='orangered', linewidth=2);
plt.text(cv1-2.2, 0.05, 'alpha=0.025', color='red', size='13')
plt.text(cv2+0.1, 0.05, 'alpha=0.025', color='red', size='13')
plt.title('Alternative Hypothesis: two-sided');<jupyter_output><empty_output><jupyter_text>To calculate the p-value when the alternative hypothesis is two-sided, we need to determine whether the test_value is to the right or left in our graph, calculate the p-value in the usual way, and then duplicate it.Let's calculate the minimum number of simulated test stats that are:- less than or equal to test_vales, - greater than test_value<jupyter_code>lowtail = np.sum(np.array(sample_distribution_Ho)<=sample_mean)
upptail = np.sum(np.array(sample_distribution_Ho) >sample_mean)
pval_t = 2*np.min([lowtail, upptail])/len(sample_distribution_Ho)
print('Two-tailed test p-value = %.4f' %pval_t)<jupyter_output>Two-tailed test p-value = 0.0012<jupyter_text>Using ecdf:<jupyter_code>if sample_mean < np.mean(sample_distribution_Ho): 
    pval_t = 2*ecdf(sample_mean)
else: pval_t = 2*(1 - ecdf(sample_mean))
print('Two-tailed test p-value = %.4f' %pval_t)<jupyter_output>Two-tailed test p-value = 0.0012<jupyter_text>We do not need to use both methods to get the p-value. We can use any of them.Let us summarize the p-value calculation in one function.<jupyter_code>def getpValue(sample_distribution, obs_value, alternative='two-sided'):    
    '''
    It returns the p-value.
    sample_distribution: the sample distribution
    obs_value:           observed value 
    alpha:               significance level
    alternative:         one of the three values: 'two-sided', 'smaller', or 'larger'    
    '''
    ecdf = ECDF(sample_distribution) 
    if alternative=='two-sided':
        if obs_value < np.mean(sample_distribution): 
            p_val = 2*ecdf(obs_value)
        else: 
            p_val = 2*(1-ecdf(obs_value)) 
    elif alternative=='smaller':
        p_val = ecdf(obs_value)
    else:
        p_val = 1-ecdf(obs_value)
    return(p_val)
# Getting the p-value with an alternative hypothesis: smaller
pval = getpValue(sample_distribution_Ho, sample_mean, alternative='smaller')
print('For alternative hypothesis smaller: p_value = %.4f' % pval)
# Getting the p-value with an alternative hypothesis: larger
pval = getpValue(sample_distribution_Ho, sample_mean, alternative='larger')
print('For alternative hypothesis larger: p_value = %.4f' % pval)
# Getting the p-value with an alternative hypothesis: two-sided (default)
pval = getpValue(sample_distribution_Ho, sample_mean)
print('For alternative hypothesis two-sided: p_value = %.4f' % pval)<jupyter_output>For alternative hypothesis two-sided: p_value = 0.0012<jupyter_text>One Sample Hypothesis Test with Bootstrapping We have:- data: a representative sample- The hypothesis we want to testWe want to define a function for hypothesis testing.  Defining bootstrap one_sample_HT We will create a function for the hypothesis test. We could receive `data` as a parameter and use it to determine its mean and the sample distribution. The calculation of sample distribution is time-consuming, and we could use the same sample distribution to test several hypotheses. So, we decided to receive it as a parameter instead of calculating it in each call to the function.<jupyter_code>def boot_one_sample_HT(sample_value, sample_distribution, population_value, alpha=0.05, alternative='two-sided'):
    '''
    Bootstrap method for one-sample hypothesis test.
    sampl_value: observed value calculated from the sample
    sampl_distribution: sample distribution calculated from the sample
    pop_value: population mean
    alpha: significance level
    alternative: one of the three values: 'two-sided', 'smaller', and 'larger'
    '''
    sigHa = {'two-sided':'!=', 'smaller':'< ', 'larger':'> '}
    print('--- Bootstrap One-Sample Hypothesis Test ---')
    print('    Ho: measure  =', np.round(population_value, 2))    
    print('    Ha: measure', sigHa[alternative], np.round(population_value, 2))    
    print('    Sample Mean = %.2f' %(sample_value))  
    sample_distribution_Ho = sample_distribution - np.mean(sample_distribution) + population_value
    sns.kdeplot(x=sample_distribution_Ho, color='lightskyblue', shade=True, alpha=0.4)   # plotting the sample distribution
    plt.axvline(x=sample_value, ymin=0, ymax= 0.03, color='black', linewidth=6)          # plotting the statistic 
    plt.title('Sampling Distribution')
    p_val = getpValue(sample_distribution_Ho, sample_value, alternative)
    if alternative=='two-sided':        # plotting the graph with two critical regions
        cv1 = np.round(np.percentile(sample_distribution_Ho, (alpha/2)*100),2)           # left critical value     
        cv2 = np.round(np.percentile(sample_distribution_Ho, 100-alpha*100),2)           # right critical value
        plt.axvline(x = cv1, ymin=0, ymax=0.5, color='orangered', linewidth=2)  
        plt.axvline(x = cv2, ymin=0, ymax=0.5, color='orangered', linewidth=2)
    elif alternative=='smaller':        # plotting the graph with left critical region
        cv1 = np.round(np.percentile(sample_distribution_Ho, alpha*100),2)               # left critical value 
        plt.axvline(x = cv1, ymin=0, ymax=0.5, color='orangered', linewidth=2)
    else:                               # plotting the graph with right critical region
        cv2 = np.round(np.percentile(sample_distribution_Ho, 100-alpha*100),2)           # right critical value
        plt.axvline(x = cv2, ymin=0, ymax=0.5, color='orangered', linewidth=2)
    print('    p-value = '+str(np.round(p_val,4)))<jupyter_output><empty_output><jupyter_text>Let us use the function `boot_one_sample_HT` with `data` for the hypothesis.- `Ho: mean = 170`- `Ha: mean != 170`<jupyter_code>boot_one_sample_HT(sample_mean, sample_distribution_mean, 170)<jupyter_output>--- Bootstrap One-Sample Hypothesis Test ---
    Ho: measure  = 170
    Ha: measure != 170
    Sample Mean = 167.54
    p-value = 0.0012<jupyter_text>p-value < 0.05, we reject Ho. We conclude the mean is different from 170.Let's analyse:- `Ho: mean = 168`- `Ha: mean != 168`<jupyter_code>boot_one_sample_HT(sample_mean, sample_distribution_mean, 168)<jupyter_output>--- Bootstrap One-Sample Hypothesis Test ---
    Ho: measure  = 168
    Ha: measure != 168
    Sample Mean = 167.54
    p-value = 0.5438<jupyter_text>There are no reasons for rejecting Ho, and we cannot affirm the mean is different from 168.Let's analyse:- `Ho: mean = 166`- `Ha: mean > 166`<jupyter_code>boot_one_sample_HT(sample_mean, sample_distribution_mean, 166, alternative='larger')<jupyter_output>--- Bootstrap One-Sample Hypothesis Test ---
    Ho: measure  = 166
    Ha: measure >  166
    Sample Mean = 167.54
    p-value = 0.0186<jupyter_text>p-value < 0.05, therefore we reject Ho and conclude the mean is greater than 166.Let's compare this result with the traditional t-test.  Testing the variance Let's run the bootstrap hypothesis test for the variance.The first step is calculate the sample distribution for the variance and the variance of our sample.<jupyter_code>sample_distribution_var = dfSamples.var()
sample_var = np.var(data)
boot_one_sample_HT(sample_var, sample_distribution_var, 20)<jupyter_output>--- Bootstrap One-Sample Hypothesis Test ---
    Ho: measure  = 20
    Ha: measure != 20
    Sample Mean = 22.29
    p-value = 0.5404<jupyter_text>There are no reasons for rejecting Ho. The sample variance is not significantly different from 20.What about variance = 30? Let's formulate the hypothesis:- `Ho: var = 30`- `Ha: var < 30`<jupyter_code>boot_one_sample_HT(sample_var, sample_distribution_var, 30, alternative='smaller')<jupyter_output>--- Bootstrap One-Sample Hypothesis Test ---
    Ho: measure  = 30
    Ha: measure <  30
    Sample Mean = 22.29
    p-value = 0.0184<jupyter_text>And indeed 22.29 < 30  Testing the skew Let's use the bootstrapping method to test the hypothesis involving the skewness of data:- `Ho: skew = 0`- `Ha: skew < 0`We stat calculating the sample distribution for the skewness.<jupyter_code>sample_distribution_skew = dfSamples.skew()
from scipy.stats import skew
sample_skew = skew(data)
boot_one_sample_HT(sample_skew, sample_distribution_skew, 0, alternative='smaller')<jupyter_output>--- Bootstrap One-Sample Hypothesis Test ---
    Ho: measure  = 0
    Ha: measure <  0
    Sample Mean = -0.43
    p-value = 0.0453
{'fraction_non_alphanumeric': 0.11531411431428096, 'fraction_numerical': 0.03188357495972893, 'mean_word_length': 5.146807784226699, 'pattern_counts': {'":': 0, '<': 108, '<?xml version=': 0, '>': 101, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '27', 'dir': 'jupyter-structured-clean-dedup', 'id': '393251', 'license': 'MIT', 'n_tokens_mistral': 6455, 'n_tokens_neox': 5781, 'n_words': 1631, 'path': '04-Inferential Statistics/05-Boot_One_Samp_Test.ipynb', 'repo_name': 'Goliath-Research/Computational-Statistics', 'size': '18003'}