<jupyter_start><jupyter_text>Table of Contents1&nbsp;&nbsp;Optimal Parameters1.1&nbsp;&nbsp;Poisson and exponential distribution - MLB no hitters2&nbsp;&nbsp;Linear regression by least squares2.1&nbsp;&nbsp;Investigating Anscombe data3&nbsp;&nbsp;Bootstrapping3.1&nbsp;&nbsp;Visualizing bootstrap samples4&nbsp;&nbsp;Pairs bootstrap4.1&nbsp;&nbsp;Bootstrapping the Linear regression parameters5&nbsp;&nbsp;Hypothesis testing5.1&nbsp;&nbsp;Permutation samples5.2&nbsp;&nbsp;Plotting the ECDF of permutation samples5.3&nbsp;&nbsp;Test statistic and p-values5.4&nbsp;&nbsp;Sample Bootstrap hypothesis testing6&nbsp;&nbsp;Bootstrap test for identical distributions6.1&nbsp;&nbsp;A two-sample bootstrap hypothesis test for difference of means.6.2&nbsp;&nbsp;A/B testing - The vote for the Civil Rights Act in 19647&nbsp;&nbsp;Hypothesis test on correlation (permutation sampling)8&nbsp;&nbsp;Do neonicotinoid have unintended consequences ?8.1&nbsp;&nbsp;EDA and comparing the empirical distribution functions8.2&nbsp;&nbsp;Bootstrap hypothesis test on bee sperm counts9&nbsp;&nbsp;Investigating the Darwin Finches dataset9.1&nbsp;&nbsp;EDA of beak depths of Darwin's finches9.2&nbsp;&nbsp;Parameter estimates of beak depths9.3&nbsp;&nbsp;Hypothesis test: Are beaks deeper in 2012?9.4&nbsp;&nbsp;Comparing evolution of beak length and beak depth9.5&nbsp;&nbsp;Linear regression9.6&nbsp;&nbsp;Beak length to depth ratio9.7&nbsp;&nbsp;EDA of heritability  Optimal Parameters<jupyter_code>import os
import pandas as pd
import numpy as np
import scipy.stats as stats
import matplotlib.pyplot as plt
import seaborn as sns
%matplotlib inline
sns.set() #Making seaborn the default styling 
data_path = '/Users/User/Desktop/Data/DataCamp-master/07-statistical-thinking-in-python-(part-2)/_datasets'
os.chdir(data_path)
os.listdir()
def ECDF(data, lab = 'ECDF'):
    #plt.figure(figsize = (8,4))
    x = np.sort(data)
    n = len(x)
    y = np.arange(1, n+1) / n
    percentiles = np.percentile(x, [25,50,75])
    plt.plot(x,y, linestyle = 'none', marker = '.' , label = lab) 
    plt.plot(percentiles, (.25,.50,.75), marker = 'D', color = 'red', linestyle = 'none')
 
    plt.xlabel('x_axis')
    plt.ylabel('Cumulative probability')
    plt.title('ECDF')
    plt.legend()
    plt.margins(0.02)
    
def ECDF_vals(data):
    x = np.sort(data)
    n = len(data)
    y = np.arange(1, n+1) / n
    return x,y<jupyter_output><empty_output><jupyter_text>Poisson and exponential distribution - MLB no hitters- If an event is described by the poisson distribution, then the time between events is described by the exponential distribution- Exponential distribution has a single parameter: tau which is the average time between eventsWe will check here whether the data from the MLB no hitters does indeed follow such distributions using hacker statistics<jupyter_code>data = pd.read_csv('mlb_nohitters.csv')
data.head()
#Extracting the number of games between no hitters
games = data['game_number'].values
nohitter_times = np.empty(294)
for i in range(293,0,-1):
    nohitter_times[i] = (games[i] - games[i-1]) -1
nohitter_times = nohitter_times[1:,] 
#Calculating mean from data and simulating the exponential distribution
mean_time = np.mean(nohitter_times)
sim_nohitter_times = np.random.exponential(mean_time,100000)
sim_nohitter_times_half = np.random.exponential(mean_time/2,100000)
sim_nohitter_times_double = np.random.exponential(mean_time*2,100000)

#Plotting the simulated and actual data CDF
plt.figure(figsize = (8,6))
ECDF(sim_nohitter_times, 'tau = mean')
ECDF(nohitter_times, 'Data')
ECDF(sim_nohitter_times_half, 'tau/2')
ECDF(sim_nohitter_times_double, '2*tau')
plt.legend()<jupyter_output><empty_output><jupyter_text>Linear regression by least squares<jupyter_code>data = pd.read_csv('female_literacy_fertility.csv')
illiteracy = 100 - data['female literacy'].values
fertility = data['fertility'].values
data.head()
plt.plot(illiteracy, fertility, marker='.', linestyle='none')<jupyter_output><empty_output><jupyter_text>Linear regression using numpy polyfit<jupyter_code># Plot the illiteracy rate versus fertility
_ = plt.plot(illiteracy, fertility, marker='.', linestyle='none')
plt.margins(0.02)
_ = plt.xlabel('percent illiterate')
_ = plt.ylabel('fertility')

# Perform a linear regression using np.polyfit(): a, b
a, b = np.polyfit(x = illiteracy, y = fertility, deg = 1)

# Print the results to the screen
print('slope =', a, 'children per woman / percent illiterate')
print('intercept =', b, 'children per woman')

# Make theoretical line to plot
x = np.array([0,100])
y = a * x + b

