<jupyter_start><jupyter_text>Inferential Statistics - Traditional x Bootstraping

I'm using this notebok to better understand inferential statistics, specially the differences between traditional and bootstraping approaches for confidence interval (CI) estimation and hypothesis testing.

A simple well-know dataset is used in this analysis. This dataset will be taken as our population, and a single Simple Random Sample (SRS) is going to be used.  Data preparation

We're not going to do a full EDA, since our goal is just to get a better undertanding on inferential statistcs.<jupyter_code># load the libraries
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
import statsmodels.stats.api as sms
from scipy import stats
# load the titanic dataset
df = sns.load_dataset('titanic')
# check the shape
df.shape
# check the datatypes
df.dtypes
# check for null values
df.isna().sum()
# check the first 5 rows
df.head()
# Take a SRS from the population
sample_size = 100
sample_df = df.sample(sample_size, replace=True)<jupyter_output><empty_output><jupyter_text>What is the 95% CI for the 'fare' price<jupyter_code># Compute the mean, std for the 'fare' from our population
# This is NOT going to be used in inference, it is just for comparision
pop_mean, pop_std = df['fare'].mean(), df['fare'].std()
pop_mean, pop_std<jupyter_output><empty_output><jupyter_text>Traditional Approach<jupyter_code># calculate the sample_mean and sample_standard_deviation
sample_mean = sample_df['fare'].mean()
standard_error = sample_df['fare'].std()/np.sqrt(sample_size)

sample_mean, standard_error
# the critical value for 95% CI using the Z-Dist
stats.norm.ppf(0.025), stats.norm.ppf(.975)
# the critical value for 95% CI using the T-Dist
stats.t.ppf(0.025, df=sample_size-1), stats.t.ppf(.975, df=sample_size-1)<jupyter_output><empty_output><jupyter_text>Since we're inferencing the population mean, and our sample size is 'small', we can't use the z-distribution, and have to use the t-distribution<jupyter_code># Calculate the critical value for 95% CI (between 2.5% and 97.5%)
critical_value = stats.t.ppf(.975, df=sample_size-1)

critical_value
# Calculate the 95% CI
ci_lower = sample_mean-standard_error*critical_value
ci_upper = sample_mean+standard_error*critical_value

ci_lower, ci_upper
# Print the results
print("Statistic:", sample_mean)
print("Confidence Interval:", ci_lower, ci_upper)
print("Margin of Error: +-", critical_value*standard_error)
# 95% CI using the statsmodels
sms.DescrStatsW(sample_df['fare']).tconfint_mean(alpha=0.05)
# 95% CI using the scipy
stats.t.interval(0.95, df=sample_size-1 ,loc=sample_mean, scale=standard_error)
# Plot the sample distribution 
plt.hist(np.random.normal(sample_mean, standard_error, 10000), bins=20, alpha=.2)
plt.axvline(x = pop_mean, color = "black", linestyle = "dotted", label='Population Mean')
plt.axvline(x = sample_mean, color = "red", linestyle = "dotted", label='Sample Mean')
plt.axvline(x = ci_lower, color = "pink", linestyle = "dotted", label='CI 95%')
plt.axvline(x = ci_upper, color = "pink", linestyle = "dotted")
plt.legend()
plt.show()<jupyter_output><empty_output><jupyter_text>Bootstraping Approach
1. Take the sample
2. Take n bootstrap samples from the sample
3. Compute the mean for each bootstrap sample
4. Compute the percentiles (.025 and .0975) of the bootstrap sample means<jupyter_code># take n bootstrap samples from the sample
bootstrap_means_list = []
for _ in range(10000):
  # take one bootstrap sample
  bootstrap_sample = sample_df.sample(sample_size, replace=True)
  # compute the statistic for each sample
  bootstrap_mean = bootstrap_sample['fare'].mean()
  # add to the list
  bootstrap_means_list.append(bootstrap_mean)

# use the list to generate a DataFrame (for ease of use)
bootstrap_means_df = pd.DataFrame({"Bootstrap":bootstrap_means_list})
# Calculate the sample distribution mean
bootstrap_mean = bootstrap_means_df['Bootstrap'].mean()
bootstrap_mean
# Calculate the 95% CI using the percentiles from the bootstrap sample means
bootstrap_ci_lower = bootstrap_means_df['Bootstrap'].quantile(.025)
bootstrap_ci_upper = bootstrap_means_df['Bootstrap'].quantile(.975)
bootstrap_ci_lower, bootstrap_ci_upper
# Print the results
print("Statistic:", bootstrap_mean)
print("Confidence Interval:", bootstrap_ci_lower, bootstrap_ci_upper)
print("Margin of Error: TODO") # What is the right way to calculate the Margin of Error (since quantiles are not symetric)
# plot the distribution from the bootstrap means 
bootstrap_means_df.plot(kind='hist', bins=20, alpha=0.2)
plt.axvline(x = pop_mean, color = "black", linestyle = "dotted", label='Population Mean')
plt.axvline(x = bootstrap_mean, color = "red", linestyle = "dotted", label='Sample Mean')
plt.axvline(x = bootstrap_ci_lower, color = "pink", linestyle = "dotted", label='CI 95%')
plt.axvline(x = bootstrap_ci_upper, color = "pink", linestyle = "dotted")
plt.legend()
plt.show()<jupyter_output><empty_output><jupyter_text>Do men pay more for the ticket than women? $$ 
H_0: \mu_{male} - \mu_{female} \le 0
$$ $$
H_1: \mu_{male} - \mu_{female} \gt 0 
$$  Traditional Approach

