<jupyter_start><jupyter_text>Hypothesis testing> How to practically compare aspects of two distributions<jupyter_code>#hide
%matplotlib inline
%load_ext autoreload
%autoreload 2
#export
import matplotlib.pyplot as plt
import matplotlib as mpl
import numpy as np
from scipy import stats
import pandas as pd

from thinkstats import basic_stats
from thinkstats import nsfg
#hide
mpl.style.use("bmh")<jupyter_output><empty_output><jupyter_text>**Basic steps for the investigation**1. calc *test statistic*: e.g. coin of interest is biased ($p \ne 0.5$)2. define *null hypothesis*: e.g. if the coin were unbiased / fair we would have $p = 0.5$)3. compute the *p-value*The core of the p-value is to estimate how probable an aspect of an observed dataset is. The aspect is the *test statistic*, like the mean occurrences of heads in a series of coin tosses. To retrieve the p-value we compute the test statistic for the observed dataset and for a lot of perturbed version of the dataset. Then we can quantify how rare the value for the test statistic we calculated for the observed, unperturbed, dataset really is.What I find needs a bit of getting used to is to figure out how to perturb the dataset properly. To apply this to various practical situations one may need quite a bit of experience. In the first example below we want to verify if two coins deviate from oneanother in their fairness. So our test statistic is the difference of the share of obtained heads. We define our hypothesis as: the two coins are equally fair. To verify this we perturb our dataset of the coin tosses of both coins by repeatedly, `n_test` times, mixing their results up and computing the test statistic. We do this because of both coins are equally fair the mixing up shouldn't matter. This way we can gauge where in the distribution of the test statistic for the mixed up datasets the test statistic of the original dataset lies. If we want to know how much more biased towards heads coin 1 is than coin 0 we only need to compute the cumulative density function over the test statistics from the mixups and find the probability (left tail only) for the test statistic of the original dataset. If we want to know about general deviation in fairness we need to consider both sides of the distribution of the test statistics from the mixups.  Synthetic data  Discrete - Bernoulli $p$ Defining generative distributions for population 0 and 1<jupyter_code>pop0_gen = stats.bernoulli(p=.5)
pop1_gen = stats.bernoulli(p=.7)<jupyter_output><empty_output><jupyter_text>Number of samples to generate for each population<jupyter_code>n0 = 100
n1 = 120<jupyter_output><empty_output><jupyter_text>Generating the populations<jupyter_code>samples0 = pop0_gen.rvs(size=n0)
samples1 = pop1_gen.rvs(size=n0)
samples0[:5], samples1[:5]<jupyter_output><empty_output><jupyter_text>We claim (hypothesis) the two coins `pop0_gen` and `pop1_gen` do not deviate from each other in their fairness → test statistic is the difference of their means<jupyter_code>#export
def get_mean_test_statistic(a:np.ndarray, b: np.ndarray) -> float: return a.mean() - b.mean()
test_statistic = get_mean_test_statistic(samples0, samples1); test_statistic<jupyter_output><empty_output><jupyter_text>To compute the p-value we need to compute the distribution of bootstrapped distributions and their test statistics. `n_tests` defines the number of bootstrappings to perform<jupyter_code>n_tests = 1000
# export
def get_boot_test_statistics(a:np.ndarray, b: np.ndarray, n_tests:int=100, 
                             test_statistic_fun:callable=get_mean_test_statistic,
                             mode="merge+shuffle"):
    "Computes bootstrapped test statistics using shuffling"
    test_statistic_boot = np.zeros(n_tests)
    if mode == "merge+shuffle":
        samples = np.hstack((a, b))
        n = len(a)
        for i in range(n_tests):
            np.random.shuffle(samples)
            _a, _b = samples[:n], samples[n:]
            test_statistic_boot[i] = test_statistic_fun(_a, _b)
    elif mode == "shuffle":
        for i in range(n_tests):
            _a, _b = np.random.choice(a, size=len(a)), np.random.choice(b, size=len(b))
            test_statistic_boot[i] = test_statistic_fun(_a, _b)
    else:
        raise NotImplementedError("`mode` needs to be 'merge+shuffle' or 'shuffle'.")
    return test_statistic_boot
%%time
test_statistic_boot = get_boot_test_statistics(samples0, samples1, n_tests=n_tests)
# export
def get_p_value(x:np.ndarray, test_statistic:float):
    "Computes the p-value using `x` as generated from `get_boot_test_statistics` and the original `test_statistic`"
    kde = stats.gaussian_kde(x)
    p_val = kde.integrate_box_1d(-np.inf, test_statistic)
    if p_val > .5:
        p_val = 1. - p_val
    return p_val
%%time
p_val = get_p_value(test_statistic_boot, test_statistic); p_val<jupyter_output>Wall time: 975 µs<jupyter_text>Visualizing where the test statistic of the observation lies relative to the bootstraps for the evaluation of the significance<jupyter_code>#export
def plot_hist_with_kde_for_hypothesis_test(test_statistic_boot:np.ndarray, p_val:float, 
                                           left_ext:float=.5, right_ext:float=.5):
    ax = basic_stats.plot_hist_with_kde(test_statistic_boot, 
                                   plot_params={"xlabel": "test statistic", "ylabel": "frequency", 
                                                "title": f"Sampled test statistics: p-value = {p_val:.4f}"},
                                   fig_params={"figsize": (7,7)}, left_ext=left_ext, right_ext=right_ext)
    ax.plot([test_statistic, test_statistic], ax.get_ylim(), marker="o", color="red", label="actual test statistic")
    ax.legend(loc="best")
    return ax