# Add regression line to your plot
_ = plt.plot(x, y)

# Draw the plot
plt.show()<jupyter_output>slope = 0.0497985480906342 children per woman / percent illiterate
intercept = 1.8880506106365562 children per woman<jupyter_text>Showing the polyfit parameters are optimal<jupyter_code>a_vals = np.linspace(0,0.1,200)
RSS = np.empty_like(a_vals)

for i,a in enumerate(a_vals):
    RSS[i] = np.sum((fertility - a*illiteracy - b)**2)

plt.plot(a_vals, RSS)
plt.plot(a_vals[np.argmin(RSS)], np.min(RSS), marker = 'o', color = 'Red')<jupyter_output><empty_output><jupyter_text>Investigating Anscombe data<jupyter_code>anscombe = sns.load_dataset("anscombe")
sets = ['I','II','III','IV']
df_dict = {sets[i]:anscombe[anscombe['dataset'] == v] for i,v in enumerate(sets)}
fig, axes = plt.subplots(1,4,figsize= (22,6), sharex = True, sharey = True)

for i,v in enumerate(df_dict):
    data = df_dict[v][['x','y']]
    a,b = np.polyfit(data['x'],data['y'],deg = 1)
    x_theor = np.array([3, 15])
    y_theor = a * x_theor + b
    axes[i].scatter(data['x'],data['y'])
    axes[i].plot(x_theor,y_theor)

plt.show()
def ECDF(data, lab = 'ECDF'):
    x = np.sort(data)
    n = len(x)
    y = np.arange(1, n+1) / n
    plt.plot(x,y, linestyle = 'none', marker = '.', color = 'gray', alpha = 0.1)  
    plt.xlabel('x_axis'); plt.ylabel('Cumulative probability'); plt.title('ECDF'); plt.legend(); plt.margins(0.02)<jupyter_output><empty_output><jupyter_text>Bootstrapping Definitions:- Bootstrap sample: - Bootstrap replicate: *Pre-processing the dataset*<jupyter_code>weather = pd.read_csv('sheffield_weather_station.csv', skiprows= 8, delim_whitespace= True, parse_dates = [0])
weather.set_index('yyyy', inplace = True, drop = False)
rainfall = weather['rain']
rainfall = rainfall.resample('AS').sum()
#rainfall = year_rain[0:len(year_rain) -1]<jupyter_output><empty_output><jupyter_text>Visualizing bootstrap samples<jupyter_code>plt.figure(figsize = (8,8))

for _ in range(100):  #Acquiring 100 bootstrap samples
    bs_sample = np.random.choice(rainfall, size = len(rainfall)) #Using np.choice function to select with replacement
    ECDF(bs_sample)

x = np.sort(rainfall)
n = len(rainfall)
y = np.arange(1, n+1) / n
plt.plot(x,y, marker = '.') 
plt.show()<jupyter_output>No handles with labels found to put in legend.
No handles with labels found to put in legend.
No handles with labels found to put in legend.
No handles with labels found to put in legend.
No handles with labels found to put in legend.
No handles with labels found to put in legend.
No handles with labels found to put in legend.
No handles with labels found to put in legend.
No handles with labels found to put in legend.
No handles with labels found to put in legend.
No handles with labels found to put in legend.
No handles with labels found to put in legend.
No handles with labels found to put in legend.
No handles with labels found to put in legend.
No handles with labels found to put in legend.
No handles with labels found to put in legend.
No handles with labels found to put in legend.
No handles with labels found to put in legend.
No handles with labels found to put in legend.
No handles with labels found to put in legend.
No handles with labels found to put in legend.
No handles wi[...]<jupyter_text>Bootstrap replicates<jupyter_code>def bootstrap_replicate_1d(data,func):
    '''Generate bootstrap replicates of 1d daya. 
    Note that np.choice performs random sampling WITH replacement'''
    #Generate a boostrap sample by choosing data with replacement randomly
    bs_sample = np.random.choice(data, len(data))
    #Generate the boostrap replicate by applying a function to compute a statistic for the bootstrap sample
    return func(bs_sample)

def draw_bs_reps(data, func, size=1):
    '''Draw bootstrap replicates'''
    bs_replicates = np.empty(size)
    
    for i in range(size):
        bs_replicates[i] = bootstrap_replicate_1d(data,func)

    return bs_replicates<jupyter_output><empty_output><jupyter_text>Bootstrap replicate of the meanCompute a bootstrap estimate of the probability distribution function of the mean annual rainfall at the Sheffield Weather Station. This is a probabilistic estimate of the mean. Plot the PDF as a histogram, and see that it is Normal.<jupyter_code>bs_mean = draw_bs_reps(rainfall,np.mean,10000)
plt.hist(bs_mean, bins = 100, normed = True);
print('Mean of bootstrap mean = ', np.mean(bs_mean))
print('Var of bootstrap mean = ', np.var(bs_mean))
print('Theoretical mean of mean = ', np.std(rainfall) / np.sqrt(len(rainfall)))

#Plotting confidence intervals using percentiles
percentiles = np.percentile(bs_mean, [2.5,97.5])
plt.axvline(x= percentiles[0], color = 'red');
plt.axvline(x= percentiles[1], color = 'red');
plt.title('PDF of the bootstrap mean with 90% confidence intervals')
plt.show()<jupyter_output>Mean of bootstrap mean =  798.090022164179
Var of bootstrap mean =  113.15178106202823
Theoretical mean of mean =  10.635458130769607<jupyter_text>Bootstrap replicates of the variance<jupyter_code>bs_var = draw_bs_reps(rainfall,np.var,10000)
plt.hist(bs_var, bins = 100, normed = True);
print('Mean of bootstrap variance = ', np.mean(bs_var))
print('Var of bootstrap variance = ', np.var(bs_var))
print('Theoretical mean of variance = ', (len(rainfall) - 1)/len(rainfall) * np.var(rainfall) ) 

