<jupyter_start><jupyter_text>In this notebook you'll create your own bootstrap function following the bootstrap algorithm (check the lecture notes!)<jupyter_code>import matplotlib
import pandas as pd
import seaborn as sns
import numpy as np
# Load the data
df = pd.read_csv('https://raw.githubusercontent.com/albanda/CE888/master/lab2%20-%20bootstrap/customers.csv')
data = df.values.T[1]
# Checking the notes from the lecture, create here your own bootstrap function:
# 1. Sample from the input array x to create an array of samples of shape (n_bootstraps, sample_size)
# Hint: Check the function random.choice() on Numpy
# 2. Calculate and save the mean of the array (this is "data_mean" that is returned by the function)
# 3. Calculate the mean from each bootstrap (i.e., row) and store it.
# (This should be an array of n_bootstraps values)
# 4. Calculate the lower and upper bounds for a 95% CI (hint: check the percentile function on Numpy)
# 5. Return data_mean, and the lower and upper bounds of your interval
def bootstrap_mean(x, sample_size, n_bootstraps):
	lists = []
	meandata = []
	for i in range(n_bootstraps):
		lists.append(np.random.choice(x, sample_size, True))
	for sample_set in lists:
		meandata.append(np.mean(sample_set))
	meandata.sort()
	data_mean = np.mean(meandata)
	lower = np.percentile(meandata, 2.5)
	upper = np.percentile(meandata, 97.5)
 # <---INSERT YOUR CODE HERE--->
	return data_mean, lower, upper

# Call your bootstrap function and plot the results
boots = []
for i in range(100, 50000, 1000):
    boot = bootstrap_mean(data, data.shape[0], i)
    boots.append([i, boot[0], "mean"])
    boots.append([i, boot[1], "lower"])
    boots.append([i, boot[2], "upper"])

df_boot = pd.DataFrame(boots, columns=['Bootstrap Iterations', 'Mean', "Value"])
sns_plot = sns.lmplot(df_boot.columns[0], df_boot.columns[1], data=df_boot, fit_reg=False, hue="Value")

sns_plot.axes[0, 0].set_ylim(0,15)
sns_plot.axes[0, 0].set_xlim(0, 100000)<jupyter_output>/usr/local/lib/python3.7/dist-packages/seaborn/_decorators.py:43: FutureWarning: Pass the following variables as keyword args: x, y. From version 0.12, the only valid positional argument will be `data`, and passing other arguments without an explicit keyword will result in an error or misinterpretation.
  FutureWarning<jupyter_text>Now, modify the bootstrap function you created above so that you can pass your desired confidence interval as a parameter.<jupyter_code>def bootstrap_stdev_ci(sample, sample_size, n_bootstraps, ci):
  lists = []
  std_data1 = []
  for i in range(n_bootstraps):
    lists.append(np.random.choice(sample, sample_size, True))
  for sample_set in lists:
    std_data1.append(np.std(sample_set))
  std_data1.sort()
  data_std = np.std(std_data1)
  lower = np.percentile(std_data1, (100-ci)/2)
  upper = np.percentile(std_data1, 100-((100-ci)/2))
  return data_std, lower, upper


def bootstrap_mean_ci(sample, sample_size, n_bootstraps, ci):
	lists = []
	mean_data1 = []
	for i in range(n_bootstraps):
		lists.append(np.random.choice(sample, sample_size, True))
	for sample in lists:
		mean_data1.append(np.std(sample))
	mean_data1.sort()
	data_mean = np.std(mean_data1)
	lower = np.percentile(mean_data1, (100-ci)/2)
	upper = np.percentile(mean_data1, 100-((100-ci)/2))
	# <---INSERT YOUR CODE HERE--->
	return data_mean, lower, upper
boots = []
for i in range(100, 50000, 1000):
    boot = bootstrap_mean_ci(data, data.shape[0], i, 80)
    boots.append([i, boot[0], "mean"])
    boots.append([i, boot[1], "lower"])
    boots.append([i, boot[2], "upper"])

