<jupyter_start><jupyter_text>Bike commute times<jupyter_code>data = pd.read_csv("data-src/bike-commute-times.csv")
#data
fig = plt.figure(figsize=(6,4))
ax1 = fig.add_axes([0.15, 0.15, 0.8, 0.8])
    
for ax in [ax1]:
    bp=ax.boxplot([data[data.Frame=="steel"].Time, data[data.Frame=="carbon"].Time], labels=["Steel\nframe", "Carbon\nframe"],
              boxprops = dict(linewidth=2, color=colors["lightgray"]),
              flierprops = dict(marker='o', markeredgecolor=colors["lightgray"], linestyle='none'),
              whiskerprops = dict(color=colors["lightgray"], linewidth=2),
              capprops = dict(color=colors["lightgray"], linewidth=2),
              medianprops = dict(color=colors["orange"], linewidth=2))
    for spine in ["bottom", "left"]:
        ax.spines[spine].set_linewidth(1)
        ax.spines[spine].set_color(colors["lightgray"])
    for spine in ["top", "right"]:
        ax.spines[spine].set_visible(False)
    ax.set_ylabel("Commuting time", size=18, color=colors["lightgray"], ha="center")
    ax.tick_params(axis="y", width=1, size=4, color=colors["lightgray"], labelcolor=colors["lightgray"], 
                   labelsize=13, pad=4)
    ax.tick_params(axis="x", width=1, size=0, labelcolor=colors["lightgray"], 
                   labelsize=16, pad=6)
    ## change color of the medians
    for color,median in zip([colors["orange"], colors["blue"]], bp['medians']):
        median.set(color=color, linewidth=2)
    
plt.savefig("/Users/Gui/Box Sync/_STATS13/_Slides/_stats13-Lectures/assets/img/lec/bike-commute-boxplot.svg", transparent=True)
print("Carbon frame:\n", "mean = {:.3f}, std = {:.3f}".format(data[data.Frame=="carbon"].Time.mean(), data[data.Frame=="carbon"].Time.std()))
print("Steel frame:\n", "mean = {:.3f}, std = {:.3f}".format(data[data.Frame=="steel"].Time.mean(), data[data.Frame=="steel"].Time.std()))<jupyter_output>Carbon frame:
 mean = 108.344, std = 6.248
Steel frame:
 mean = 107.809, std = 4.892<jupyter_text>Simulation null hypothesis<jupyter_code>n_simul = 10000
res_diff_mean = np.zeros(n_simul) #store the data
res_diff_median = np.zeros(n_simul) #store the data

n_bike_carbon = len(data[data.Frame=="carbon"])
n_bike_steel = len(data[data.Frame=="steel"])

all_times = data.Time.values.copy()

#each simulation is the shuffling of the full population and the
#calculation of the difference in proportion of shifts with death
for i in range(n_simul):
    np.random.shuffle(all_times)
    carbon_times = all_times[:n_bike_carbon]
    steel_times = all_times[n_bike_carbon:]
    res_diff_mean[i] = np.mean(carbon_times)-np.mean(steel_times)
    res_diff_median[i] = np.median(carbon_times)-np.median(steel_times)<jupyter_output><empty_output><jupyter_text>Simulation steps<jupyter_code>x1, y1 = st.to_dotplot(res_diff_mean[:120], kind="bins", scale=0.05, nbins=20)

#randomize order of each rows of data points
rows_yval = np.unique(y1)
idx_by_rows = [list(np.where(y1 == val)[0]) for val in rows_yval]
for i in range(len(idx_by_rows)):
    np.random.shuffle(idx_by_rows[i])
shuffled_idx = np.concatenate(idx_by_rows)

x1 = x1[shuffled_idx]
y1 = y1[shuffled_idx]

fig, ax1 = plt.subplots(figsize=(6, 4))
ax2 = ax1.twinx()
ax3 = ax1.twinx()
ax4 = ax1.twinx()
ax5 = ax1.twinx()
ax6 = ax1.twinx()
ax7 = ax1.twinx()
ax8 = ax1.twinx()

for ax in [ax1]:
    for spine in ["bottom"]:
        ax.spines[spine].set_linewidth(1)
        ax.spines[spine].set_color(colors["lightgray"])
    for spine in ["top", "right", "left"]:
        ax.spines[spine].set_visible(False)
    ax.set_yticks([])
    ax.tick_params(axis="both", width=1, size=4, color=colors["lightgray"], labelcolor=colors["lightgray"], 
                   labelsize=13, pad=4)
    ax.set_xlabel(r'$\bar{x}_{\mathrm{carbon}}-\bar{x}_{\mathrm{steel}}$ (min)', size=18, color=colors["lightgray"])
    ax.set_ylim(-0.02, 1.2)
    ax.set_xlim(res_diff_mean.min()-1, res_diff_mean.max()+1)
    ax.ticklabel_format(style='sci', axis='x', scilimits=(0,5))

for ax in [ax2, ax3, ax4, ax5, ax6, ax7, ax8]:
    ax.set_ylim(ax1.get_ylim())
    ax.axis("off")
    
for ax in [ax2]:
    ax.scatter(x1[0], y1[0], s=50, color=colors["blue"])

for ax in [ax3]:
    ax.scatter(x1[1], y1[1], s=50, color=colors["blue"])

for ax in [ax4]:
    ax.scatter(x1[2], y1[2], s=50, color=colors["blue"])

for ax in [ax5]:
    ax.scatter(x1[3:25], y1[3:25], color=colors["blue"])

for ax in [ax6]:
    ax.scatter(x1[25:50], y1[25:50], color=colors["blue"])

for ax in [ax7]:
    ax.scatter(x1[50:80], y1[50:80], color=colors["blue"])

for ax in [ax8]:
    ax.scatter(x1[80:], y1[80:], color=colors["blue"])

plt.tight_layout()

plt.savefig("/Users/Gui/Box Sync/_STATS13/_Slides/_stats13-Lectures/assets/img/lec/bike-commute-simulation-diff-mean.svg", transparent=True)
x1, y1 = st.to_dotplot(res_diff_median[:120], kind="bins", scale=0.05, nbins=20)