#Plotting confidence intervals using percentiles
percentiles_var = np.percentile(bs_var, [2.5,97.5])
plt.axvline(x= percentiles_var[0], color = 'red');
plt.axvline(x= percentiles_var[1], color = 'red');
plt.title('PDF of the bootstrap mean with 90% confidence intervals')
plt.show()<jupyter_output>Mean of bootstrap variance =  15036.225406357986
Var of bootstrap variance =  2927884.9032053766
Theoretical mean of variance =  15044.024963629994<jupyter_text>Confidence interval on the rate of no-hittersConsider again the inter-no-hitter intervals for the modern era of baseball. Generate 10,000 bootstrap replicates of the optimal parameter τ. Plot a histogram of your replicates and report a 95% confidence interval.<jupyter_code>tau = np.mean(nohitter_times)
bs_tau = draw_bs_reps(nohitter_times, np.mean,10000)
plt.hist(bs_tau, bins = 100, normed = True);
print('Mean of bootstrap tau = ', np.mean(bs_tau))
print('Var of bootstrap tau = ', np.var(bs_tau))
print('Theoretical mean of variance = ', tau ) 

#Plotting confidence intervals using percentiles
percentiles_tau = np.percentile(bs_tau, [2.5,97.5])
plt.axvline(x= percentiles_tau[0], color = 'red');
plt.axvline(x= percentiles_tau[1], color = 'red');
plt.title('PDF of the bootstrap mean with 95% confidence intervals')
plt.show()
percentiles_tau<jupyter_output><empty_output><jupyter_text>Pairs bootstrap Bootstrapping the Linear regression parameters<jupyter_code>def draw_bs_pairs_linreg(x, y, size=1):
    """Perform pairs bootstrap for linear regression."""

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

    # Initialize replicates: bs_slope_reps, bs_intercept_reps
    bs_slope_reps = np.empty(size)
    bs_intercept_reps = np.empty(size)

    # Generate replicates
    for i in range(size):
        bs_inds = np.random.choice(inds, size = len(inds))
        bs_x, bs_y = x[bs_inds], y[bs_inds]
        bs_slope_reps[i], bs_intercept_reps[i] = np.polyfit(bs_x, bs_y, 1)

    return bs_slope_reps, bs_intercept_reps

# Generate replicates of slope and intercept using pairs bootstrap
bs_slope_reps, bs_intercept_reps = draw_bs_pairs_linreg(illiteracy, fertility, 10000)

# Compute and print 95% CI for slope
print(np.percentile(bs_slope_reps,[2.5, 97.5]))

# Plot the histogram
_ = plt.hist(bs_slope_reps, bins=50, normed=True)
_ = plt.xlabel('slope')
_ = plt.ylabel('PDF')
plt.show()
# Generate array of x-values for bootstrap lines: x
x = np.array([0, 100])

# Plot the bootstrap lines
for i in range(1000):
    _ = plt.plot(x, bs_slope_reps[i] * x + bs_intercept_reps[i],
                 linewidth=0.5, alpha=0.01, color='red')

# Plot the data
_ = plt.plot(illiteracy, fertility, marker='.', linestyle='none')

# Label axes, set the margins, and show the plot
_ = plt.xlabel('illiteracy')
_ = plt.ylabel('fertility')
plt.margins(0.02)
plt.show()<jupyter_output><empty_output><jupyter_text>Hypothesis testing  Permutation samples<jupyter_code>def permutation_sample(data1,data2):
    '''Generate a permutation sample from two data sets'''
    data = np.concatenate((data1,data2))
    permutated_data = np.random.permutation(data)
    perm_sample1 = permutated_data[:len(data1)]
    perm_sample2 = permutated_data[len(data1):]
    return perm_sample1, perm_sample2   
rain_july = np.array([ 66.2,  39.7,  76.4,  26.5,  11.2,  61.8,   6.1,  48.4,  89.2,       104. ,  34. ,  60.6,  57.1,  79.1,  90.9,  32.3,  63.8,  78.2,     27.5,  43.4,  30.1,  17.3,  77.5,  44.9,  92.2,  39.6,  79.4, 66.1,  53.5,  98.5,  20.8,  55.5,  39.6,  56. ,  65.1,  14.8,
13.2,  88.1,   8.4,  32.1,  19.6,  40.4,   2.2,  77.5, 105.4, 77.2,  38. ,  27.1, 111.8,  17.2,  26.7,  23.3,  77.2,  87.2,        27.7,  50.6,  60.3,  15.1,   6. ,  29.4,  39.3,  56.3,  80.4, 85.3,  68.4,  72.5,  13.3,  28.4,  14.7,  37.4,  49.5,  57.2,        85.9,  82.1,  31.8, 126.6,  30.7,  41.4,  33.9,  13.5,  99.1, 70.2,  91.8,  61.3,  13.7,  54.9,  62.5,  24.2,  69.4,  83.1, 44. ,  48.5,  11.9,  16.6,  66.4,  90. ,  34.9, 132.8,  33.4, 225. ,   7.6,  40.9,  76.5,  48. , 140. ,  55.9,  54.1,  46.4,        68.6,  52.2, 108.3,  14.6,  11.3,  29.8, 130.9, 152.4,  61. , 46.6,  43.9,  30.9, 111.1,  68.5,  42.2,   9.8, 285.6,  56.7,        168.2,  41.2,  47.8, 166.6,  37.8,  45.4,  43.2])

