bootstrapping/bootstrapping_fill_in.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Day 3, Morning Session 2, MD/MC workshop\n",
    "### Dr. Michael Shirts, CU Boulder"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# load up some modules.\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import scipy\n",
    "import scipy.stats as stats\n",
    "import pandas as pd\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Bootstrapping for complicated uncertainties.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There is a difference in statistics between a **population** and a **sample**. The population is all the possible observations out there. For instance if I were an epidemiologist, this might be all the people in the U.S. or all the children in Massachusetts.\n",
    "\n",
    "\n",
    "<img width=\"304\" height=\"232\" alt=\"Image result for population versus sample\" src=\"http://sphweb.bumc.bu.edu/otlt/MPH-Modules/BS/BS704_BiostatisticsBasics/Sampling3.jpg\">\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "But this can also be applied to measurements in the lab as well such as possible voltage values from this battery. If you ran an infinite number of experients, you would get an infinite number of measurements, but either in this case of the whole US population case: accessing the actual population is virtually impossible.\n",
    "\n",
    "<img width=\"304\" height=\"223\" alt=\"Image result for thermometer gif\" src=\"https://media.giphy.com/media/26FL3uMhARSAvIZZS/giphy.gif\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As a scientist, you only have access to a sample. Part of designing an experiment is choosing how big your sample should be.\n",
    "\n",
    "But a key problem is: if you change your sample, it could change your sample mean and sample variance significantly. So the question is, how can you understand the variation in the **population** by only looking at your **sample**.\n",
    "\n",
    "Let me repeat: the question, which the bootstrap is trying to answer, is to understand the variation in the  **population** by understanding the  variation at your **sample**.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's say we were in the lab, and we were making a calibration curve, say absorbance against concentration using UV-vis measurements. You generate the following plot:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def generate_data(xlim=[0,5],ylim=[0,100],npoints=20,noise=1.0,seed=None,noise_model=\"even\"):\n",
    "    '''\n",
    "    This function just generates some random data that might look like at UV-vis output.\n",
    "    Inputs: range of data in the x and y direction, the number of points, and the level of noise.\n",
    "    Add random number seed to allow repetability\n",
    "    \n",
    "    We have two types of noise; even noise, and the other exponential noise model.\n",
    "    '''\n",
    "    np.random.seed(seed) # so we can control the random noise\n",
    "    x = np.linspace(xlim[0], xlim[1], npoints)\n",
    "    y_raw = np.linspace(ylim[0], ylim[1], npoints)\n",
    "    if noise_model == \"even\":\n",
    "        y = y_raw + noise*np.random.rand(npoints)\n",
    "    else:\n",
    "        y = y_raw + np.exp(x*np.random.rand(npoints)) - 1\n",
    "    return x, y "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "xmin = 0\n",
    "xmax = 5\n",
    "\n",
    "ymin = 0\n",
    "ymax = 100\n",
    "\n",
    "npoints = 50\n",
    "\n",
    "# in the absence of noise, the slope is 25 and the intercept is 0, \n",
    "# but because of the heteroscedasticity, the linear fit will be somewhat biased.\n",
    "\n",
    "x, y = generate_data(xlim=[xmin,xmax],ylim=[ymin,ymax],noise=1.0,npoints=npoints,noise_model=\"exponenential\", seed=1)\n",
    "\n",
    "results = stats.linregress(x,y)\n",
    "slope = results.slope\n",
    "intercept = results.intercept\n",
    "y_fit = intercept + np.array([xmin,xmax])*slope\n",
    "original_data = pd.DataFrame({'x': x, 'y':y})\n",
    "original_data.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we plot the data and the fit."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(8,6))\n",
    "ax.scatter(original_data['x'], original_data['y'], alpha=0.5, label='raw data')\n",
    "ax.plot(np.array([0,5]), y_fit, linewidth=4, label='linear fit')\n",
    "ax.set_xlim([0, 5])\n",
    "ax.set_ylim([0, 200])\n",
    "ax.legend(loc=\"lower right\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "But uh oh. You can see that some are those points have error getting a little large. You want to get some indication of how your calibration curve would change if you repeated the experiment again. \n",
    "\n",
    "Instead of actually repeating the experiment, you decide to use bootstrapping, and *estimate the variation within your sample as a replacement for the variation within your population.*\n",
    "\n",
    "To do this, you will generate \"new\" samples, picking from the data **with replacement**, pretending the data you collected is **entirely representative** of the actual population) **with replacement** and perform a regression on this resampling of your sample.\n",
    "**With replacement**, means you can draw the same point multiple times, is an important statement. We are assuming our population is much, much, bigger than our sample, but has the exact same distribution as our sample.  So when we withdraw the value at has x=1.456, there are still infinte more points with x=1.456 remaining. So the next point you pick is just as likely to be x=1.456 as it was the first time you drew it.\n",
    "\n",
    "You will generate some code to do this. First, though, we will see what happens if we resample from the original population, then we will see what bootstrapping can do. \n",
    "\n",
    "So first, we are drawing from the _population_, because we are regenerating the entire dataset, including with new amounts of noise."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def draw_from_population(data,npoints=20,ndraws=10):\n",
    "    '''\n",
    "    Draw from the population using the generate_data function,\n",
    "    and determining the slope and intercept for each process. \n",
    "    '''\n",
    "    fig, ax = plt.subplots(figsize=(8,6))\n",
    "    ax.scatter(data['x'], data['y'], alpha=0.5, label='raw data')\n",
    "\n",
    "    xmin = data['x'].min()\n",
    "    xmax = data['x'].max()\n",
    "    \n",
    "    # initializing outputs\n",
    "    slope = np.zeros(ndraws)\n",
    "    intercept = np.zeros(ndraws)\n",
    "\n",
    "    for i in range(ndraws):\n",
    "        # vvvv - these are the lines to change for bootstrap\n",
    "        #x, y = generate_data(xlim=[xmin,xmax],ylim=[ymin,ymax],npoints=npoints,noise=10.0) # don't set the seed, so we get new random data.\n",
    "        x, y = generate_data(xlim=[xmin,xmax],ylim=[ymin,ymax],npoints=npoints,noise=1.0,noise_model=\"exponenential\") # don't set the seed, so we get new random data.\n",
    "        sample = pd.DataFrame({'x': x, 'y':y})\n",
    "        # ^^^^ \n",
    "        results = stats.linregress(sample['x'].values,sample['y'])\n",
    "        slope[i] = results.slope\n",
    "        intercept[i] = results.intercept\n",
    "        y_fit = intercept[i] + np.array([0,5])*slope[i]\n",
    "        ax.plot(np.array([xmin,xmax]), y_fit, linewidth=2, color='b', alpha=0.2)\n",
    "    return slope, intercept"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "draw_from_population(original_data,npoints=20,ndraws=10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can repeat this, like, 5000 times, generate a distribution of slopes and intercepts, and use the distributions to calculate confidence intervals on the slope and intercept:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "slopes, intercepts = draw_from_population(original_data,npoints=50,ndraws=5000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "plt.hist(slopes,bins=30)\n",
    "plt.title('slope values')\n",
    "plt.figure()\n",
    "plt.hist(intercepts,bins=30)\n",
    "plt.title('intercept values')\n",
    "ci = 95\n",
    "\n",
    "print('mean slope: {:.2f} +/- {:.2f} ({:2d}% CI: {:.2f} - {:.2f})'.format(np.mean(slopes), \n",
    "                                                       np.std(slopes,ddof=1), ci,\n",
    "                                                       np.percentile(slopes, ((100-ci)/2)),\n",
    "                                                       np.percentile(slopes, ((100+ci)/2))))\n",
    "\n",
    "print('mean intercept: {:.2f} +/- {:.2f} ({:2d}% CI: {:.2f} - {:.2f})'.format(np.mean(intercepts), \n",
    "                                                       np.std(intercepts,ddof=1), ci,\n",
    "                                                       np.percentile(intercepts, ((100-ci)/2)),\n",
    "                                                       np.percentile(intercepts, ((100+ci)/2))))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's see how these values compare to the ones estimated by the standard linear estimator."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "results = stats.linregress(original_data['x'],original_data['y'])\n",
    "print(f\"slope + standard error of slope = {results.slope:.2f} +/- {results.stderr:.2f}\")\n",
    "print(f\"intercept + standard error of intercept = {results.intercept:.2f} +/- {results.intercept_stderr:.2f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "They seem statistically consistent!  Both the constant term and the x term are within the error bars; or equivalently, the values are well within the distributions. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Now, to the bootstrap\n",
    "\n",
    "To pull a sample with replacement, we will take advantage of `pandas` sampling capabilities:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "original_data.head()\n",
    "subsample = original_data.sample(n=original_data.x.count(), replace=True)\n",
    "print(subsample)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that some samples occur more than once, and some don't appear at all! That is what with replacement means!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Hacking Time**: That code draws *one* bootstrap sample. Copy the `bootstrap_from_population` function from above to create a new function `bootstrap_from_data` to plot the lines from _multiple bootstrap samples_ of the original dataset instead of multiple draws from the population."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that bootstrap means will not be exactly the same as the population means, as each draw from the population will be different, and if we determine the bootstrap means more precisely, they will converge to the value for that sample, which should be well WITHIN the distribution of populations means (it will be have the distribution of the population distribution of means).  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Once you have the code working, then run the line below - you should hopefully get something similar to the draw from population!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "slopes_boot, intercepts_boot = draw_from_bootstrap(original_data,ndraws=5000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Hacking time, Part 2**: Plot the histograms of the distribution of slopes and intercepts of the population and the bootstrap samples and compare (just adapt the code from plotting the population distributions)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'plt' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "Input \u001b[0;32mIn [1]\u001b[0m, in \u001b[0;36m<cell line: 1>\u001b[0;34m()\u001b[0m\n\u001b[0;32m----> 1\u001b[0m \u001b[43mplt\u001b[49m\u001b[38;5;241m.\u001b[39mfigure()\n\u001b[1;32m      2\u001b[0m plt\u001b[38;5;241m.\u001b[39mhist(slopes,bins\u001b[38;5;241m=\u001b[39m\u001b[38;5;241m30\u001b[39m,alpha\u001b[38;5;241m=\u001b[39m\u001b[38;5;241m0.4\u001b[39m,density\u001b[38;5;241m=\u001b[39m\u001b[38;5;28;01mTrue\u001b[39;00m,label\u001b[38;5;241m=\u001b[39m\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mpopulation\u001b[39m\u001b[38;5;124m'\u001b[39m)\n\u001b[1;32m      3\u001b[0m plt\u001b[38;5;241m.\u001b[39mhist(slopes_boot,bins\u001b[38;5;241m=\u001b[39m\u001b[38;5;241m30\u001b[39m,alpha\u001b[38;5;241m=\u001b[39m\u001b[38;5;241m0.4\u001b[39m,color\u001b[38;5;241m=\u001b[39m\u001b[38;5;124m'\u001b[39m\u001b[38;5;124morange\u001b[39m\u001b[38;5;124m'\u001b[39m,density\u001b[38;5;241m=\u001b[39m\u001b[38;5;28;01mTrue\u001b[39;00m,label\u001b[38;5;241m=\u001b[39m\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mbootstrap\u001b[39m\u001b[38;5;124m'\u001b[39m)\n",
      "\u001b[0;31mNameError\u001b[0m: name 'plt' is not defined"
     ]
    }
   ],
   "source": [
    "plt.figure()\n",
    "plt.hist(slopes,bins=30,alpha=0.4,density=True,label='population')\n",
    "plt.hist(slopes_boot,bins=30,alpha=0.4,color='orange',density=True,label='bootstrap')\n",
    "plt.legend()\n",
    "plt.title('slope values')\n",
    "ci = 95\n",
    "print('mean slope (bootstrap): {:.2f} +/- {:.2f} ({:2d}% CI: {:.2f} - {:.2f})'.format(np.mean(slopes_boot), \n",
    "                                                       np.std(slopes_boot,ddof=1),ci,\n",
    "                                                       np.percentile(slopes_boot, ((100-ci)/2)),\n",
    "                                                       np.percentile(slopes_boot, ((100+ci)/2))))\n",
    "plt.show()\n",
    "\n",
    "plt.figure()\n",
    "plt.hist(intercepts,bins=30,alpha=0.4,density=True,label='population')\n",
    "plt.hist(intercepts_boot,bins=30,alpha=0.4,color='orange',density=True,label='bootstrap')\n",
    "plt.title('intercept values')\n",
    "plt.legend()\n",
    "ci = 95\n",
    "print('mean intercept (bootstrap): {:.2f} +/- {:.2f} ({:2d}% CI: {:.2f} - {:.2f})'.format(np.mean(intercepts_boot), \n",
    "                                                       np.std(intercepts_boot,ddof=1),ci,\n",
    "                                                       np.percentile(intercepts_boot, ((100-ci)/2)),\n",
    "                                                       np.percentile(intercepts_boot, ((100+ci)/2))))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Bootstrapping for complicated function distributions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Bootstrapping is generalizable to any statistical quantities you may be interested in. One useful example is error propagation. If you know the error in raw datasets, and you want to know how that error will impact downstream calculations, you can use bootstrapping to do it **without having to do the calculus** involved in standard error propagation, or if the error propagation is really impossible to do. \n",
    "\n",
    "For example, try obtaining the distribution of the log of the absolute value of the (intercept/slope) in the linear fit above. What is the average and the distribution of this quantity? (Why this quantity? No reason, really, except that it would be REALLY HARD to get the distribution of errors using any standard error propagation). Compare this to what one would get using bootstrapping.\n",
    "\n",
    "Similarly, bootstrapping can be used to estimate error in ANY fitting procedure extending beyond linear regression."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "vals = np.sin(np.log(np.fabs(slopes/intercepts)))\n",
    "vals_boot = np.sin(np.log(np.fabs(slopes_boot/intercepts_boot)))\n",
    "plt.hist(vals,bins=30,alpha=0.4, density=True,label='population')\n",
    "plt.hist(vals_boot,bins=30,alpha=0.4,density=True,color='orange',label='bootstrap')\n",
    "plt.title('f(slope/intercepts) values')\n",
    "plt.legend()\n",
    "ci = 95\n",
    "\n",
    "print('mean sin(log(abs(slope/intercept))): {:.2f} ({:2d}% CI: {:.2f} - {:.2f})'.format(np.mean(vals), ci,\n",
    "                                                       np.percentile(vals, ((100-ci)/2)),\n",
    "                                                       np.percentile(vals, ((100+ci)/2))))\n",
    "print('mean sin(log(abs(slope/intercept))): (Bootstrap): {:.2f} ({:2d}% CI: {:.2f} - {:.2f})'.format(np.mean(vals_boot), ci,\n",
    "                                                       np.percentile(vals_boot, ((100-ci)/2)),\n",
    "                                                       np.percentile(vals_boot, ((100+ci)/2))))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Bootstrapping for computing self-diffusion"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "in 3D, one can estimate the self-diffusion coefficient with the following formula:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$\\langle (x(t)-x(t+\\tau))^2\\rangle = 6D\\tau$ \n",
    "\n",
    "This is for 3D, the 6 = 2 times the number of dimensions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The algorithm is then:\n",
    "\n",
    "1. Calculate the mean square displacement of a particle as a function of time\n",
    "2. Average all the square displacements.\n",
    "3. Fit the result function to a line with intercept = 0 to find the slope.\n",
    "4. Divide the slope by 6."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Load in some data. Each entry in the array is a square displacement of a the oxygen atom of a water molecule. Note that the interval between steps is 0.5 ps, so the total trajectory is 1 ns long. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "msd_data = np.load(\"manytraj.npy\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "nparticles, length = np.shape(msd_data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for i in range(nparticles):\n",
    "    plt.plot(msd_data[i,:],'b',alpha=0.02)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# we can average over all particles\n",
    "avemsd = np.mean(msd_data,axis=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(avemsd)\n",
    "plt.xlabel(\"frames\")\n",
    "plt.ylabel(r\"$\\langle MSD \\rangle$\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We want to minimize $\\sum_i (a\\tau_i-y(\\tau_i))^2$.  A linear fit, but with no intercept."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# frames are from 0 to 1000, 0.5 ps each\n",
    "# this is the function we want to minimize to:\n",
    "def func(a):\n",
    "    return a*np.linspace(0,1000,2001)-avemsd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "slope = scipy.optimize.leastsq(func,0.1)[0][0]\n",
    "D = result/6\n",
    "print(D) # this is in nm^2 / ps. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# to get it into cm^2 / s, this is \n",
    "(D / (10**7 * 10**7))*10**12"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise\n",
    "\n",
    "Perform a bootstrap error analysis by constructing bootstrap samples over _particles_, that are then used to compute the average MSD.\n",
    "\n",
    "1. Plot 500 bootstrapped average MSD's.\n",
    "2. Show the 95% confidence interval region of bootstrapped average MSD trajectories.\n",
    "3. Plot the distribution of diffusion constants  Is it Gaussian? If so what is the standard error in the calculation? "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.10.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}