#randomize order of each rows of data points
rows_yval = np.unique(y1)
idx_by_rows = [list(np.where(y1 == val)[0]) for val in rows_yval]
for i in range(len(idx_by_rows)):
    np.random.shuffle(idx_by_rows[i])
shuffled_idx = np.concatenate(idx_by_rows)

x1 = x1[shuffled_idx]
y1 = y1[shuffled_idx]

fig, ax1 = plt.subplots(figsize=(6, 4))
ax2 = ax1.twinx()
ax3 = ax1.twinx()
ax4 = ax1.twinx()
ax5 = ax1.twinx()
ax6 = ax1.twinx()
ax7 = ax1.twinx()
ax8 = ax1.twinx()

for ax in [ax1]:
    for spine in ["bottom"]:
        ax.spines[spine].set_linewidth(1)
        ax.spines[spine].set_color(colors["lightgray"])
    for spine in ["top", "right", "left"]:
        ax.spines[spine].set_visible(False)
    ax.set_yticks([])
    ax.tick_params(axis="both", width=1, size=4, color=colors["lightgray"], labelcolor=colors["lightgray"], 
                   labelsize=13, pad=4)
    ax.set_xlabel('$\mathrm{m}_{\mathrm{carbon}}-\mathrm{m}_{\mathrm{steel}}$ (min)', size=18, color=colors["lightgray"])
    ax.set_ylim(-0.02, 1.2)
    ax.set_xlim(res_diff_median.min()-1, res_diff_median.max()+1)
    ax.ticklabel_format(style='sci', axis='x', scilimits=(0,5))

for ax in [ax2, ax3, ax4, ax5, ax6, ax7, ax8]:
    ax.set_ylim(ax1.get_ylim())
    ax.axis("off")
    
for ax in [ax2]:
    ax.scatter(x1[0], y1[0], s=50, color=colors["blue"])

for ax in [ax3]:
    ax.scatter(x1[1], y1[1], s=50, color=colors["blue"])

for ax in [ax4]:
    ax.scatter(x1[2], y1[2], s=50, color=colors["blue"])

for ax in [ax5]:
    ax.scatter(x1[3:25], y1[3:25], color=colors["blue"])

for ax in [ax6]:
    ax.scatter(x1[25:50], y1[25:50], color=colors["blue"])

for ax in [ax7]:
    ax.scatter(x1[50:80], y1[50:80], color=colors["blue"])

for ax in [ax8]:
    ax.scatter(x1[80:], y1[80:], color=colors["blue"])

plt.tight_layout()

plt.savefig("/Users/Gui/Box Sync/_STATS13/_Slides/_stats13-Lectures/assets/img/lec/bike-commute-simulation-diff-median.svg", transparent=True)<jupyter_output><empty_output><jupyter_text>Full simulation (10000)<jupyter_code>#Observed diff
bike_diff_mean = data[data.Frame=="carbon"].Time.mean()-data[data.Frame=="steel"].Time.mean()
bike_diff_mean_se = np.sqrt((data[data.Frame=="carbon"].Time.std()**2/n_bike_carbon)+(data[data.Frame=="steel"].Time.std()**2/n_bike_steel))

#t distribution
tdist_bike_diff = stats.t(df=len(data)-2, loc=0, scale=bike_diff_mean_se)
xt_bike_diff = np.linspace(-6.5, 6.5, 1000)
yt_bike_diff = tdist_bike_diff.pdf(xt_bike_diff)

fig = plt.figure(figsize=(6,4))
ax1 = fig.add_axes([0.1, 0.15, 0.82, 0.75])
ax2 = ax1.twinx()
ax3 = ax1.twinx()
ax4 = ax1.twinx()
ax5 = ax1.twinx()
    
for ax in [ax1]:
    ax.hist(res_diff_mean, bins="auto", color=colors["blue"])
    for spine in ["bottom"]:
        ax.spines[spine].set_linewidth(1)
        ax.spines[spine].set_color(colors["lightgray"])
    for spine in ["top", "right", "left"]:
        ax.spines[spine].set_visible(False)
    ax.set_yticks([])
    ax.tick_params(axis="both", width=1, size=4, color=colors["lightgray"], labelcolor=colors["lightgray"], 
                   labelsize=13, pad=4)
    ax.set_xlabel(r'$\bar{x}_{\mathrm{carbon}}-\bar{x}_{\mathrm{steel}}$ (min)', size=18, color=colors["lightgray"])
    ax.set_xlim(xt_bike_diff.min(), xt_bike_diff.max())
    ax.set_ylim(0)
    
for ax in [ax2]:
    ax.set_ylim(ax1.get_ylim())
    ax.axvline(bike_diff_mean, color=colors["red"], ls="--")
    ax.text(bike_diff_mean, ax.get_ylim()[1]*1.01, "Initial observation\n" r"$\bar{{x}}_{{\mathrm{{carbon}}}}-\bar{{x}}_{{\mathrm{{steel}}}}={:.2f}$".format(bike_diff_mean), color=colors["red"], size=13, ha="left")
    ax.axis("off")
    
for ax in [ax3]:
    ax.set_ylim(ax1.get_ylim())
    ax.axvline(0-bike_diff_mean, color=colors["red"], ls="--")
    ax.text(bike_diff_mean+3.5, ax.get_ylim()[1]*0.5, "{:.0f} simulations\n$\geq${:.2f}".format(np.sum(res_diff_mean>=bike_diff_mean), bike_diff_mean), color=colors["lightgray"], size=13, ha="center")
    ax.text(bike_diff_mean-4.5, ax.get_ylim()[1]*0.5, "{:.0f} simulations\n$\leq${:.2f}".format(np.sum(res_diff_mean<=0-bike_diff_mean), 0-bike_diff_mean), color=colors["lightgray"], size=13, ha="center")
    n_extreme = np.sum(res_diff_mean<=0-bike_diff_mean)+np.sum(res_diff_mean>=bike_diff_mean)
    ax.text(ax.get_xlim()[0], ax.get_ylim()[1]*1, "2-tail p-value\np = {:.3f}".format(n_extreme/10000), color=colors["lightgray"], size=13, ha="left")
    ax.axis("off")