rain_november = np.array([ 83.6,  30.9,  62.2,  37. ,  41. , 160.2,  18.2, 122.4,  71.3,         44.2,  49.1,  37.6, 114.5,  28.8,  82.5,  71.9,  50.7,  67.7, 112. ,  63.6,  42.8,  57.2,  99.1,  86.4,  84.4,  38.1,  17.7, 102.2, 101.3,  58. ,  82. , 101.4,  81.4, 100.1,  54.6,  39.6, 57.5,  29.2,  48.8,  37.3, 115.4,  55.6,  62. ,  95. ,  84.2, 118.1, 153.2,  83.4, 104.7,  59. ,  46.4,  50. , 147.6,  76.8,
        59.9, 101.8, 136.6, 173. ,  92.5,  37. ,  59.8, 142.1,   9.9, 158.2,  72.6,  28. , 112.9, 119.3, 199.2,  50.7,  44. , 170.7, 67.2,  21.4,  61.3,  15.6, 106. , 116.2,  42.3,  38.5, 132.5, 40.8, 147.5,  93.9,  71.4,  87.3, 163.7, 141.4,  62.6,  84.9, 28.8, 121.1,  28.6,  32.4, 112. ,  50. ,  96.9,  81.8,  70.4, 117.5,  41.2, 124.9,  78.2,  93. ,  53.5,  50.5,  42.6,  47.9, 73.1, 129.1,  56.9, 103.3,  60.5, 134.3,  93.1,  49.5,  48.2, 167.9,  27. , 111.1,  55.4,  36.2,  57.4,  66.8,  58.3,  60. , 161.6, 112.7,  37.4, 110.6,  56.6,  95.8, 126.8])<jupyter_output><empty_output><jupyter_text>Plotting the ECDF of permutation samples<jupyter_code>plt.figure(figsize = (8,8));

for i in range(100):
    perm_sample1, perm_sample2 = permutation_sample(rain_july, rain_november)
    x_1, y_1 = ECDF_vals(perm_sample1)
    x_2, y_2 = ECDF_vals(perm_sample2)
    
    plt.plot(x_1,y_1, marker = '.', linestyle = 'None', color = 'red', alpha = 0.02);
    plt.plot(x_2,y_2, marker = '.', linestyle = 'None', color = 'blue', alpha = 0.02);
    

x_1, y_1 = ECDF_vals(rain_july)
x_2, y_2 = ECDF_vals(rain_november)
    
plt.plot(x_1,y_1, marker = '.', linestyle = 'None', color = 'red');
plt.plot(x_2,y_2, marker = '.', linestyle = 'None', color = 'blue');
plt.margins(0.02); plt.xlabel('monthly rainfall (mm)'); plt.ylabel('ECDF'); plt.show()<jupyter_output><empty_output><jupyter_text>Test statistic and p-values Writing a function to draw permutation replicates<jupyter_code>def draw_perm_reps(data_1, data_2, func, size=1):
    '''Generate multiple permutation replicates. 
    Here func is a function that takes two arrays as arguments'''
    
    perm_replicates = np.empty(size)
    
    for i in range(size):
        perm_sample_1, perm_sample_2 = permutation_sample(data_1, data_2)
        perm_replicates[i] = func(perm_sample_1,perm_sample_2)

    return perm_replicates    
frog_tongue = pd.read_csv('frog_tongue.csv', skiprows=14)
df = frog_tongue[frog_tongue['ID'].isin(['II', 'IV'])][['ID','impact force (mN)']].copy()
df['impact force (N)'] = df['impact force (mN)'] / 1000 
df.drop(axis = 1, labels='impact force (mN)', inplace = True)
force_a = df[df['ID'] == 'II']['impact force (N)'].copy().values
force_b = df[df['ID'] == 'IV']['impact force (N)'].copy().values
sns.swarmplot(x = 'ID', y = 'impact force (N)', data=df)
display(df.groupby('ID').describe())
print('The difference between the means is :' , round(0.70735 - 0.4191,2))<jupyter_output><empty_output><jupyter_text>Hypothesis test of difference of meansCompute the probability of getting at least a 0.29 N difference in mean strike force under the hypothesis that the distributions of strike forces for the two frogs are identical. We use a permutation test with a test statistic of the difference of means to test this hypothesis. - H0 : The two samples have the same mean- H1: Mean force a is greater than mean force b<jupyter_code>def diff_of_means(data1,data2):
    '''Compute the difference of means'''
    return np.mean(data1) - np.mean(data2)
means_rep = draw_perm_reps(force_a, force_b, diff_of_means, 100000)
plt.hist(means_rep, bins = 100, normed = True);

#Plotting 95 % confidence intervals using percentiles
percentiles_diff = np.percentile(means_rep, [2.5,97.5])
plt.axvline(x= percentiles_diff[0], color = 'red');
plt.axvline(x= percentiles_diff[1], color = 'red');
plt.title('PDF of the difference of means with 95% confidence intervals')
plt.show()

print('Empirical difference of means = ', diff_of_means(force_a, force_b)) 
print('95% confidence intervals ', percentiles_diff, 'With range = ', abs(percentiles_diff[0] - percentiles_diff[1]) ) 
#Compute the p_value
empirical_diff_means = diff_of_means(force_a, force_b)