The traditional approach to this type of hypothesis is to perform a t-test for the means of two independent samples, a.k.a. Two Sample T-Test.<jupyter_code># define the significance level for our hypotesis tests
alpha = 0.05
# separate the gender's fare's in two dataframes
male_fares_df = sample_df[sample_df['sex']=='male']['fare']
female_fares_df = sample_df[sample_df['sex']=='female']['fare']
# Performe the levene test to check if the variances are statistically equal or not,
_, pvalue_levene = stats.levene(male_fares_df, female_fares_df)

# use levene's result to decide if the variances are equal or not 
equal_variances = pvalue_levene < alpha

equal_variances
# Use scypy.stats.ttest_ind to perform the two sample t-test
_, pvalue_ttest = stats.ttest_ind(male_fares_df, female_fares_df, equal_var=equal_variances)
_, pvalue_ttest
# Since the ttest_ind H1 is for the differences on both ends of a t-distribution,
# and we're only interested in one end, we're going to divide the p-value by 2
pvalue_ttest /= 2
pvalue_ttest
# Use the t-test two sample result to evaluate H0
if (pvalue_ttest < alpha):
  print("We have evidence to reject the null hypothesis. Men pay more for tickets than women.")
else:
  print("We don't have evidence to reject the null hypothesis. Men don't pay more for tickets than women.")<jupyter_output>We don't have evidence to reject the null hypothesis. Men don't pay more for tickets than women.<jupyter_text>Bootstrap Approach

1. Take the sample
2. Take n bootstrap samples from the sample
3. Compute the mean for each bootstrap sample
4. Simulate the distribution under H0
5. Compare the booststrap sample distribution and the H0 distribution<jupyter_code>bootstrap_means_male, bootstrap_means_female = [], []

# take n bootstrap samples from the sample
for _ in range(10000):
  # take one bootstrap sample
  bootstrap_sample = sample_df.sample(sample_size, replace=True)
  # compute the statistic for each sample and add to the list
  bootstrap_means_male.append(bootstrap_sample[bootstrap_sample['sex']=='male']['fare'].mean())
  bootstrap_means_female.append(bootstrap_sample[bootstrap_sample['sex']=='female']['fare'].mean())

# use the lists to generate a DataFrame (for ease of use)
bootstrap_means_df = pd.DataFrame({'male':bootstrap_means_male, 'female':bootstrap_means_female})

# compute the mean diff between the two 'lists'
bootstrap_means_df['diff'] = bootstrap_means_df['male']-bootstrap_means_df['female']
bootstrap_means_df.head()
print(bootstrap_means_df['male'].mean())
print(bootstrap_means_df['female'].mean())
# Plot both (male and female) bootstraps distributions of the means 
bootstrap_means_df[['male', 'female']].plot(kind='hist', bins=20, alpha=.2)
# Plot the diff distribution 
bootstrap_means_df['diff'].plot(kind='hist', bins=20, alpha=.2)
# Simulate the distribution under the null hypothesis
diff_means_h0 = []
for _ in range(10000):
  diff_means_h0.append(np.random.normal(0, bootstrap_means_df['diff'].std(), size=sample_size).mean())

bootstrap_means_df['diff_h0'] = diff_means_h0
# Plot the null distribution 
bootstrap_means_df['diff_h0'].plot(kind='hist', bins=20, alpha=.2)
# Plot both (diff and h0) bootstraps distributions of the means 
bootstrap_means_df[['diff', 'diff_h0']].plot(kind='hist', bins=200, alpha=.2)
# Out of N Bootstrap executions, how many times we observed the boot sample mean
# Greater than the Null Hipótesis
(bootstrap_means_df['diff'] >= bootstrap_means_df['diff_h0']).value_counts()
# Compute the bootstrap equivalent of the t-statistic
(bootstrap_means_df['diff'].mean()-0)/bootstrap_means_df['diff'].std()
pvalue_bootstrap = (bootstrap_means_df['diff'] > bootstrap_means_df['diff_h0']).mean()
pvalue_bootstrap
# Use the pvalue to evaluate H0
if (pvalue_bootstrap < alpha):
  print("We have evidence to reject the null hypothesis. Men pay more for tickets than women.")
else:
  print("We don't have evidence to reject the null hypothesis. Men pay as much as women for tickets.")<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.10366706097429831, 'fraction_numerical': 0.015377997634154211, 'mean_word_length': 5.671449067431851, 'pattern_counts': {'":': 1, '<': 31, '<?xml version=': 0, '>': 30, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '7', 'dir': 'jupyter-structured-clean-dedup', 'id': '520734', 'license': 'MIT', 'n_tokens_mistral': 3164, 'n_tokens_neox': 2897, 'n_words': 941, 'path': 'inferential_statistics_traditional_x_bootstraping.ipynb', 'repo_name': 'thiagofuruchima/neptune', 'size': '9299'}