for ax in [ax4]:
    lower_lim = 0-bike_diff_mean
    upper_lim = bike_diff_mean
    ax.fill_between(xt_bike_diff, yt_bike_diff,
                       color=colors["orange"], alpha=0.7)
    ax.set_ylim(0)
    ax.axis("off")

for ax in [ax5]:
    lower_lim = 0-bike_diff_mean
    upper_lim = bike_diff_mean
    ax.fill_between(xt_bike_diff[(xt_bike_diff>lower_lim) & (xt_bike_diff<upper_lim)], yt_bike_diff[(xt_bike_diff>lower_lim) & (xt_bike_diff<upper_lim)],
                       color=colors["orange"], alpha=0.7)
    ax.fill_between(xt_bike_diff[(xt_bike_diff<=lower_lim)], yt_bike_diff[(xt_bike_diff<=lower_lim)],
                       color=colors["red"], alpha=0.7)
    ax.fill_between(xt_bike_diff[(xt_bike_diff>=upper_lim)], yt_bike_diff[(xt_bike_diff>=upper_lim)],
                       color=colors["red"], alpha=0.7)
    ax.set_ylim(0)
    norm_pvalue_2tail = tdist_bike_diff.sf(0.53/bike_diff_mean_se)*2
    ax.text(ax.get_xlim()[0], ax.get_ylim()[1]*0.76, "2-tail p-value\n(normal approximation)\np = {:.3f}".format(norm_pvalue_2tail), color=colors["lightgray"], size=13, ha="left")
    ax.text(0.65, ax.get_ylim()[1]*0.91, "$\pi_0+t*\sigma$", color=colors["lightgray"], size=13, ha="left")
    ax.axis("off")


plt.savefig("/Users/Gui/Box Sync/_STATS13/_Slides/_stats13-Lectures/assets/img/lec/bike-commute-simulation-diff-mean-all.svg", transparent=True)
t_stat = (108.34-107.81)/np.sqrt(6.25**2/26+4.89**2/30)
t_stat
#Observed diff
bike_diff_median = data[data.Frame=="carbon"].Time.median()-data[data.Frame=="steel"].Time.median()

fig = plt.figure(figsize=(6,4))
ax1 = fig.add_axes([0.1, 0.15, 0.82, 0.75])
ax2 = ax1.twinx()
ax3 = ax1.twinx()
    
for ax in [ax1]:
    ax.hist(res_diff_median, bins="auto", color=colors["blue"])
    for spine in ["bottom"]:
        ax.spines[spine].set_linewidth(1)
        ax.spines[spine].set_color(colors["lightgray"])
    for spine in ["top", "right", "left"]:
        ax.spines[spine].set_visible(False)
    ax.set_yticks([])
    ax.tick_params(axis="both", width=1, size=4, color=colors["lightgray"], labelcolor=colors["lightgray"], 
                   labelsize=13, pad=4)
    ax.set_xlabel('$\mathrm{m}_{\mathrm{carbon}}-\mathrm{m}_{\mathrm{steel}}$ (min)', size=18, color=colors["lightgray"])
    ax.set_xlim(-7, 7)
    ax.set_ylim(0)
    
for ax in [ax2]:
    ax.set_ylim(ax1.get_ylim())
    ax.axvline(bike_diff_median, color=colors["orange"], ls="--")
    ax.text(bike_diff_median, ax.get_ylim()[1]*1.01, "Initial observation\n" "$\mathrm{{m}}_{{\mathrm{{carbon}}}}-\mathrm{{m}}_{{\mathrm{{steel}}}}={:.2f}$".format(bike_diff_median), color=colors["orange"], size=13, ha="left")
    ax.axis("off")
    
for ax in [ax3]:
    ax.set_ylim(ax1.get_ylim())
    ax.axvline(0-bike_diff_median, color=colors["orange"], ls="--")
    ax.text(bike_diff_median+3.5, ax.get_ylim()[1]*0.5, "{:.0f} simulations\n$\geq${:.2f}".format(np.sum(res_diff_median>=bike_diff_median), bike_diff_median), color=colors["red"], size=13, ha="center")
    ax.text(bike_diff_median-4.5, ax.get_ylim()[1]*0.5, "{:.0f} simulations\n$\leq${:.2f}".format(np.sum(res_diff_median<=0-bike_diff_median), 0-bike_diff_median), color=colors["red"], size=13, ha="center")
    n_extreme = np.sum(res_diff_median<=0-bike_diff_median)+np.sum(res_diff_median>=bike_diff_median)
    ax.text(ax.get_xlim()[0], ax.get_ylim()[1]*1, "2-tail p-value\np = {:.3f}".format(n_extreme/10000), color=colors["orange"], size=13, ha="left")
    ax.axis("off")


plt.savefig("/Users/Gui/Box Sync/_STATS13/_Slides/_stats13-Lectures/assets/img/lec/bike-commute-simulation-diff-median-all.svg", transparent=True)<jupyter_output><empty_output><jupyter_text>Simulation 95% confidence intervals<jupyter_code>n_simul = 10000
res_median = np.zeros(n_simul) #store the data

n_bike_carbon = len(data[data.Frame=="carbon"])
n_bike_steel = len(data[data.Frame=="steel"])

#each simulation is the difference in median time from isolated bootstrap samples
for i in range(n_simul):
    bootstrap_carbon = np.random.choice(data[data.Frame=="carbon"].Time, size=n_bike_carbon)
    bootstrap_steel = np.random.choice(data[data.Frame=="steel"].Time, size=n_bike_steel)
    res_median[i] = np.median(bootstrap_carbon)-np.median(bootstrap_steel)
#Observed diff
bike_diff_median = data[data.Frame=="carbon"].Time.median()-data[data.Frame=="steel"].Time.median()

fig = plt.figure(figsize=(6,4))
ax1 = fig.add_axes([0.1, 0.15, 0.82, 0.75])
ax2 = ax1.twinx()
ax3 = ax1.twinx()
    