p_value = sum(means_rep >= empirical_diff_means ) / len(means_rep)
p_value * 100<jupyter_output><empty_output><jupyter_text>Given that the p value is below 0.06, we can reject the null hypothesis with a confidence level of 94%  Sample Bootstrap hypothesis testing Another juvenile frog was studied, Frog C, and you want to see if Frog B and Frog C have similar impact forces. Unfortunately, you do not have Frog C's impact forces available, but you know they have a mean of 0.55 N. Because you don't have the original data, you cannot do a permutation test, and you cannot assess the hypothesis that the forces from Frog B and Frog C come from the same distribution. You will therefore test another, less restrictive hypothesis: The mean strike force of Frog B is equal to that of Frog C.To set up the bootstrap hypothesis test, you will take the mean as our test statistic. Remember, your goal is to calculate the probability of getting a mean impact force less than or equal to what was observed for Frog B if the hypothesis that the true mean of Frog B's impact forces is equal to that of Frog C is true. You first translate all of the data of Frog B such that the mean is 0.55 N. This involves adding the mean force of Frog C and subtracting the mean force of Frog B from each measurement of Frog B. This leaves other properties of Frog B's distribution, such as the variance, unchanged. - H0: mean frog B = 0.55N- H1: mean frog B < 0.55NObserved mean of frog b = 0.4191000000000001<jupyter_code>#Under H0 the mean is .55 so we translate force_b to match H0 
translated_force_b = force_b - np.mean(force_b) + 0.55

#Perform bootstrapping and replicate calculation
bs_replicates = draw_bs_reps(translated_force_b, np.mean, 50000)

# Describe statistics
stats.describe(bs_replicates)
fig, (ax1, ax2) = plt.subplots(1,2, figsize = (16,6))
#Distribution plot
ax1 = sns.distplot(bs_replicates, ax = ax1)
ax1.axvline(0.4191, ymin = 0, ymax = 0.0455, color='r')

#QQ plot 
ax2 = stats.probplot(bs_replicates, dist = "norm", plot = plt)
#Resulting distribution is clearly normal
#Calculating the p value of the observed mean: 0.4191000000000001
P = (bs_replicates <= np.mean(force_b)).sum()/bs_replicates.size
print("P value is: ", P)<jupyter_output>P value is:  0.00592<jupyter_text>Bootstrap test for identical distributions Perform a bootstrap test of the hypothesis that Frog A and Frog B have identical distributions of impact forces using the difference of means test statistic. Note that one-sample bootstrap hypothesis test are impossible to do with permutation.**Guidance: **To do this test on two arrays with n1 and n2 entries, we do a very similar procedure as a permutation test. We concatenate the arrays, generate a bootstrap sample from it, and take the first n1 entries of the bootstrap sample as belonging to the first data set and the last n2 as belonging to the second. We then compute the test statistic, e.g., the difference of means, to get a bootstrap replicate. The p-value is the number of bootstrap replicates for which the test statistic is less than what was observed.<jupyter_code># Compute difference of mean impact force from experiment: empirical_diff_means
empirical_diff_means = np.mean(force_a) - np.mean(force_b)

# Concatenate forces: forces_concat
forces_concat = np.concatenate([force_a , force_b])

# Initialize bootstrap replicates: bs_replicates
bs_replicates = np.empty(10000)

for i in range(10000):
    # Generate bootstrap sample
    bs_sample = np.random.choice(forces_concat, size=len(forces_concat))
    
    # Compute replicate
    bs_replicates[i] = diff_of_means(bs_sample[:len(force_a)], bs_sample[len(force_b):])

p = (bs_replicates>=empirical_diff_means).sum()/len(bs_replicates)
fig, (ax1, ax2) = plt.subplots(1,2, figsize = (16,6))
#Distribution plot
ax1 = sns.distplot(bs_replicates, ax = ax1)
ax1.axvline(empirical_diff_means, ymin = 0, ymax = 0.0455, color='r')

#QQ plot 
ax2 = stats.probplot(bs_replicates, dist = "norm", plot = plt)
#Resulting distribution is clearly normal<jupyter_output><empty_output><jupyter_text>A two-sample bootstrap hypothesis test for difference of means. You performed a one-sample bootstrap hypothesis test, which is impossible to do with permutation. Testing the hypothesis that two samples have the same distribution may be done with a bootstrap test, but a permutation test is preferred because it is more accurate (exact, in fact). But therein lies the limit of a permutation test; it is not very versatile. We now want to test the hypothesis that Frog A and Frog B have the same mean impact force, but not necessarily the same distribution. This, too, is impossible with a permutation test.To do the two-sample bootstrap test, we shift both arrays to have the same mean, since we are simulating the hypothesis that their means are, in fact, equal. We then draw bootstrap samples out of the shifted arrays and compute the difference in means. This constitutes a bootstrap replicate, and we generate many of them. The p-value is the fraction of replicates with a difference in means greater than or equal to what was observed. - H0: mean frog A = mean frog B - H1: mean frog A $\neq$ mean frog B<jupyter_code>#Mean of all forces
mean_force = np.mean(forces_concat)

#Shift both samples under H0 
force_a_shifted = force_a - np.mean(force_a) + mean_force
force_b_shifted = force_b - np.mean(force_b) + mean_force

#Compute 10,000 bootstrap replicates from shifted arrays
bs_replicates_a = draw_bs_reps(force_a_shifted, np.mean, 10000)
bs_replicates_b = draw_bs_reps(force_b_shifted, np.mean, 10000)