ax = plot_hist_with_kde_for_hypothesis_test(test_statistic_boot, p_val, left_ext=.5, right_ext=.5)
plt.show()<jupyter_output><empty_output><jupyter_text>Discrete observables - $n$-sided dice $p$s Comparing fairness of dice using the $\chi^2$ test statistic<jupyter_code>pop0_gen = stats.rv_discrete(values=(np.arange(5), [.2, .2, .2, .2, .2]))
pop1_gen = stats.rv_discrete(values=(np.arange(5), [.2, .25, .15, .2, .2]))

n0 = 1000
n1 = 120

samples0 = pop0_gen.rvs(size=n0)
samples1 = pop1_gen.rvs(size=n0)
samples0[:5], samples1[:5]
pd.Series(samples0).value_counts().sort_index()
#export
def get_chi2_test_statistic(a:np.ndarray, b: np.ndarray) -> float: 
    "`b` is the reference"
    freq_e = pd.Series(b).value_counts(normalize=True).sort_index()
    freq_o = pd.Series(a).value_counts(normalize=True).sort_index()
    return ((freq_o - freq_e)**2/freq_e).sum()
test_statistic = get_chi2_test_statistic(samples1, samples0); test_statistic
#hide
n_tests = 1000
%%time
test_statistic_boot = get_boot_test_statistics(samples1, samples0, n_tests=n_tests, 
                                               test_statistic_fun=get_chi2_test_statistic)
%%time
p_val = get_p_value(test_statistic_boot, test_statistic); p_val
ax = plot_hist_with_kde_for_hypothesis_test(test_statistic_boot, p_val, left_ext=.5, right_ext=.5)
plt.show()<jupyter_output><empty_output><jupyter_text>Continuous - correlation Testing how significant the ranked Spearman correlation of two variables is<jupyter_code>x_noise = stats.norm(loc=0, scale=1.)
y_noise = stats.norm(loc=0, scale=1.)

n = 100

samples0 = np.linspace(-1,1,n) + x_noise.rvs(size=n)
samples1 = samples0 + y_noise.rvs(size=n)
samples0[:5], samples1[:5]
plt.scatter(samples0, samples1)
plt.show()
# export
def get_spearmanr_test_statistic(a:np.ndarray, b:np.ndarray) -> float:
    "Computes the ranked Spearman correlation"
    assert len(a) == len(b)
    return stats.spearmanr(a, b).correlation
test_statistic = get_spearmanr_test_statistic(samples1, samples0); test_statistic
#hide
n_tests = 1000
%%time
test_statistic_boot = get_boot_test_statistics(samples1, samples0, n_tests=n_tests, 
                                               test_statistic_fun=get_spearmanr_test_statistic)
%%time
p_val = get_p_value(test_statistic_boot, test_statistic); p_val
ax = plot_hist_with_kde_for_hypothesis_test(test_statistic_boot, p_val, left_ext=.5, right_ext=.5)
plt.show()<jupyter_output><empty_output><jupyter_text>NSFG data Computing if the weight of newborns to first time mothers is significantly different to mothers with 1+ children using `get_mean_test_statistic`.<jupyter_code>df = pd.read_csv(nsfg.clean_file_path)
df.head()<jupyter_output><empty_output><jupyter_text>Grouping samples into first newborn (`first`) and 2nd or higher baby count (`other`)<jupyter_code>live = df.loc[df["outcome"]==1, :]
first, other = live.loc[live["pregordr"]==1, :], live.loc[live["pregordr"]!=1, :]
len(first), len(other)<jupyter_output><empty_output><jupyter_text>Collecting the weights and cleaning<jupyter_code>col = "totalwgt_kg"  # prglngth = pregnancy length, totalwgt_kg = total weight of the newborn
samples0, samples1 = first[col].dropna().values, other[col].dropna().values<jupyter_output>[3.9972795 4.139027  3.8838815 3.4302895 3.5436875] [3.572037  3.175144  2.8066005 4.3374735 3.798833 ]<jupyter_text>Computing the test statistic over the original dataset<jupyter_code>test_statistic = get_mean_test_statistic(samples0, samples1); test_statistic<jupyter_output><empty_output><jupyter_text>Mixing up the dataset computing test statistics<jupyter_code>%%time
test_statistic_boot = get_boot_test_statistics(samples0, samples1, n_tests=n_tests)
%%time
p_val = get_p_value(test_statistic_boot, test_statistic); p_val
ax = plot_hist_with_kde_for_hypothesis_test(test_statistic_boot, p_val, left_ext=.5, right_ext=.5)
plt.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.09585921325051759, 'fraction_numerical': 0.02360248447204969, 'mean_word_length': 4.85160508782556, 'pattern_counts': {'":': 6, '<': 55, '<?xml version=': 0, '>': 60, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 10, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '14', 'dir': 'jupyter-structured-clean-dedup', 'id': '457765', 'license': 'Apache-2.0', 'n_tokens_mistral': 3354, 'n_tokens_neox': 3166, 'n_words': 885, 'path': 'notebooks/02_hypothesis_testing.ipynb', 'repo_name': 'eschmidt42/thinkstats', 'size': '9660'}