for ax in [ax1]:
    ax.hist(res_median, bins="auto", color=colors["blue"])
    for spine in ["bottom"]:
        ax.spines[spine].set_linewidth(1)
        ax.spines[spine].set_color(colors["lightgray"])
    for spine in ["top", "right", "left"]:
        ax.spines[spine].set_visible(False)
    ax.set_yticks([])
    ax.tick_params(axis="both", width=1, size=4, color=colors["lightgray"], labelcolor=colors["lightgray"], 
                   labelsize=13, pad=4)
    ax.set_xlabel('$\mathrm{m}_{\mathrm{carbon}}-\mathrm{m}_{\mathrm{steel}}$ (min)', size=18, color=colors["lightgray"])
    ax.set_xlim(-7, 7)
    ax.set_ylim(0)
    
for ax in [ax2]:
    ax.set_ylim(ax1.get_ylim())
    ax.axvline(bike_diff_median, color=colors["orange"], ls="--")
    ax.text(bike_diff_median, ax.get_ylim()[1]*1.01, "Initial observation\n" "$\mathrm{{m}}_{{\mathrm{{carbon}}}}-\mathrm{{m}}_{{\mathrm{{steel}}}}={:.2f}$".format(bike_diff_median), color=colors["orange"], size=13, ha="center")
    ax.axis("off")

for ax in [ax3]:
    ax.axvline(np.percentile(res_median, 2.5), ymax=0.55, color=colors["red"], lw=2)
    ax.axvline(np.percentile(res_median, 97.5), ymax=0.55, color=colors["red"], lw=2)
    ax.text(np.percentile(res_median, 2.5), ax.get_ylim()[1]*0.6, "2.5$^{{th}}$\npercentile\n{:.2f}".format(np.percentile(res_median, 2.5)), color=colors["red"], size=15, ha="center")
    ax.text(np.percentile(res_median, 97.5), ax.get_ylim()[1]*0.6, "97.5$^{{th}}$\npercentile\n{:.2f}".format(np.percentile(res_median, 97.5)), color=colors["red"], size=15, ha="center")
    ax.axis("off")


plt.savefig("/Users/Gui/Box Sync/_STATS13/_Slides/_stats13-Lectures/assets/img/lec/bike-commute-simulation-diff-median-resampling-ci95.svg", transparent=True)<jupyter_output><empty_output><jupyter_text>t test<jupyter_code>fig = plt.figure(figsize=(12,4))
ax1 = fig.add_axes([0.065, 0.155, 0.2, 0.65])
ax2 = fig.add_axes([0.315, 0.155, 0.31, 0.65])
ax3 = ax2.twinx()
ax4 = fig.add_axes([0.645, 0.155, 0.3, 0.65])

for ax in [ax1]:
    bp = ax.boxplot([data[data.Frame=="steel"].Time, data[data.Frame=="carbon"].Time], labels=["Steel\nframe", "Carbon\nframe"],
              boxprops = dict(linewidth=2, color=colors["lightgray"]),
              flierprops = dict(marker='o', markeredgecolor=colors["lightgray"], linestyle='none'),
              whiskerprops = dict(color=colors["lightgray"], linewidth=2),
              capprops = dict(color=colors["lightgray"], linewidth=2),
              medianprops = dict(color=colors["orange"], linewidth=2))
    for spine in ["bottom", "left"]:
        ax.spines[spine].set_linewidth(1)
        ax.spines[spine].set_color(colors["lightgray"])
    for spine in ["top", "right"]:
        ax.spines[spine].set_visible(False)
    ax.set_ylabel("Commuting time", size=18, color=colors["lightgray"], ha="center")
    ax.tick_params(axis="y", width=1, size=4, color=colors["lightgray"], labelcolor=colors["lightgray"], 
                   labelsize=13, pad=4)
    ax.tick_params(axis="x", width=1, size=0, labelcolor=colors["lightgray"], 
                   labelsize=13, pad=6)
    ## change color of the medians
    for color,median in zip([colors["orange"], colors["blue"]], bp['medians']):
        median.set(color=color, linewidth=2)
    ax.set_title("Boxplot\nof the samples", color=colors["lightgray"], size=16, ha="center")
        

for ax in [ax2]:
    ax.hist(data[data.Frame=="steel"].Time, bins="auto", color=colors["orange"], alpha=0.7, label="Steel")
    ax.hist(data[data.Frame=="carbon"].Time, bins="auto", color=colors["blue"], alpha=0.7, label="Carbon")
    for spine in ["bottom", "left"]:
        ax.spines[spine].set_linewidth(1)
        ax.spines[spine].set_color(colors["lightgray"])
    for spine in ["top", "right"]:
        ax.spines[spine].set_visible(False)
    ax.set_ylabel("Frequency", size=18, color=colors["lightgray"], ha="center")
    ax.tick_params(axis="y", width=1, size=4, color=colors["lightgray"], labelcolor=colors["lightgray"], 
                   labelsize=13, pad=4)
    ax.tick_params(axis="x", width=1, size=4, color=colors["lightgray"], labelcolor=colors["lightgray"], 
                   labelsize=13, pad=6)
    ax.set_xlabel("Commuting time", size=18, color=colors["lightgray"], ha="center")
    #ax.set_ylim(0)
    ax.set_xlim(data.Time.min()-7, data.Time.max()+7)
    legend = ax.legend(prop=dict(size=14), loc="upper right", frameon=False, facecolor="none",
                      bbox_to_anchor=(1, 1.05))
    for text in legend.get_texts():
        text.set_color(colors["lightgray"])
    ax.set_title("Underlying distributions\nof the populations the\nsamples are coming from", color=colors["lightgray"], size=16, ha="center")
    
for ax in [ax3]:
    ##norm distributions
    x_norm = np.linspace(ax.get_xlim()[0], ax.get_xlim()[1], 100)
    norm_steel = stats.norm(loc=data[data.Frame=="steel"].Time.mean(), scale=data[data.Frame=="steel"].Time.std())
    norm_carbon = stats.norm(loc=data[data.Frame=="carbon"].Time.mean(), scale=data[data.Frame=="carbon"].Time.std())
    y_norm_steel = norm_steel.pdf(x_norm)
    y_norm_carbon = norm_carbon.pdf(x_norm)
    ax.plot(x_norm, y_norm_steel, color=colors["orange"])
    ax.plot(x_norm, y_norm_carbon, color=colors["blue"])
    ax.set_ylim(0)
    ax.axis("off")
    