# Get replicates of difference of means: bs_replicates
bs_replicates = (bs_replicates_a - bs_replicates_b)

# Compute and print p-value: p
p = np.sum(bs_replicates >= empirical_diff_means) / len(bs_replicates)
print('p-value =', p)                           
fig, (ax1, ax2) = plt.subplots(1,2, figsize = (16,6))
#Distribution plot
ax1 = sns.distplot(bs_replicates, ax = ax1)
ax1.axvline(empirical_diff_means, ymin = 0, ymax = 0.0455, color='r')

#QQ plot 
ax2 = stats.probplot(bs_replicates, dist = "norm", plot = plt)
#Resulting distribution is clearly normal<jupyter_output><empty_output><jupyter_text>A/B testing - The vote for the Civil Rights Act in 1964 The Civil Rights Act of 1964 was one of the most important pieces of legislation ever passed in the USA. Excluding "present" and "abstain" votes, 153 House Democrats and 136 Republicans voted yay. However, 91 Democrats and 35 Republicans voted nay. Did party affiliation make a difference in the vote?<jupyter_code># Construct arrays of data: dems, reps
dems = np.array([True] * 153 + [False] * 91)
reps = np.array([True] * 136 + [False] * 35)

def frac_yay_dems(dems, reps):
    """Compute fraction of Democrat yay votes."""
    frac = np.sum(dems) / len(dems)
    return frac

# Acquire permutation samples: perm_replicates
perm_replicates = draw_perm_reps(dems, reps, frac_yay_dems, 10000)

# Compute and print p-value: p
p = np.sum(perm_replicates <= 153/244) / len(perm_replicates)
print('p-value =', p)
fig, (ax1, ax2) = plt.subplots(1,2, figsize = (16,6))
#Distribution plot
ax1 = sns.distplot(perm_replicates, ax = ax1)

#QQ plot 
ax2 = stats.probplot(perm_replicates, dist = "norm", plot = plt)
#Resulting distribution is clearly normal<jupyter_output><empty_output><jupyter_text>Hypothesis test on correlation (permutation sampling)Testing the correlation between illiteracy and fertility rate in the set of 162 countries. Permute the illiteracy values but leave the fertility values fixed. This simulates the hypothesis that they are totally independent of each other. For each permutation, compute the Pearson correlation coefficient and assess how many of your permutation replicates have a Pearson correlation coefficient greater than the observed one.<jupyter_code>np.corrcoef(illiteracy,fertility)[0,1]
# Compute observed correlation: r_obs
r_obs = np.corrcoef(illiteracy,fertility)[0,1]

# Initialize permutation replicates: perm_replicates
perm_replicates = np.empty(10000)

# Draw replicates
for i in range(10000):
    # Permute illiteracy measurments: illiteracy_permuted
    illiteracy_permuted = np.random.permutation(illiteracy)

    # Compute Pearson correlation
    perm_replicates[i] = np.corrcoef(illiteracy_permuted,fertility)[0,1]

# Compute p-value: p
p = np.sum(perm_replicates >= r_obs) / len(perm_replicates)
print('p-val =', p)
plt.plot(illiteracy, fertility, marker='.', linestyle='none')<jupyter_output><empty_output><jupyter_text>Do neonicotinoid have unintended consequences ?<jupyter_code>bee_sperm = pd.read_csv("bee_sperm.csv", skiprows = 3)
display(bee_sperm.tail())
control = bee_sperm[bee_sperm["Treatment"] == "Control"]["Alive Sperm Millions"].values
treated = bee_sperm[bee_sperm["Treatment"] == "Pesticide"]["Alive Sperm Millions"].values

def ECDF(data):
    x = np.sort(data)
    n = len(data)
    y = np.arange(1, n+1) / n
    return x,y<jupyter_output><empty_output><jupyter_text>EDA and comparing the empirical distribution functions<jupyter_code># Compute x,y values for ECDFs
x_control, y_control = ECDF(control)
x_treated, y_treated = ECDF(treated)

# Plot the ECDFs
plt.plot(x_control, y_control, marker='.', linestyle='none')
plt.plot(x_treated, y_treated, marker='.', linestyle='none')

# Set the margins
plt.margins(0.02)

# Add a legend
plt.legend(('control', 'treated'), loc='lower right')

# Label axes and show plot
plt.xlabel('millions of alive sperm per mL')
plt.ylabel('ECDF')
plt.show()<jupyter_output><empty_output><jupyter_text>Bootstrap hypothesis test on bee sperm countsTest the following hypothesis: On average, male bees treated with neonicotinoid insecticide have the same number of active sperm per milliliter of semen than do untreated male bees.<jupyter_code># Compute the difference in mean sperm count: diff_means
diff_means = diff_of_means(control,treated)

# Compute mean of pooled data: mean_count
mean_count = np.mean(np.concatenate((control,treated)))

# Generate shifted data sets
control_shifted = control - np.mean(control) + mean_count
treated_shifted = treated - np.mean(treated) + mean_count

# Generate bootstrap replicates
bs_reps_control = draw_bs_reps(control_shifted, np.mean, size=10000)
bs_reps_treated = draw_bs_reps(treated_shifted, np.mean, size=10000)

# Get replicates of difference of means: bs_replicates
bs_replicates = bs_reps_control - bs_reps_treated

# Compute and print p-value: p
p = np.sum(bs_replicates >= np.mean(control) - np.mean(treated)) \
            / len(bs_replicates)