df_boot = pd.DataFrame(boots, columns=['Boostrap Iterations', 'Mean', "Value"])
sns_plot = sns.lmplot(df_boot.columns[0], df_boot.columns[1], data=df_boot, fit_reg=False, hue="Value")

sns_plot.axes[0, 0].set_ylim(0, 15)
sns_plot.axes[0, 0].set_xlim(0, 100000)

#sns_plot.savefig("bootstrap_confidence_80.pdf", bbox_inches='tight')

data_mean, lower, upper = bootstrap_mean_ci(data, data.shape[0], 10000, 90)
print('lower =',lower)
print('upper =',upper)

data_mean, lower, upper = bootstrap_mean_ci(data, data.shape[0], 10000, 80)
print('lower =',lower)
print('upper =',upper)

data_std, lower, upper = bootstrap_stdev_ci(data, data.shape[0], 10000, 85)
print('lower =',lower)
print('upper =',upper)

data_std, lower, upper = bootstrap_stdev_ci(data, data.shape[0], 10000, 99)
print('lower =',lower)
print('upper =',upper)<jupyter_output>lower = 3.043770772547368
upper = 7.149996432138399
lower = 3.6146668285883847
upper = 6.782518020778474
lower = 3.320468371191274
upper = 6.956497476665764
lower = 2.1957693536676266
upper = 7.910237802712791<jupyter_text>Vehicles datasetNow let's work on a different dataset, which is stored in the vehicles.csv file.<jupyter_code># Load and visualise the vehicles dataset
# To load the dataset: https://neptune.ai/blog/google-colab-dealing-with-files (check section "Load individual files directly from GitHub")


# Note that the current and new fleets are in different columns and have different lengths, so bear this in mind when you're plotting.
# You can create separate scatterplots for the two fleets, as you would with the histograms, 
# or plot them both in one plot (but not one against the other).
# <---INSERT YOUR CODE HERE--->
df = pd.read_csv('https://raw.githubusercontent.com/albanda/CE888/master/lab2%20-%20bootstrap/vehicles.csv')

current = df.values.T[0][np.logical_not(np.isnan(df.values.T[0]))]
new = df.values.T[1][np.logical_not(np.isnan(df.values.T[1]))]
# Note: you can add more cells as needed to organise your code and your plots<jupyter_output><empty_output><jupyter_text>Compare the two fleetsThe business analysts come up a comparison algorithm that requires the upper and lower bounds for the mean in order to say which fleet is better.1. Calculate the mean of both samples.2. Using the bootstrap function that you created:    - Construct the 95% CI of the mean of the current fleet.    - Construct the 95% CI of the mean of the new fleet.    - Are they comparable? (i.e., is one better than the other?) -- you can do this with a permutation test (check the lecture notes!)<jupyter_code># <---INSERT YOUR CODE HERE--->
boot_current = []
for i in range(100, 50000, 1000):
    cboot = bootstrap_mean_ci(current, data.shape[0], i, 95)
    boot_current.append([i, cboot[0], "mean"])
    boot_current.append([i, cboot[1], "lower"])
    boot_current.append([i, cboot[2], "upper"])
df_boot_c = pd.DataFrame(boot_current, columns=['Boostrap Iterations', 'Mean', "Value"])
sns_plot_c = sns.lmplot(df_boot_c.columns[0], df_boot_c.columns[1], data=df_boot_c, fit_reg=False, hue="Value")

sns_plot_c.axes[0, 0].set_ylim(0,25)
sns_plot_c.axes[0, 0].set_xlim(0, 100000)

boot_new = []
for j in range(100, 50000, 1000):
    nboot = bootstrap_mean_ci(new, data.shape[0], j, 95)
    boot_new.append([j, nboot[0], "mean"])
    boot_new.append([j, nboot[1], "lower"])
    boot_new.append([j, nboot[2], "upper"])

df_boot_n = pd.DataFrame(boot_new, columns=['Boostrap Iterations', 'Mean', "Value"])
sns_plot_n = sns.lmplot(df_boot_n.columns[0], df_boot_n.columns[1], data=df_boot_n, fit_reg=False, hue="Value")