for ax in [ax4]:
    for spine in ["bottom"]:
        ax.spines[spine].set_linewidth(1)
        ax.spines[spine].set_color(colors["lightgray"])
    for spine in ["top", "right", "left"]:
        ax.spines[spine].set_visible(False)
    ax.set_yticks([])
    ax.tick_params(axis="x", width=1, size=4, color=colors["lightgray"], labelcolor=colors["lightgray"], 
                   labelsize=13, pad=4)
    ##t distributions
    bike_diff_mean_se = np.sqrt((data[data.Frame=="carbon"].Time.std()**2/len(data[data.Frame=="carbon"]))+(data[data.Frame=="steel"].Time.std()**2/len(data[data.Frame=="steel"])))
    x_t_diff = np.linspace(-6.5, 6.5, 100)
    t_diff = stats.t(df=np.min([len(data[data.Frame=="carbon"])-1, len(data[data.Frame=="steel"])-1]), loc=0, scale=bike_diff_mean_se)
    y_t_diff = t_diff.pdf(x_t_diff)
    ax.fill_between(x_t_diff, y_t_diff, color=colors["red"], alpha=0.8)
    ax.set_ylim(0)
    ax.set_xlabel(r"$\bar{x}_{\mathrm{carbon}}-\bar{x}_{\mathrm{steel}}$", size=18, color=colors["lightgray"], ha="center")
    ax.set_title("Sampling distribution of the difference\nin means of two samples from the\nunderlying distributions if $\mu_\mathrm{carbon}=\mu_\mathrm{steel}$", color=colors["lightgray"], size=16, ha="center")
    ax.text(2.5, ax.get_ylim()[1]*0.8, r"$\mu=0$", color=colors["lightgray"], size=13, ha="left")
    ax.text(2.5, ax.get_ylim()[1]*0.65, r"$\sigma=\sqrt{\frac{s_1^2}{n_1}+\frac{s_2^2}{n_2}}$", color=colors["lightgray"], size=13, ha="left")
    ax.text(7.75, ax.get_ylim()[1]*0.45, "df=$min(n_1-1;$" "\n" "$n_2-1)$", color=colors["lightgray"], size=13, ha="right")


    
plt.savefig("/Users/Gui/Box Sync/_STATS13/_Slides/_stats13-Lectures/assets/img/lec/bike-commute-ttest-model.svg", transparent=True)
fig = plt.figure(figsize=(12,4))
ax1 = fig.add_axes([0.065, 0.155, 0.2, 0.65])
ax2 = fig.add_axes([0.315, 0.155, 0.34, 0.65])
ax3 = ax2.twinx()
ax4 = fig.add_axes([0.675, 0.155, 0.3, 0.65])
    
for ax in [ax1]:
    bp = ax.boxplot([data[data.Frame=="steel"].Time, data[data.Frame=="carbon"].Time], labels=["Steel\nframe", "Carbon\nframe"],
              boxprops = dict(linewidth=2, color=colors["lightgray"]),
              flierprops = dict(marker='o', markeredgecolor=colors["lightgray"], linestyle='none'),
              whiskerprops = dict(color=colors["lightgray"], linewidth=2),
              capprops = dict(color=colors["lightgray"], linewidth=2),
              medianprops = dict(color=colors["orange"], linewidth=2))
    for spine in ["bottom", "left"]:
        ax.spines[spine].set_linewidth(1)
        ax.spines[spine].set_color(colors["lightgray"])
    for spine in ["top", "right"]:
        ax.spines[spine].set_visible(False)
    ax.set_ylabel("Commuting time", size=18, color=colors["lightgray"], ha="center")
    ax.tick_params(axis="y", width=1, size=4, color=colors["lightgray"], labelcolor=colors["lightgray"], 
                   labelsize=13, pad=4)
    ax.tick_params(axis="x", width=1, size=0, labelcolor=colors["lightgray"], 
                   labelsize=16, pad=6)
    ## change color of the medians
    for color,median in zip([colors["orange"], colors["blue"]], bp['medians']):
        median.set(color=color, linewidth=2)
    ax.set_title("Boxplot\nof the samples", color=colors["lightgray"], size=16, ha="center")
        

for ax in [ax2]:
    ax.hist([data[data.Frame=="steel"].Time, data[data.Frame=="carbon"].Time], bins=12, stacked=True, color=[colors["orange"], colors["blue"]], alpha=0.9, label=["Steel", "Carbon"])
    for spine in ["bottom", "left"]:
        ax.spines[spine].set_linewidth(1)
        ax.spines[spine].set_color(colors["lightgray"])
    for spine in ["top", "right"]:
        ax.spines[spine].set_visible(False)
    ax.set_ylabel("Frequency", size=18, color=colors["lightgray"], ha="center")
    ax.tick_params(axis="y", width=1, size=4, color=colors["lightgray"], labelcolor=colors["lightgray"], 
                   labelsize=13, pad=4)
    ax.tick_params(axis="x", width=1, size=4, color=colors["lightgray"], labelcolor=colors["lightgray"], 
                   labelsize=13, pad=6)
    ax.set_xlabel("Commuting time", size=18, color=colors["lightgray"], ha="center")
    #ax.set_ylim(0)
    ax.set_xlim(data.Time.min()-7, data.Time.max()+7)
    legend = ax.legend(prop=dict(size=14), loc="upper right", frameon=False, facecolor="none",
                      bbox_to_anchor=(1, 1.05))
    for text in legend.get_texts():
        text.set_color(colors["lightgray"])
    ax.set_title("Distribution of the\nunderlying population\nthe samples are coming from", color=colors["lightgray"], size=16, ha="center")
    
for ax in [ax3]:
    ##norm distributions
    x_norm = np.linspace(ax.get_xlim()[0], ax.get_xlim()[1], 100)
    norm_merged = stats.norm(loc=(data[data.Frame=="steel"].Time.mean()+data[data.Frame=="carbon"].Time.mean())/2, scale=data.Time.std())
    y_norm_merged = norm_steel.pdf(x_norm)
    ax.plot(x_norm, y_norm_merged, color=colors["lightgray"], lw=2)
    ax.set_ylim(0)
    ax.axis("off")
    