print('p-value =', p)

fig, (ax1, ax2) = plt.subplots(1,2, figsize = (16,6))
#Distribution plot
ax1 = sns.distplot(bs_replicates, ax = ax1)

#QQ plot 
ax2 = stats.probplot(bs_replicates, dist = "norm", plot = plt)
#Resulting distribution is clearly normal<jupyter_output><empty_output><jupyter_text>Investigating the Darwin Finches dataset EDA of beak depths of Darwin's finchesFor your first foray into the Darwin finch data, you will study how the beak depth (the distance, top to bottom, of a closed beak) of the finch species Geospiza scandens has changed over time. The Grants have noticed some changes of beak geometry depending on the types of seeds available on the island, and they also noticed that there was some interbreeding with another major species on Daphne Major, Geospiza fortis. These effects can lead to changes in the species over time.In the next few problems, you will look at the beak depth of G. scandens on Daphne Major in 1975 and in 2012. To start with, let's plot all of the beak depth measurements in 1975 and 2012 in a bee swarm plot.<jupyter_code>finch_beaks_1975 = pd.read_csv("finch_beaks_1975.csv")
finch_beaks_2012 = pd.read_csv("finch_beaks_2012.csv")
finch_beaks_1975.rename(index=str, columns={"Beak length, mm": "blength","Beak depth, mm":"bdepth" }, inplace = True)
df = pd.concat([finch_beaks_1975,finch_beaks_2012], keys = ['1975','2012'], names = ["Year"]).reset_index(level = 0, )
df.head()
scandens_bdepth_2012 = finch_beaks_2012[finch_beaks_2012["species"]=="scandens"]["bdepth"].values
scandens_bdepth_1975 = finch_beaks_1975[finch_beaks_1975["species"]=="scandens"]["bdepth"].values
scandens_blength_2012 = finch_beaks_2012[finch_beaks_2012["species"]=="scandens"]["blength"].values
scandens_blength_1975 = finch_beaks_1975[finch_beaks_1975["species"]=="scandens"]["blength"].values
# Create bee swarm plot
_ = sns.swarmplot(x = 'Year', y = 'bdepth', data=df[df["species"]=="scandens"])

# Label the axes
_ = plt.xlabel('year')
_ = plt.ylabel('beak depth (mm)')

# Show the plot
plt.show()<jupyter_output><empty_output><jupyter_text>ECDF of beak depths<jupyter_code># Compute ECDFs
x_1975, y_1975 = ECDF(scandens_bdepth_1975)
x_2012, y_2012 = ECDF(scandens_bdepth_2012)

# Plot the ECDFs
_ = plt.plot(x_1975, y_1975, marker='.', linestyle='none')
_ = plt.plot(x_2012, y_2012, marker='.', linestyle='none')

# Set margins
plt.margins(0.02)

# Add axis labels and legend
_ = plt.xlabel('beak depth (mm)')
_ = plt.ylabel('ECDF')
_ = plt.legend(('1975', '2012'), loc='lower right')

# Show the plot
plt.show()<jupyter_output><empty_output><jupyter_text>Parameter estimates of beak depths Estimate the *difference* of the mean beak depth of the G. scandens samples from 1975 and 2012 and report a 95% confidence interval.<jupyter_code># Compute the difference of the sample means: mean_diff
mean_diff = diff_of_means(scandens_bdepth_2012, scandens_bdepth_1975)

# Get bootstrap replicates of means
bs_replicates_1975 = draw_bs_reps(scandens_bdepth_1975, np.mean, 10000)
bs_replicates_2012 = draw_bs_reps(scandens_bdepth_2012, np.mean, 10000)

# Compute samples of difference of means: bs_diff_replicates
bs_diff_replicates = bs_replicates_2012 - bs_replicates_1975

# Compute 95% confidence interval: conf_int
conf_int = np.percentile(bs_diff_replicates, [2.5,97.5])

# Print the results
print('difference of means =', mean_diff, 'mm')
print('95% confidence interval =', conf_int, 'mm')<jupyter_output>difference of means = 0.22622047244094645 mm
95% confidence interval = [0.05943936 0.39024783] mm<jupyter_text>Hypothesis test: Are beaks deeper in 2012?Your plot of the ECDF and determination of the confidence interval make it pretty clear that the beaks of G. scandens on Daphne Major have gotten deeper. But is it possible that this effect is just due to random chance? In other words, what is the probability that we would get the observed difference in mean beak depth if the means were the same?Be careful! The hypothesis we are testing is not that the beak depths come from the same distribution. For that we could use a permutation test. The hypothesis is that the means are equal. To perform this hypothesis test, we need to shift the two data sets so that they have the same mean and then use bootstrap sampling to compute the difference of means.<jupyter_code># Compute mean of combined data set: combined_mean
combined_mean = np.mean(np.concatenate((scandens_bdepth_1975, scandens_bdepth_2012)))

# Shift the samples
bd_1975_shifted = scandens_bdepth_1975 - np.mean(scandens_bdepth_1975) + combined_mean
bd_2012_shifted = scandens_bdepth_2012 - np.mean(scandens_bdepth_2012) + combined_mean

# Get bootstrap replicates of shifted data sets
bs_replicates_1975 = draw_bs_reps(bd_1975_shifted, np.mean, 10000)
bs_replicates_2012 = draw_bs_reps(bd_2012_shifted, np.mean, 10000)

# Compute replicates of difference of means: bs_diff_replicates
bs_diff_replicates = bs_replicates_2012 - bs_replicates_1975