sns_plot_n.axes[0, 0].set_ylim(0,)
sns_plot_n.axes[0, 0].set_xlim(0, 100000)
# Create your own function for a permutation test here (you will need it for the lab quiz!):
def permut_test(sample1, sample2, n_permutations):
    """
    sample1: 1D array
    sample2: 1D array (note that the size of the two arrays can be different)
    n_permutations: number of permutations to calculate the p-value
    """
    concat = np.concatenate((sample1, sample2))
    mean_o, mean_n = np.mean(sample1),np.mean(sample2)
    t_obs = mean_n - mean_o
    ctr = 0
    for i in range(n_permutations):
      perm = np.random.permutation(concat)
      po = perm[:len(sample1)]
      pn = perm[len(sample1):]
      mean_po, mean_pn = np.mean(po), np.mean(pn)
      t_perm = mean_pn - mean_po
      if t_perm > t_obs:
        ctr += 1
    pvalue = ctr / n_permutations
    return pvalue
print(permut_test(current,new,30000))

print(bootstrap_stdev_ci(new, len(new), 10000, 73))
print(bootstrap_mean_ci(new, new.shape[0], 10000, 80))
print(bootstrap_mean_ci(data, data.shape[0], 10000, 95))
print(bootstrap_stdev_ci(current, len(current), 10000, 95))
dem_share_PA = [60.08, 40.64, 36.07, 41.21, 31.04, 43.78, 44.08, 46.85, 44.71, 46.15, 63.10, 52.20, 43.18, 40.24, 39.92, 47.87, 37.77, 40.11, 49.85, 48.61, 38.62, 54.25, 34.84, 47.75, 43.82, 55.97, 58.23, 42.97, 42.38, 36.11, 37.53, 42.65, 50.96, 47.43, 56.24, 45.60, 46.39, 35.22, 48.56, 32.97, 57.88, 36.05, 37.72, 50.36, 32.12, 41.55, 54.66, 57.81, 54.58, 32.88, 54.37, 40.45, 47.61, 60.49, 43.11, 27.32, 44.03, 33.56, 37.26, 54.64, 43.12, 25.34, 49.79, 83.56, 40.09, 60.81, 49.81]
dem_share_OH = [56.94, 50.46, 65.99, 45.88, 42.23, 45.26, 57.01, 53.61, 59.10, 61.48, 43.43, 44.69, 54.59, 48.36, 45.89, 48.62, 43.92, 38.23, 28.79, 63.57, 38.07, 40.18, 43.05, 41.56, 42.49, 36.06, 52.76, 46.07, 39.43, 39.26, 47.47, 27.92, 38.01, 45.45, 29.07, 28.94, 51.28, 50.10, 39.84, 36.43, 35.71, 31.47, 47.01, 40.10, 48.76, 31.56, 39.86, 45.31, 35.47, 51.38, 46.33, 48.73, 41.77, 41.32, 48.46, 53.14, 34.01, 54.74, 40.67, 38.96, 46.29, 38.25, 6.80, 31.75, 46.33, 44.90, 33.57, 38.10, 39.67, 40.47, 49.44, 37.62, 36.71, 46.73, 42.20, 53.16, 52.40, 58.36, 68.02, 38.53, 34.58, 69.64, 60.50, 53.53, 36.54, 49.58, 41.97, 38.11]

print(len(dem_share_OH))
print(len(dem_share_PA))
 
print(bootstrap_mean_ci(dem_share_PA, len(dem_share_PA), 25000, 90))
print(bootstrap_mean_ci(dem_share_OH, len(dem_share_OH), 25000, 90))

print('p_value = ',permut_test(dem_share_OH,dem_share_PA,20000))<jupyter_output>88
67
(1.1671801288192998, 7.825365858522503, 11.627339895370902)
(0.9865224884255013, 8.157692737251708, 11.397315797477411)
p_value =  0.23705
StarcoderdataJupyterStructuredCleanDedup