for ax in [ax4]:
    for spine in ["bottom"]:
        ax.spines[spine].set_linewidth(1)
        ax.spines[spine].set_color(colors["lightgray"])
    for spine in ["top", "right", "left"]:
        ax.spines[spine].set_visible(False)
    ax.set_yticks([])
    ax.tick_params(axis="x", width=1, size=4, color=colors["lightgray"], labelcolor=colors["lightgray"], 
                   labelsize=13, pad=4)
    ##t distributions
    bike_diff_mean_se = np.sqrt((data[data.Frame=="carbon"].Time.std()**2/len(data[data.Frame=="carbon"]))+(data[data.Frame=="steel"].Time.std()**2/len(data[data.Frame=="steel"])))
    x_t_diff = np.linspace(-6.5, 6.5, 100)
    t_diff = stats.t(df=len(data)-2, loc=0, scale=bike_diff_mean_se)
    y_t_diff = t_diff.pdf(x_t_diff)
    ax.fill_between(x_t_diff, y_t_diff, color=colors["red"], alpha=0.8)
    ax.set_ylim(0)
    ax.set_xlabel(r"$\bar{x}_{\mathrm{carbon}}-\bar{x}_{\mathrm{steel}}$", size=18, color=colors["lightgray"], ha="center")
    ax.set_title("Sampling distribution of\nthe difference in means of 2 samples\nfrom the distribution", color=colors["lightgray"], size=16, ha="center")
    ax.text(2.5, ax.get_ylim()[1]*0.8, r"$\mu=0$", color=colors["lightgray"], size=13, ha="left")
    ax.text(2.5, ax.get_ylim()[1]*0.65, r"$\sigma=\sqrt{\frac{s_1^2(n_1-1)+s_2^2(n_2-1)}{n_1+n_2-2}}$", color=colors["lightgray"], size=13, ha="left")
    ax.text(2.5, ax.get_ylim()[1]*0.5, "df=$(n_1+n_2)-2$", color=colors["lightgray"], size=13, ha="left")
    
#plt.savefig("/Users/Gui/Box Sync/_STATS13/_Slides/_stats13-Lectures/assets/img/lec/bike-commute-pooled-ttest-model.svg", transparent=True)
t_diff.sf(t_stat)*2
np.sqrt((6.25**2/26)+(4.89**2/30))<jupyter_output><empty_output><jupyter_text>Paired data<jupyter_code>data = pd.read_csv("data-Lock5DataCSV/Wetsuits.csv")
data
fig = plt.figure(figsize=(6,4))
ax1 = fig.add_axes([0.15, 0.15, 0.8, 0.8])
ax2 = ax1.twinx()
ax3 = ax1.twinx()

for ax in [ax1]:
    ax.plot(np.zeros(len(data)), data.NoWetsuit, "o", ms=10, mfc="none", mec=colors["orange"], mew=2)
    ax.plot(np.zeros(len(data))+1, data.Wetsuit, "o", ms=10, mfc="none", mec=colors["blue"], mew=2)
    for spine in ["bottom", "left"]:
        ax.spines[spine].set_linewidth(1)
        ax.spines[spine].set_color(colors["lightgray"])
    for spine in ["top", "right"]:
        ax.spines[spine].set_visible(False)
    ax.tick_params(axis="y", width=1, size=4, color=colors["lightgray"], labelcolor=colors["lightgray"], 
                   labelsize=13, pad=4)
    ax.tick_params(axis="x", width=1, size=0, labelcolor=colors["lightgray"], 
                   labelsize=16, pad=6)
    ax.set_ylabel("Maximum Velocity (m/s)", size=16, color=colors["lightgray"], ha="center")
    ax.set_xticks([0, 1])
    ax.set_xticklabels(["No Wetsuit", "Wetsuit"])
    ax.set_xlim(-0.4, 1.4)
    ax.set_ylim(1.05, 1.8)

for ax in [ax2]:
    for i,m in zip([0, 1], data[["NoWetsuit", "Wetsuit"]].mean()):
        ax.plot(i, m, "_", ms=50, color=colors["red"], mew=2)
    ax.set_ylim(ax1.get_ylim())
    ax.axis("off")

for ax in [ax3]:
    ax.axis("off")
    for row in range(len(data)):
        ax.plot([0, 1], data.ix[row][["NoWetsuit", "Wetsuit"]], "-", color=colors["lightgray"])
    ax.set_ylim(ax1.get_ylim())
    
plt.savefig("/Users/Gui/Box Sync/_STATS13/_Slides/_stats13-Lectures/assets/img/lec/wetsuit-data.svg", transparent=True)<jupyter_output><empty_output><jupyter_text>NHST - not paired<jupyter_code>n_simul = 10000

res_diff = np.zeros(n_simul)

population = data[["NoWetsuit", "Wetsuit"]].values.flatten()

for i in range(n_simul):
    np.random.shuffle(population)
    sample1 = population[:len(data.Wetsuit)]
    sample2 = population[len(data.Wetsuit):]
    res_diff[i] = np.mean(sample2)-np.mean(sample1)

fig = plt.figure(figsize=(6,4))
ax1 = fig.add_axes([0.01, 0.15, 0.82, 0.75])
ax2 = ax1.twinx()
ax3 = ax1.twinx()
    
for ax in [ax1]:
    ax.hist(res_diff, bins="auto", color=colors["blue"])
    for spine in ["bottom"]:
        ax.spines[spine].set_linewidth(1)
        ax.spines[spine].set_color(colors["lightgray"])
    for spine in ["top", "right", "left"]:
        ax.spines[spine].set_visible(False)
    ax.set_yticks([])
    ax.tick_params(axis="both", width=1, size=4, color=colors["lightgray"], labelcolor=colors["lightgray"], 
                   labelsize=13, pad=4)
    ax.set_xlabel(r'$\bar{x}_{\mathrm{Wetsuit}}-\bar{x}_{\mathrm{No Wetsuit}}$ (m/s)', size=18, color=colors["lightgray"])
    ax.set_ylim(0)