# Compute the p-value
p = np.sum(bs_diff_replicates >= mean_diff) / len(bs_diff_replicates)

# Print p-value
print('p =', p)<jupyter_output>p = 0.003<jupyter_text>Comparing evolution of beak length and beak depth EDA of beak length and depth<jupyter_code># Make scatter plot of 1975 data
_ = plt.plot(scandens_blength_1975, scandens_bdepth_1975, marker='.',
             linestyle='none', alpha = 0.5, color = 'blue')

# Make scatter plot of 2012 data
_ = plt.plot(scandens_blength_2012, scandens_bdepth_2012, marker='.',
             linestyle='none', alpha = 0.5, color = 'red')

# Label axes and make legend
_ = plt.xlabel('beak length (mm)')
_ = plt.ylabel('beak depth (mm)')
_ = plt.legend(('1975', '2012'), loc='upper left')

# Show the plot
plt.show()<jupyter_output><empty_output><jupyter_text>Linear regression<jupyter_code># Compute the linear regressions
slope_1975, intercept_1975 = np.polyfit(x = scandens_blength_1975, y = scandens_bdepth_1975, deg = 1)
slope_2012, intercept_2012 = np.polyfit(x = scandens_blength_2012, y = scandens_bdepth_2012, deg = 1)

# Perform pairs bootstrap for the linear regressions
bs_slope_reps_1975, bs_intercept_reps_1975 = draw_bs_pairs_linreg(scandens_blength_1975,scandens_bdepth_1975, 1000 )
bs_slope_reps_2012, bs_intercept_reps_2012 = draw_bs_pairs_linreg(scandens_blength_2012,scandens_bdepth_2012, 1000 )

# Compute confidence intervals of slopes
slope_conf_int_1975 = np.percentile(bs_slope_reps_1975, [2.5,97.5])
slope_conf_int_2012 = np.percentile(bs_slope_reps_2012, [2.5,97.5])
intercept_conf_int_1975 = np.percentile(bs_intercept_reps_1975, [2.5,97.5])

intercept_conf_int_2012 = np.percentile(bs_intercept_reps_2012, [2.5,97.5])


# Print the results
print('1975: slope =', slope_1975,
      'conf int =', slope_conf_int_1975)
print('1975: intercept =', intercept_1975,
      'conf int =', intercept_conf_int_1975)
print('2012: slope =', slope_2012,
      'conf int =', slope_conf_int_2012)
print('2012: intercept =', intercept_2012,
      'conf int =', intercept_conf_int_2012)
# Make scatter plot of 1975 data
_ = plt.plot(scandens_blength_1975, scandens_bdepth_1975, marker='.',
             linestyle='none', color='blue', alpha=0.5)

# Make scatter plot of 2012 data
_ = plt.plot(scandens_blength_2012, scandens_bdepth_2012, marker='.',
             linestyle='none', color='red', alpha=0.5)

# Label axes and make legend
_ = plt.xlabel('beak length (mm)')
_ = plt.ylabel('beak depth (mm)')
_ = plt.legend(('1975', '2012'), loc='upper left')

# Generate x-values for bootstrap lines: x
x = np.array([10, 17])

# Plot the bootstrap lines
for i in range(100):
    plt.plot(x, + x * bs_slope_reps_1975[i] +  bs_intercept_reps_1975[i],
             linewidth=0.5, alpha=0.2, color='blue')
    plt.plot(x, + x * bs_slope_reps_2012[i] +  bs_intercept_reps_2012[i],
             linewidth=0.5, alpha=0.2, color='red')

# Draw the plot again
plt.show()<jupyter_output><empty_output><jupyter_text>Beak length to depth ratio<jupyter_code># Compute length-to-depth ratios
ratio_1975 = scandens_blength_1975 / scandens_bdepth_1975
ratio_2012 = scandens_blength_2012 / scandens_bdepth_2012

# Compute means
mean_ratio_1975 = np.mean(ratio_1975)
mean_ratio_2012 = np.mean(ratio_2012)

# Generate bootstrap replicates of the means
bs_replicates_1975 = draw_bs_reps(ratio_1975, np.mean,10000)
bs_replicates_2012 = draw_bs_reps(ratio_2012, np.mean,10000)

# Compute the 99% confidence intervals
conf_int_1975 = np.percentile(bs_replicates_1975,[.5,99.5])
conf_int_2012 = np.percentile(bs_replicates_2012,[.5,99.5])

# Print the results
print('1975: mean ratio =', mean_ratio_1975,
      'conf int =', conf_int_1975)
print('2012: mean ratio =', mean_ratio_2012,
      'conf int =', conf_int_2012)<jupyter_output>1975: mean ratio = 1.5788823771858533 conf int = [1.55695077 1.6015921 ]
2012: mean ratio = 1.4658342276847767 conf int = [1.44451022 1.48850976]
{'fraction_non_alphanumeric': 0.1154998276913289, 'fraction_numerical': 0.0626408292023434, 'mean_word_length': 4.937047529115517, 'pattern_counts': {'":': 2, '<': 124, '<?xml version=': 0, '>': 127, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 2, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '32', 'dir': 'jupyter-structured-clean-dedup', 'id': '364185', 'license': 'MIT', 'n_tokens_mistral': 15103, 'n_tokens_neox': 13174, 'n_words': 3855, 'path': 'Stats_Thinking_2_Bootstrap_permutation_AB.ipynb', 'repo_name': 'draperkm/Data-Portfolio', 'size': '37723'}