Statistical Thinking in Python II.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Parameter estimation by optimization"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "** How often do we get no-hitters?**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Seed random number generator\n",
    "np.random.seed(42)\n",
    "\n",
    "# Compute mean no-hitter time: tau\n",
    "tau = np.mean(nohitter_times)\n",
    "\n",
    "# Draw out of an exponential distribution with parameter tau: inter_nohitter_time\n",
    "inter_nohitter_time = np.random.exponential(tau, 100000)\n",
    "\n",
    "# Plot the PDF and label axes\n",
    "_ = plt.hist(inter_nohitter_time,\n",
    "             bins=50, normed=True, histtype='step')\n",
    "_ = plt.xlabel('Games between no-hitters')\n",
    "_ = plt.ylabel('PDF')\n",
    "\n",
    "# Show the plot\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Overlay the theoretical CDF with the ECDF from the data**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It looks like no-hitters in the modern era of Major League Baseball are Exponentially distributed. Based on the story of the Exponential distribution, this suggests that they are a random process; when a no-hitter will happen is independent of when the last no-hitter was."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Create an ECDF from real data: x, y\n",
    "x, y = ecdf(nohitter_times)\n",
    "\n",
    "# Create a CDF from theoretical samples: x_theor, y_theor\n",
    "x_theor, y_theor = ecdf(inter_nohitter_time)\n",
    "\n",
    "# Overlay the plots\n",
    "plt.plot(x_theor, y_theor)\n",
    "plt.plot(x, y, marker='.', linestyle='none')\n",
    "\n",
    "# Margins and axis labels\n",
    "plt.margins(0.02)\n",
    "plt.xlabel('Games between no-hitters')\n",
    "plt.ylabel('CDF')\n",
    "\n",
    "# Show the plot\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "** Choosing the optimal parameter **"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Plot the theoretical CDFs\n",
    "plt.plot(x_theor, y_theor)\n",
    "plt.plot(x, y, marker='.', linestyle='none')\n",
    "plt.margins(0.02)\n",
    "plt.xlabel('Games between no-hitters')\n",
    "plt.ylabel('CDF')\n",
    "\n",
    "# Take samples with half tau: samples_half\n",
    "samples_half = np.random.exponential(tau/2,10000)\n",
    "\n",
    "# Take samples with double tau: samples_double\n",
    "samples_double = np.random.exponential(tau*2,10000)\n",
    "\n",
    "# Generate CDFs from these samples\n",
    "x_half, y_half = ecdf(samples_half)\n",
    "x_double, y_double = ecdf(samples_double)\n",
    "\n",
    "# Plot these CDFs as lines\n",
    "_ = plt.plot(x_half, y_half)\n",
    "_ = plt.plot(x_double, y_double)\n",
    "\n",
    "# Show the plot\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**EDA of literacy/fertility data**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Plot the illiteracy rate versus fertility\n",
    "_ = plt.plot(illiteracy,fertility, marker='.', linestyle='none')\n",
    "\n",
    "# Set the margins and label axes\n",
    "plt.margins(0.02)\n",
    "_ = plt.xlabel('percent illiterate')\n",
    "_ = plt.ylabel('fertility')\n",
    "\n",
    "# Show the plot\n",
    "plt.show()\n",
    "\n",
    "# Show the Pearson correlation coefficient\n",
    "print(pearson_r(illiteracy, fertility))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Linear regression**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The slope tells us how the fertility rate varies with illiteracy. We can find the best fit line using np.polyfit()."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Plot the illiteracy rate versus fertility\n",
    "_ = plt.plot(illiteracy, fertility, marker='.', linestyle='none')\n",
    "plt.margins(0.02)\n",
    "_ = plt.xlabel('percent illiterate')\n",
    "_ = plt.ylabel('fertility')\n",
    "\n",
    "# Perform a linear regression using np.polyfit(): a, b\n",
    "a, b = np.polyfit(illiteracy, fertility,1)\n",
    "\n",
    "# Print the results to the screen\n",
    "print('slope =', a, 'children per woman / percent illiterate')\n",
    "print('intercept =', b, 'children per woman')\n",
    "\n",
    "# Make theoretical line to plot\n",
    "x = np.array([0,100])\n",
    "y = a * x + b\n",
    "\n",
    "# Add regression line to your plot\n",
    "_ = plt.plot(x, y)\n",
    "_ = plt.xlabel('illeteracy')\n",
    "_ = plt.ylabel('fertility')\n",
    "\n",
    "# Draw the plot\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**How is it optimal?**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice that the minimum on the plot, that is the value of the slope that gives the minimum sum of the square of the residuals, is the same value you got when performing the regression."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Specify slopes to consider: a_vals\n",
    "a_vals = np.linspace(0, 0.1, 200)\n",
    "\n",
    "# Initialize sum of square of residuals: rss\n",
    "rss = np.empty_like(a_vals)\n",
    "\n",
    "# Compute sum of square of residuals for each value of a_vals\n",
    "for i, a in enumerate(a_vals):\n",
    "    rss[i] = np.sum((fertility - a*illiteracy - b)**2)\n",
    "\n",
    "# Plot the RSS\n",
    "plt.plot(a_vals, rss, '-')\n",
    "plt.xlabel('slope (children per woman / percent illiterate)')\n",
    "plt.ylabel('sum of square of residuals')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Linear regression on appropriate Anscombe data**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Perform linear regression: a, b\n",
    "a, b = np.polyfit(x,y,deg=1)\n",
    "\n",
    "# Print the slope and intercept\n",
    "print(a, b)\n",
    "\n",
    "# Generate theoretical x and y data: x_theor, y_theor\n",
    "x_theor = np.array([3, 15])\n",
    "y_theor = a * x_theor + b\n",
    "\n",
    "# Plot the Anscombe data and theoretical line\n",
    "_ = plt.plot(x,y,marker='.', linestyle='none')\n",
    "_ = plt.plot(x_theor, y_theor)\n",
    "\n",
    "# Label the axes\n",
    "plt.xlabel('x')\n",
    "plt.ylabel('y')\n",
    "\n",
    "# Show the plot\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Linear regression on all Anscombe data**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Iterate through x,y pairs\n",
    "for x, y in zip(anscombe_x,anscombe_y):\n",
    "    # Compute the slope and intercept: a, b\n",
    "    a, b = np.polyfit(x,y,deg=1)\n",
    "\n",
    "    # Print the result\n",
    "    print('slope:', a, 'intercept:', b)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Bootstrap confidence intervals"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Bootstrap is a replicate data generated by randomly selecting similar data to the original dataset."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "**Visualizing bootstrap samples **"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "for x in range(50):\n",
    "    # Generate bootstrap sample: bs_sample\n",
    "    bs_sample = np.random.choice(rainfall, size=len(rainfall))\n",
    "\n",
    "    # Compute and plot ECDF from bootstrap sample\n",
    "    x, y = ecdf(bs_sample)\n",
    "    _ = plt.plot(x, y, marker='.', linestyle='none',\n",
    "                 color='gray', alpha=0.1)\n",
    "\n",
    "# Compute and plot ECDF from original data\n",
    "x, y = ecdf(rainfall)\n",
    "_ = plt.plot(x, y, marker='.')\n",
    "\n",
    "# Make margins and label axes\n",
    "plt.margins(0.02)\n",
    "_ = plt.xlabel('yearly rainfall (mm)')\n",
    "_ = plt.ylabel('ECDF')\n",
    "\n",
    "# Show the plot\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "** Generating many bootstrap replicates**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def bootstrap_replicate_1d(data, func):\n",
    "    return func(np.random.choice(data, size=len(data)))\n",
    "\n",
    "def draw_bs_reps(data, func, size=1):\n",
    "    \"\"\"Draw bootstrap replicates.\"\"\"\n",
    "\n",
    "    # Initialize array of replicates: bs_replicates\n",
    "    bs_replicates = np.empty(size)\n",
    "\n",
    "    # Generate replicates\n",
    "    for i in range(size):\n",
    "        bs_replicates[i] = bootstrap_replicate_1d(data, func)\n",
    "\n",
    "    return bs_replicates"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Bootstrap replicates of the mean and the SEM**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Take 10,000 bootstrap replicates of the mean: bs_replicates\n",
    "bs_replicates = draw_bs_reps(rainfall, np.mean, 10000)\n",
    "\n",
    "# Compute and print SEM\n",
    "sem = np.std(rainfall)/ np.sqrt(len(rainfall))\n",
    "print(sem)\n",
    "\n",
    "# Compute and print standard deviation of bootstrap replicates\n",
    "bs_std = np.std(bs_replicates)\n",
    "print(bs_std)\n",
    "\n",
    "# Make a histogram of the results\n",
    "_ = plt.hist(bs_replicates, bins=50, normed=True)\n",
    "_ = plt.xlabel('mean annual rainfall (mm)')\n",
    "_ = plt.ylabel('PDF')\n",
    "\n",
    "# Show the plot\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "### computing confidence interval\n",
    "conf_interval = np.percentile(bs_replicates, [2.5, 97.5]) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "** Bootstrap replicates of other statistics** "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Generate 10,000 bootstrap replicates of the variance: bs_replicates\n",
    "bs_replicates = draw_bs_reps(rainfall, np.var, 10000)\n",
    "\n",
    "# Put the variance in units of square centimeters\n",
    "bs_replicates= bs_replicates/100\n",
    "\n",
    "# Make a histogram of the results\n",
    "_ = plt.hist(bs_replicates, normed=True, bins=50)\n",
    "_ = plt.xlabel('variance of annual rainfall (sq. cm)')\n",
    "_ = plt.ylabel('PDF')\n",
    "\n",
    "# Show the plot\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**A function to do pairs bootstrap **"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As discussed in the video, pairs bootstrap involves resampling pairs of data. Each collection of pairs fit with a line, in this case using np.polyfit(). We do this again and again, getting bootstrap replicates of the parameter values. To have a useful tool for doing pairs bootstrap, you will write a function to perform pairs bootstrap on a set of x,y data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def draw_bs_pairs_linreg(x, y, size=1):\n",
    "    \"\"\"Perform pairs bootstrap for linear regression.\"\"\"\n",
    "\n",
    "    # Set up array of indices to sample from: inds\n",
    "    inds = np.arange(len(x))\n",
    "\n",
    "    # Initialize replicates: bs_slope_reps, bs_intercept_reps\n",
    "    bs_slope_reps = np.empty(size)\n",
    "    bs_intercept_reps = np.empty(size)\n",
    "\n",
    "    # Generate replicates\n",
    "    for i in range(size):\n",
    "        bs_inds = np.random.choice(inds, size=len(inds))\n",
    "        bs_x, bs_y = x[bs_inds], y[bs_inds]\n",
    "        bs_slope_reps[i], bs_intercept_reps[i] = np.polyfit(bs_x, bs_y, 1)\n",
    "\n",
    "    return bs_slope_reps, bs_intercept_reps"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Pairs bootstrap of literacy/fertility data**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Generate replicates of slope and intercept using pairs bootstrap\n",
    "bs_slope_reps, bs_intercept_reps = draw_bs_pairs_linreg(illiteracy, fertility, size=1000)\n",
    "\n",
    "# Compute and print 95% CI for slope\n",
    "print(np.percentile(bs_slope_reps, [2.5, 97.5]))\n",
    "\n",
    "# Plot the histogram\n",
    "_ = plt.hist(bs_slope_reps, bins=50, normed=True)\n",
    "_ = plt.xlabel('slope')\n",
    "_ = plt.ylabel('PDF')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Plotting bootstrap regressions**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Generate array of x-values for bootstrap lines: x\n",
    "x = np.array([0,100])\n",
    "\n",
    "# Plot the bootstrap lines\n",
    "for i in range(100):\n",
    "    _ = plt.plot(x, bs_slope_reps[i]*x + bs_intercept_reps[i],\n",
    "                 linewidth=0.5, alpha=0.2, color='red')\n",
    "\n",
    "# Plot the data\n",
    "_ = plt.plot()\n",
    "\n",
    "# Label axes, set the margins, and show the plot\n",
    "_ = plt.xlabel('illiteracy')\n",
    "_ = plt.ylabel('fertility')\n",
    "plt.margins(0.02)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Hypothesis testing "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "** Generating a permutation sample **"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def permutation_sample(data1, data2):\n",
    "    \"\"\"Generate a permutation sample from two data sets.\"\"\"\n",
    "\n",
    "    # Concatenate the data sets: data\n",
    "    data = np.concatenate((data1, data2))\n",
    "\n",
    "    # Permute the concatenated array: permuted_data\n",
    "    permuted_data = np.random.permutation(data)\n",
    "\n",
    "    # Split the permuted array into two: perm_sample_1, perm_sample_2\n",
    "    perm_sample_1 = permuted_data[:len(data1)]\n",
    "    perm_sample_2 = permuted_data[len(data1):]\n",
    "\n",
    "    return perm_sample_1, perm_sample_2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "** Vizualizing permutation samplting **"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "for x in range(50):\n",
    "    # Generate permutation samples\n",
    "    perm_sample_1, perm_sample_2 = permutation_sample(rain_july, rain_november)\n",
    "\n",
    "\n",
    "    # Compute ECDFs\n",
    "    x_1, y_1 = ecdf(perm_sample_1)\n",
    "    x_2, y_2 = ecdf(perm_sample_2)\n",
    "\n",
    "    # Plot ECDFs of permutation sample\n",
    "    _ = plt.plot(x_1, y_1, marker='.', linestyle='none',\n",
    "                 color='red', alpha=0.02)\n",
    "    _ = plt.plot(x_2, y_2, marker='.', linestyle='none',\n",
    "                 color='blue', alpha=0.02)\n",
    "\n",
    "# Create and plot ECDFs from original data\n",
    "x_1, y_1 = ecdf(rain_july)\n",
    "x_2, y_2 = ecdf(rain_november)\n",
    "_ = plt.plot(x_1, y_1, marker='.', linestyle='none', color='red')\n",
    "_ = plt.plot(x_2, y_2, marker='.', linestyle='none', color='blue')\n",
    "\n",
    "# Label axes, set margin, and show plot\n",
    "plt.margins(0.02)\n",
    "_ = plt.xlabel('monthly rainfall (mm)')\n",
    "_ = plt.ylabel('ECDF')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Generating permutation replicates**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def draw_perm_reps(data_1, data_2, func, size=1):\n",
    "    \"\"\"Generate multiple permutation replicates.\"\"\"\n",
    "\n",
    "    # Initialize array of replicates: perm_replicates\n",
    "    perm_replicates = np.empty(size)\n",
    "\n",
    "    for i in range(size):\n",
    "        # Generate permutation sample\n",
    "        perm_sample_1, perm_sample_2 = permutation_sample(data_1, data_2)\n",
    "\n",
    "        # Compute the test statistic\n",
    "        perm_replicates[i] = func(perm_sample_1, perm_sample_2)\n",
    "\n",
    "    return perm_replicates"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "** Permutation test on frog data **"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The average strike force of Frog A was 0.71 Newtons (N), and that of Frog B was 0.42 N for a difference of 0.29 N. It is possible the frogs strike with the same force and this observed difference was by chance. You will compute 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.\n",
    "\n",
    "For your convenience, the data has been stored in the arrays force_a and force_b."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def diff_of_means(data_1, data_2):\n",
    "    \"\"\"Difference in means of two arrays.\"\"\"\n",
    "\n",
    "    # The difference of means of data_1, data_2: diff\n",
    "    diff = np.mean(data_1) - np.mean(data_2)\n",
    "\n",
    "    return diff\n",
    "\n",
    "# Compute difference of mean impact force from experiment: empirical_diff_means\n",
    "empirical_diff_means = diff_of_means(force_a, force_b)\n",
    "\n",
    "# Draw 10,000 permutation replicates: perm_replicates\n",
    "perm_replicates = draw_perm_reps(force_a, force_b,\n",
    "                                 diff_of_means, size=10000)\n",
    "\n",
    "# Compute p-value: p\n",
    "p = np.sum(perm_replicates >= empirical_diff_means) / len(perm_replicates)\n",
    "\n",
    "# Print the result\n",
    "print('p-value =', p)\n",
    "\n",
    "'''\n",
    "The p-value tells you that there is about a 0.6% chance that you would get the difference of means\n",
    "observed in the experiment if frogs were exactly the same. \n",
    "A p-value below 0.01 is typically said to be \"statistically significannt,\", \n",
    "but: warning! warning! warning! You have computed a p-value; it is a number. \n",
    "I encourage you not to distill it to a yes-or-no phrase. \n",
    "p = 0.006 and p = 0.000000006 are both said to be \"statistically significant,\" \n",
    "but they are definitely not the same!\n",
    "'''"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "** Bootstrap hypothesis test **"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "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.\n",
    "\n",
    "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."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Make an array of translated impact forces: translated_force_b\n",
    "translated_force_b = (force_b + 0.55) - np.mean(force_b)\n",
    "\n",
    "# Take bootstrap replicates of Frog B's translated impact forces: bs_replicates\n",
    "bs_replicates = draw_bs_reps(translated_force_b, np.mean, 10000)\n",
    "\n",
    "# Compute fraction of replicates that are less than the observed Frog B force: p\n",
    "p = np.sum(bs_replicates <= np.mean(force_b)) / 10000\n",
    "\n",
    "# Print the p-value\n",
    "print('p = ', p)\n",
    "\n",
    "'''\n",
    " The low p-value suggests that the null hypothesis that Frog B and Frog C have the same mean impact force is false.\n",
    "'''"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**A bootstrap test for identical distributions**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can even test the same hypothesis that we tested with a permutation test: that the Frog A and Frog B have identically distributed impact forces. 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.\n",
    "\n",
    "Now, you will 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.\n",
    "\n",
    "The two arrays are available to you as force_a and force_b."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Compute difference of mean impact force from experiment: empirical_diff_means\n",
    "empirical_diff_means = diff_of_means(force_a, force_b)\n",
    "\n",
    "# Concatenate forces: forces_concat\n",
    "forces_concat = np.concatenate((force_a, force_b))\n",
    "\n",
    "# Initialize bootstrap replicates: bs_replicates\n",
    "bs_replicates = np.empty(10000)\n",
    "\n",
    "for i in range(10000):\n",
    "    # Generate bootstrap sample\n",
    "    bs_sample = np.random.choice(forces_concat, size=len(forces_concat))\n",
    "    \n",
    "    # Compute replicate\n",
    "    bs_replicates[i] = diff_of_means(bs_sample[:len(force_a)],\n",
    "                                     bs_sample[len(force_a):])\n",
    "\n",
    "# Compute and print p-value: p\n",
    "p = np.sum(bs_replicates > empirical_diff_means) / 10000\n",
    "print('p-value =', p)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "** A two-sample bootstrap hypothesis test for difference of means.**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "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.\n",
    "\n",
    "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.\n",
    "\n",
    "The objects forces_concat and empirical_diff_means are already in your namespace."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Compute mean of all forces: mean_force\n",
    "mean_force = np.mean(forces_concat)\n",
    "\n",
    "# Generate shifted arrays\n",
    "force_a_shifted = force_a - np.mean(force_a) + mean_force\n",
    "force_b_shifted = force_b - np.mean(force_b) + mean_force\n",
    "\n",
    "# Compute 10,000 bootstrap replicates from shifted arrays\n",
    "bs_replicates_a = draw_bs_reps(force_a_shifted, np.mean, 10000)\n",
    "bs_replicates_b = draw_bs_reps(force_b_shifted, np.mean, 10000)\n",
    "\n",
    "# Get replicates of difference of means: bs_replicates\n",
    "bs_replicates = bs_replicates_a - bs_replicates_b\n",
    "\n",
    "# Compute and print p-value: p\n",
    "p = np.sum(bs_replicates > empirical_diff_means) / 10000\n",
    "print('p-value =', p)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Hypothesis test examples "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "** A/B testing - The vote for the Civil Rights Act in 1964**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Construct arrays of data: dems, reps\n",
    "dems = np.array([True] * 153 + [False] * 91)\n",
    "reps = np.array([True] * 136 + [False] * 35)\n",
    "\n",
    "def frac_yay_dems(dems, reps):\n",
    "    \"\"\"Compute fraction of Democrat yay votes.\"\"\"\n",
    "    frac = np.sum(dems) / len(dems)\n",
    "    return frac\n",
    "\n",
    "# Acquire permutation samples: perm_replicates\n",
    "perm_replicates = draw_perm_reps(dems, reps, frac_yay_dems, 10000)\n",
    "\n",
    "# Compute and print p-value: p\n",
    "p = np.sum(perm_replicates <= 153/244) / len(perm_replicates)\n",
    "print('p-value =', p)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**A/B testing - time on website analog  **"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Compute the observed difference in mean inter-no-hitter times: nht_diff_obs\n",
    "nht_diff_obs = diff_of_means(nht_dead, nht_live)\n",
    "\n",
    "# Acquire 10,000 permutation replicates of difference in mean no-hitter time: perm_replicates\n",
    "perm_replicates = draw_perm_reps(nht_dead, nht_live, diff_of_means, 10000)\n",
    "\n",
    "\n",
    "# Compute and print the p-value: p\n",
    "p = np.sum(perm_replicates<=nht_diff_obs)/len(perm_replicates)\n",
    "print('p-val =',p)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Simulating a null hypothesis concerning correlation**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The observed correlation between female illiteracy and fertility in the data set of 162 countries may just be by chance; the fertility of a given country may actually be totally independent of its illiteracy. You will test this null hypothesis in the next exercise.\n",
    "\n",
    "To do the test, you need to simulate the data assuming the null hypothesis is true. Of the following choices, which is the best way to to do it?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Do a permutation test: Permute the illiteracy values, \n",
    "# but leave the fertility values fixed to generate a new set of (illiteracy, fertility) data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Hypothesis test on Pearson correlation**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The observed correlation between female illiteracy and fertility may just be by chance; the fertility of a given country may actually be totally independent of its illiteracy. You will test this hypothesis. To do so, 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.\n",
    "\n",
    "The function pearson_r() that you wrote in the prequel to this course for computing the Pearson correlation coefficient is already in your name space."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Compute observed correlation: r_obs\n",
    "r_obs = pearson_r(illiteracy, fertility)\n",
    "\n",
    "# Initialize permutation replicates: perm_replicates\n",
    "perm_replicates = np.empty(10000)\n",
    "\n",
    "# Draw replicates\n",
    "for i in range(10000):\n",
    "    # Permute illiteracy measurments: illiteracy_permuted\n",
    "    illiteracy_permuted = np.random.permutation(illiteracy)\n",
    "\n",
    "    # Compute Pearson correlation\n",
    "    perm_replicates[i] = pearson_r(illiteracy_permuted, fertility)\n",
    "\n",
    "# Compute p-value: p\n",
    "p = np.sum(perm_replicates >= r_obs) / len(perm_replicates)\n",
    "print('p-val =', p)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Do neonicotinoid insecticides have unintended consequences?**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As a final exercise in hypothesis testing before we put everything together in our case study in the next chapter, you will investigate the effects of neonicotinoid insecticides on bee reproduction. These insecticides are very widely used in the United States to combat aphids and other pests that damage plants.\n",
    "\n",
    "In a recent study, Straub, et al. (Proc. Roy. Soc. B, 2016) investigated the effects of neonicotinoids on the sperm of pollinating bees. In this and the next exercise, you will study how the pesticide treatment affected the count of live sperm per half milliliter of semen.\n",
    "\n",
    "First, we will do EDA, as usual. Plot ECDFs of the alive sperm count for untreated bees (stored in the Numpy array control) and bees treated with pesticide (stored in the Numpy array treated)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Compute x,y values for ECDFs\n",
    "x_control, y_control = ecdf(control)\n",
    "x_treated, y_treated = ecdf(treated)\n",
    "\n",
    "# Plot the ECDFs\n",
    "plt.plot(x_control, y_control, marker='.', linestyle='none')\n",
    "plt.plot(x_treated, y_treated, marker='.', linestyle='none')\n",
    "\n",
    "# Set the margins\n",
    "plt.margins(0.02)\n",
    "\n",
    "# Add a legend\n",
    "plt.legend(('control', 'treated'), loc='lower right')\n",
    "\n",
    "# Label axes and show plot\n",
    "plt.xlabel('millions of alive sperm per mL')\n",
    "plt.ylabel('ECDF')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Bootstrap hypothesis test on bee sperm counts**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, you will test 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. You will use the difference of means as your test statistic.\n",
    "\n",
    "For your reference, the call signature for the draw_bs_reps() function you wrote in chapter 2 is draw_bs_reps(data, func, size=1)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Compute the difference in mean sperm count: diff_means\n",
    "diff_means = diff_of_means(control, treated)\n",
    "\n",
    "# Compute mean of pooled data: mean_count\n",
    "mean_count = np.mean(np.concatenate((control, treated)))\n",
    "\n",
    "# Generate shifted data sets\n",
    "control_shifted = control - np.mean(control) + mean_count\n",
    "treated_shifted = treated - np.mean(treated) + mean_count\n",
    "\n",
    "# Generate bootstrap replicates\n",
    "bs_reps_control = draw_bs_reps(control_shifted,\n",
    "                       np.mean, size=10000)\n",
    "bs_reps_treated = draw_bs_reps(treated_shifted,\n",
    "                       np.mean, size=10000)\n",
    "\n",
    "# Get replicates of difference of means: bs_replicates\n",
    "bs_replicates = bs_reps_control-bs_reps_treated\n",
    "\n",
    "# Compute and print p-value: p\n",
    "p = np.sum(bs_replicates >= np.mean(control) - np.mean(treated)) \\\n",
    "            / len(bs_replicates)\n",
    "print('p-value =', p)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Case study"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "** EDA **"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Create bee swarm plot\n",
    "_ = sns.swarmplot(x='year', y='beak_depth', data=df)\n",
    "\n",
    "# Label the axes\n",
    "_ = plt.xlabel('year')\n",
    "_ = plt.ylabel('beak depth (mm)')\n",
    "\n",
    "# Show the plot\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "** ECDF of beak depths **"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Compute ECDFs\n",
    "x_1975, y_1975 = ecdf(bd_1975)\n",
    "x_2012, y_2012 = ecdf(bd_2012)\n",
    "\n",
    "# Plot the ECDFs\n",
    "_ = plt.plot(x_1975, y_1975, marker='.', linestyle='none')\n",
    "_ = plt.plot(x_2012, y_2012, marker='.', linestyle='none')\n",
    "\n",
    "# Set margins\n",
    "plt.margins(0.02)\n",
    "\n",
    "# Add axis labels and legend\n",
    "_ = plt.ylabel('beak depth (mm)')\n",
    "_ = plt.xlabel('ECDF')\n",
    "_ = plt.legend(('1975', '2012'), loc='lower right')\n",
    "\n",
    "# Show the plot\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**ECDF of beak depths**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Compute ECDFs\n",
    "x_1975, y_1975 = ecdf(bd_1975)\n",
    "x_2012, y_2012 = ecdf(bd_2012)\n",
    "\n",
    "# Plot the ECDFs\n",
    "_ = plt.plot(x_1975, y_1975, marker='.', linestyle='none')\n",
    "_ = plt.plot(x_2012, y_2012, marker='.', linestyle='none')\n",
    "\n",
    "# Set margins\n",
    "plt.margins(0.02)\n",
    "\n",
    "# Add axis labels and legend\n",
    "_ = plt.ylabel('ECDF')\n",
    "_ = plt.xlabel('beak depth (mm)')\n",
    "_ = plt.legend(('1975', '2012'), loc='lower right')\n",
    "\n",
    "# Show the plot\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "**Parameter estimates of beak depths**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Compute the difference of the sample means: mean_diff\n",
    "mean_diff = np.mean(bd_2012) - np.mean(bd_1975)\n",
    "\n",
    "# Get bootstrap replicates of means\n",
    "bs_replicates_1975 = draw_bs_reps(bd_1975, np.mean, 10000)\n",
    "bs_replicates_2012 = draw_bs_reps(bd_2012, np.mean, 10000)\n",
    "\n",
    "# Compute samples of difference of means: bs_diff_replicates\n",
    "bs_diff_replicates = bs_replicates_2012 - bs_replicates_1975 \n",
    "# Compute 95% confidence interval: conf_int\n",
    "conf_int = np.percentile(bs_diff_replicates, [2.5, 97.5]) \n",
    "\n",
    "# Print the results\n",
    "print('difference of means =', mean_diff, 'mm')\n",
    "print('95% confidence interval =', conf_int, 'mm')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Hypothesis test: Are beaks deeper in 2012?**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Compute mean of combined data set: combined_mean\n",
    "combined_mean = np.mean(np.concatenate((bd_1975, bd_2012)))\n",
    "\n",
    "# Shift the samples\n",
    "bd_1975_shifted = bd_1975 - np.mean(bd_1975) + combined_mean\n",
    "bd_2012_shifted = bd_2012 - np.mean(bd_2012) + combined_mean\n",
    "\n",
    "# Get bootstrap replicates of shifted data sets\n",
    "bs_replicates_1975 = draw_bs_reps(bd_1975_shifted, np.mean, 10000)\n",
    "bs_replicates_2012 = draw_bs_reps(bd_2012_shifted, np.mean, 10000)\n",
    "\n",
    "\n",
    "# Compute replicates of difference of means: bs_diff_replicates\n",
    "bs_diff_replicates = bs_replicates_2012 - bs_replicates_1975\n",
    "\n",
    "# Compute the p-value\n",
    "p = np.sum(bs_diff_replicates >= mean_diff) / len(bs_diff_replicates)\n",
    "\n",
    "# Print p-value\n",
    "print('p =', p)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "** EDA of beak length and depth **"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Make scatter plot of 1975 data\n",
    "_ = plt.plot(bl_1975, bd_1975, marker='.',\n",
    "             linestyle='none', color='blue', alpha=0.5)\n",
    "\n",
    "# Make scatter plot of 2012 data\n",
    "_ = plt.plot(bl_2012, bd_2012, marker='.',\n",
    "             linestyle='none', color='red', alpha=0.5)\n",
    "\n",
    "# Label axes and make legend\n",
    "_ = plt.xlabel('beak length (mm)')\n",
    "_ = plt.ylabel('beak depth (mm)')\n",
    "_ = plt.legend(('1975', '2012'), loc='upper left')\n",
    "\n",
    "# Show the plot\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Linear regressions**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Compute the linear regressions\n",
    "slope_1975, intercept_1975 = np.polyfit(bl_1975, bd_1975, 1)\n",
    "slope_2012, intercept_2012 = np.polyfit(bl_2012, bd_2012, 1)\n",
    "\n",
    "# Perform pairs bootstrap for the linear regressions\n",
    "bs_slope_reps_1975, bs_intercept_reps_1975 = \\\n",
    "        draw_bs_pairs_linreg(bl_1975, bd_1975, 1000)\n",
    "bs_slope_reps_2012, bs_intercept_reps_2012 = \\\n",
    "        draw_bs_pairs_linreg(bl_2012, bd_2012, 1000)\n",
    "\n",
    "# Compute confidence intervals of slopes\n",
    "slope_conf_int_1975 = np.percentile(bs_slope_reps_1975, [2.5, 97.5])\n",
    "slope_conf_int_2012 = np.percentile(bs_slope_reps_2012, [2.5, 97.5])\n",
    "intercept_conf_int_1975 = np.percentile(\n",
    "                            bs_intercept_reps_1975, [2.5, 97.5])\n",
    "intercept_conf_int_2012 = np.percentile(\n",
    "                            bs_intercept_reps_2012, [2.5, 97.5])\n",
    "\n",
    "# Print the results\n",
    "print('1975: slope =', slope_1975,\n",
    "      'conf int =', slope_conf_int_1975)\n",
    "print('1975: intercept =', intercept_1975,\n",
    "      'conf int =', intercept_conf_int_1975)\n",
    "print('2012: slope =', slope_2012,\n",
    "      'conf int =', slope_conf_int_2012)\n",
    "print('2012: intercept =', intercept_2012,\n",
    "      'conf int =', intercept_conf_int_2012)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Displaying the linear regression results**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Make scatter plot of 1975 data\n",
    "_ = plt.plot(bl_1975, bd_1975, marker='.',\n",
    "             linestyle='none', color='blue', alpha=0.5)\n",
    "\n",
    "# Make scatter plot of 2012 data\n",
    "_ = plt.plot(bl_2012, bd_2012, marker='.',\n",
    "             linestyle='none', color='red', alpha=0.5)\n",
    "\n",
    "# Label axes and make legend\n",
    "_ = plt.xlabel('beak length (mm)')\n",
    "_ = plt.ylabel('beak depth (mm)')\n",
    "_ = plt.legend(('1975', '2012'), loc='upper left')\n",
    "\n",
    "# Generate x-values for bootstrap lines: x\n",
    "x = np.array([10, 17])\n",
    "\n",
    "# Plot the bootstrap lines\n",
    "for i in range(100):\n",
    "    plt.plot(x, bs_slope_reps_1975[i] * x + bs_intercept_reps_1975[i],\n",
    "             linewidth=0.5, alpha=0.2, color='blue')\n",
    "    plt.plot(x, bs_slope_reps_2012[i] * x + bs_intercept_reps_2012[i],\n",
    "             linewidth=0.5, alpha=0.2, color='red')\n",
    "\n",
    "# Draw the plot again\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Beak length to depth ratio**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Compute length-to-depth ratios\n",
    "ratio_1975 = bl_1975 / bd_1975\n",
    "ratio_2012 = bl_2012 / bd_2012\n",
    "\n",
    "# Compute means\n",
    "mean_ratio_1975 = np.mean(ratio_1975)\n",
    "mean_ratio_2012 = np.mean(ratio_2012)\n",
    "\n",
    "# Generate bootstrap replicates of the means\n",
    "bs_replicates_1975 = draw_bs_reps(ratio_1975, np.mean, size=10000)\n",
    "bs_replicates_2012 = draw_bs_reps(ratio_2012, np.mean, size=10000)\n",
    "\n",
    "# Compute the 99% confidence intervals\n",
    "conf_int_1975 = np.percentile(bs_replicates_1975, [0.5, 99.5])\n",
    "conf_int_2012 = np.percentile(bs_replicates_2012, [0.5, 99.5])\n",
    "\n",
    "# Print the results\n",
    "print('1975: mean ratio =', mean_ratio_1975,\n",
    "      'conf int =', conf_int_1975)\n",
    "print('2012: mean ratio =', mean_ratio_2012,\n",
    "      'conf int =', conf_int_2012)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Correlation of offspring and parental data**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def draw_bs_pairs(x, y, func, size=1):\n",
    "    \"\"\"Perform pairs bootstrap for single statistic.\"\"\"\n",
    "\n",
    "    # Set up array of indices to sample from: inds\n",
    "    inds = np.arange(len(x))\n",
    "\n",
    "    # Initialize replicates: bs_replicates\n",
    "    bs_replicates = np.empty(size)\n",
    "\n",
    "    # Generate replicates\n",
    "    for i in range(size):\n",
    "        bs_inds = np.random.choice(inds, len(inds))\n",
    "        bs_x, bs_y = x[bs_inds], y[bs_inds]\n",
    "        bs_replicates[i] = func(bs_x, bs_y)\n",
    "\n",
    "    return bs_replicates"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Pearson correlation of offspring and parental data**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Compute the Pearson correlation coefficients\n",
    "r_scandens = pearson_r(bd_parent_scandens, bd_offspring_scandens)\n",
    "r_fortis = pearson_r(bd_parent_fortis, bd_offspring_fortis)\n",
    "\n",
    "# Acquire 1000 bootstrap replicates of Pearson r\n",
    "bs_replicates_scandens = draw_bs_pairs(\n",
    "        bd_parent_scandens, bd_offspring_scandens, pearson_r, size=1000)\n",
    "bs_replicates_fortis = draw_bs_pairs(\n",
    "        bd_parent_fortis, bd_offspring_fortis, pearson_r, size=1000)\n",
    "\n",
    "# Compute 95% confidence intervals\n",
    "conf_int_scandens = np.percentile(bs_replicates_scandens, [2.5, 97.5])\n",
    "conf_int_fortis = np.percentile(bs_replicates_fortis, [2.5, 97.5])\n",
    "\n",
    "# Print results\n",
    "print('G. scandens:', r_scandens, conf_int_scandens)\n",
    "print('G. fortis:', r_fortis, conf_int_fortis)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Measuring heritability**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def heritability(parents, offspring):\n",
    "    \"\"\"Compute the heritability from parent and offspring samples.\"\"\"\n",
    "    covariance_matrix = np.cov(parents, offspring)\n",
    "    return covariance_matrix[0,1] / covariance_matrix[0,0]\n",
    "\n",
    "# Compute the heritability\n",
    "heritability_scandens = heritability(bd_parent_scandens,\n",
    "                                     bd_offspring_scandens)\n",
    "heritability_fortis = heritability(bd_parent_fortis,\n",
    "                                   bd_offspring_fortis)\n",
    "\n",
    "# Acquire 1000 bootstrap replicates of heritability\n",
    "replicates_scandens = draw_bs_pairs(\n",
    "        bd_parent_scandens, bd_offspring_scandens, heritability, size=1000)\n",
    "replicates_fortis = draw_bs_pairs(\n",
    "        bd_parent_fortis, bd_offspring_fortis, heritability, size=1000)\n",
    "\n",
    "# Compute 95% confidence intervals\n",
    "conf_int_scandens = np.percentile(replicates_scandens, [2.5, 97.5])\n",
    "conf_int_fortis = np.percentile(replicates_fortis, [2.5, 97.5])\n",
    "\n",
    "# Print results\n",
    "print('G. scandens:', heritability_scandens, conf_int_scandens)\n",
    "print('G. fortis:', heritability_fortis, conf_int_fortis)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Is beak depth heritable at all in *G. scandens*?**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Initialize array of replicates: perm_replicates\n",
    "perm_replicates = np.empty(10000)\n",
    "\n",
    "# Draw replicates\n",
    "for i in range(10000):\n",
    "    # Permute parent beak depths\n",
    "    bd_parent_permuted = np.random.permutation(bd_parent_scandens)\n",
    "    perm_replicates[i] = heritability(bd_parent_permuted,\n",
    "                                      bd_offspring_scandens)\n",
    "\n",
    "# Compute p-value: p\n",
    "p = np.sum(perm_replicates >= heritability_scandens) / len(perm_replicates)\n",
    "\n",
    "# Print the p-value\n",
    "print('p-val =', p)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}