for ax in [ax2]:
    ax.set_ylim(ax1.get_ylim())
    ax.axvline(np.diff(data[["NoWetsuit", "Wetsuit"]].mean()), color=colors["orange"], ls="--")
    ax.text(np.diff(data[["NoWetsuit", "Wetsuit"]].mean()), ax.get_ylim()[1]*1.01, "Initial observation\n" "{:.3f}".format(np.diff(data[["NoWetsuit", "Wetsuit"]].mean())[0]), color=colors["orange"], size=13, ha="center")
    ax.axis("off")
    
for ax in [ax3]:
    ax.set_ylim(ax1.get_ylim())
    n_above = np.sum(res_diff>=np.diff(data[["NoWetsuit", "Wetsuit"]].mean()))
    ax.text(np.diff(data[["NoWetsuit", "Wetsuit"]].mean())+0.01, ax.get_ylim()[1]*0.7, "{:.0f} simulations with\n" r'$\bar{{x}}_{{\mathrm{{Wetsuit}}}}-\bar{{x}}_{{\mathrm{{No Wetsuit}}}}\geq$' "{:.3f}".format(n_above, np.diff(data[["NoWetsuit", "Wetsuit"]].mean())[0]), color=colors["lightgray"], size=13, ha="left")
    ax.text(np.diff(data[["NoWetsuit", "Wetsuit"]].mean())+0.01, ax.get_ylim()[1]*0.5, "1-tail p-value\np = {:.3f}".format(n_above/10000), color=colors["lightgray"], size=13, ha="left")
    ax.axis("off")
    
plt.savefig("/Users/Gui/Box Sync/_STATS13/_Slides/_stats13-Lectures/assets/img/lec/wetsuit-resampling-not-paired.svg", transparent=True)

#for next graph
lims = ax1.get_xlim()<jupyter_output><empty_output><jupyter_text>NHST - paired<jupyter_code>n_simul = 10000

res_diff_paired = np.zeros(n_simul)

for i in range(n_simul):
    labels = np.random.choice([-1, 1], size=len(data))
    res_diff_paired[i] = np.mean((data.Wetsuit-data.NoWetsuit)*labels)
fig = plt.figure(figsize=(6,4))
ax1 = fig.add_axes([0.01, 0.15, 0.82, 0.75])
ax2 = ax1.twinx()
ax3 = ax1.twinx()
    
for ax in [ax1]:
    ax.hist(res_diff_paired, bins="auto", color=colors["red"])
    for spine in ["bottom"]:
        ax.spines[spine].set_linewidth(1)
        ax.spines[spine].set_color(colors["lightgray"])
    for spine in ["top", "right", "left"]:
        ax.spines[spine].set_visible(False)
    ax.set_yticks([])
    ax.tick_params(axis="both", width=1, size=4, color=colors["lightgray"], labelcolor=colors["lightgray"], 
                   labelsize=13, pad=4)
    ax.set_xlabel(r'$\bar{x}_{\mathrm{Wetsuit}}-\bar{x}_{\mathrm{No Wetsuit}}$ (m/s)', size=18, color=colors["lightgray"])
    ax.set_ylim(0)
    ax.set_xlim(lims) #from previous graph
    
for ax in [ax2]:
    ax.set_ylim(ax1.get_ylim())
    ax.axvline(np.diff(data[["NoWetsuit", "Wetsuit"]].mean()), color=colors["orange"], ls="--")
    ax.text(np.diff(data[["NoWetsuit", "Wetsuit"]].mean()), ax.get_ylim()[1]*1.01, "Initial observation\n" "{:.3f}".format(np.diff(data[["NoWetsuit", "Wetsuit"]].mean())[0]), color=colors["orange"], size=13, ha="center")
    ax.axis("off")
    
for ax in [ax3]:
    ax.set_ylim(ax1.get_ylim())
    n_above = np.sum(res_diff_paired>=np.diff(data[["Wetsuit", "NoWetsuit"]].mean()))
    ax.text(np.diff(data[["NoWetsuit", "Wetsuit"]].mean())+0.01, ax.get_ylim()[1]*0.7, "{:.0f} simulations with\n" r'$\bar{{x}}_{{\mathrm{{Wetsuit}}}}-\bar{{x}}_{{\mathrm{{No Wetsuit}}}}\geq$' "{:.3f}".format(n_above, np.diff(data[["NoWetsuit", "Wetsuit"]].mean())[0]), color=colors["lightgray"], size=13, ha="left")
    ax.text(np.diff(data[["NoWetsuit", "Wetsuit"]].mean())+0.01, ax.get_ylim()[1]*0.5, "1-tail p-value\np<0.0001", color=colors["lightgray"], size=13, ha="left")
    ax.axis("off")
    
plt.savefig("/Users/Gui/Box Sync/_STATS13/_Slides/_stats13-Lectures/assets/img/lec/wetsuit-resampling-paired.svg", transparent=True)<jupyter_output><empty_output><jupyter_text>non paired vs paired analysis<jupyter_code>fig = plt.figure(figsize=(6,4))
ax1 = fig.add_axes([0.01, 0.15, 0.82, 0.75])
ax2 = ax1.twinx()
ax3 = ax1.twinx()
    
for ax in [ax1]:
    ax.hist(res_diff, bins="auto", color=colors["blue"], alpha=0.8)
    ax.hist(res_diff_paired, bins="auto", color=colors["red"])
    for spine in ["bottom"]:
        ax.spines[spine].set_linewidth(1)
        ax.spines[spine].set_color(colors["lightgray"])
    for spine in ["top", "right", "left"]:
        ax.spines[spine].set_visible(False)
    ax.set_yticks([])
    ax.tick_params(axis="both", width=1, size=4, color=colors["lightgray"], labelcolor=colors["lightgray"], 
                   labelsize=13, pad=4)
    ax.set_xlabel(r'$\bar{x}_{\mathrm{Wetsuit}}-\bar{x}_{\mathrm{No Wetsuit}}$ (m/s)', size=18, color=colors["lightgray"])
    ax.set_ylim(0)
    ax.text(ax.get_xlim()[0], ax.get_ylim()[1]*1.12, "Null hypothesis sampling\ndistribution when data\nare considered:", color=colors["lightgray"], size=13, ha="left", va="top")
    ax.text(ax.get_xlim()[0], ax.get_ylim()[1]*0.85, "- Independant", color=colors["blue"], size=15, ha="left")
    ax.text(ax.get_xlim()[0], ax.get_ylim()[1]*0.75, "- Paired", color=colors["red"], size=15, ha="left")
    ax.set_xlim(lims) #from previous graph
    
for ax in [ax2]:
    ax.set_ylim(ax1.get_ylim())
    ax.axvline(np.diff(data[["NoWetsuit", "Wetsuit"]].mean()), color=colors["orange"], ls="--")
    ax.text(np.diff(data[["NoWetsuit", "Wetsuit"]].mean()), ax.get_ylim()[1]*1.01, "Initial observation\n" "{:.3f}".format(np.diff(data[["NoWetsuit", "Wetsuit"]].mean())[0]), color=colors["orange"], size=13, ha="center")
    ax.axis("off")
    
for ax in [ax3]:
    ax.set_ylim(ax1.get_ylim())
    n_above = np.sum(res_diff_paired>=np.diff(data[["NoWetsuit", "Wetsuit"]].mean()))
    ax.text(np.diff(data[["NoWetsuit", "Wetsuit"]].mean())+0.01, ax.get_ylim()[1]*0.7, "{:.0f} simulations with\n" r'$\bar{{x}}_{{\mathrm{{Wetsuit}}}}-\bar{{x}}_{{\mathrm{{No Wetsuit}}}}\geq$' "{:.3f}".format(n_above, np.diff(data[["NoWetsuit", "Wetsuit"]].mean())[0]), color=colors["lightgray"], size=13, ha="left")
    ax.text(np.diff(data[["NoWetsuit", "Wetsuit"]].mean())+0.01, ax.get_ylim()[1]*0.5, "1-tail p-value\np<0.0001", color=colors["lightgray"], size=13, ha="left")
    ax.axis("off")
    
plt.savefig("/Users/Gui/Box Sync/_STATS13/_Slides/_stats13-Lectures/assets/img/lec/wetsuit-resampling-paired-vs-non-paired.svg", transparent=True)<jupyter_output><empty_output><jupyter_text>95% confidence intervals<jupyter_code>diffs = data.Wetsuit-data.NoWetsuit


n_simul = 10000
res_diff_ci = np.zeros(n_simul)

for i in range(n_simul):
    bootsample = np.random.choice(diffs, size=len(diffs))
    res_diff_ci[i] = np.mean(bootsample)

np.percentile(res_diff_ci, [2.5, 97.5])
#Observed diff
observed_diff = np.mean(diffs)

fig = plt.figure(figsize=(6,4))
ax1 = fig.add_axes([0.1, 0.15, 0.82, 0.75])
ax2 = ax1.twinx()
ax3 = ax1.twinx()
    
for ax in [ax1]:
    ax.hist(res_diff_ci, bins="auto", color=colors["blue"])
    for spine in ["bottom"]:
        ax.spines[spine].set_linewidth(1)
        ax.spines[spine].set_color(colors["lightgray"])
    for spine in ["top", "right", "left"]:
        ax.spines[spine].set_visible(False)
    ax.set_yticks([])
    ax.tick_params(axis="both", width=1, size=4, color=colors["lightgray"], labelcolor=colors["lightgray"], 
                   labelsize=13, pad=4)
    ax.set_xlabel(r'$\bar{x}_{\mathrm{Wetsuit}}-\bar{x}_{\mathrm{No Wetsuit}}$', size=18, color=colors["lightgray"])
    #ax.set_xlim(-7, 7)
    ax.set_ylim(0)
    
for ax in [ax2]:
    ax.set_ylim(ax1.get_ylim())
    ax.axvline(observed_diff, color=colors["orange"], ls="--")
    ax.text(observed_diff, ax.get_ylim()[1]*1.01, "Initial observation\n" r"$\bar{x}_{\mathrm{Wetsuit}}-\bar{x}_{\mathrm{No Wetsuit}}=0.078$", color=colors["orange"], size=13, ha="center")
    ax.axis("off")

for ax in [ax3]:
    ax.axvline(np.percentile(res_diff_ci, 2.5), ymax=0.55, color=colors["red"], lw=2)
    ax.axvline(np.percentile(res_diff_ci, 97.5), ymax=0.55, color=colors["red"], lw=2)
    ax.text(np.percentile(res_diff_ci, 2.5), ax.get_ylim()[1]*0.6, "2.5$^{{th}}$\npercentile\n{:.3f}".format(np.percentile(res_diff_ci, 2.5)), color=colors["red"], size=15, ha="center")
    ax.text(np.percentile(res_diff_ci, 97.5), ax.get_ylim()[1]*0.6, "97.5$^{{th}}$\npercentile\n{:.3f}".format(np.percentile(res_diff_ci, 97.5)), color=colors["red"], size=15, ha="center")
    ax.axis("off")


plt.savefig("/Users/Gui/Box Sync/_STATS13/_Slides/_stats13-Lectures/assets/img/lec/wetsuit-resampling-paired-ci95.svg", transparent=True)
np.diff(data[["Wetsuit", "NoWetsuit"]].mean())[0]
plt.plot(data.Wetsuit, data.NoWetsuit, "o")
data.describe()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.21071409118143689, 'fraction_numerical': 0.03870036494362438, 'mean_word_length': 4.678781317661614, 'pattern_counts': {'":': 0, '<': 54, '<?xml version=': 0, '>': 55, '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': '11', 'dir': 'jupyter-structured-clean-dedup', 'id': '648360', 'license': 'MIT', 'n_tokens_mistral': 15427, 'n_tokens_neox': 14626, 'n_words': 2073, 'path': 'notebook-fig/figs-slides-6-comparing-two-means.ipynb', 'repo_name': 'gcalmettes/stats13-slides', 'size': '36718'}