Scientific_Python-2.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "nbpresent": {
     "id": "81e1b977-fa86-44d7-b5f0-a6f1cca78ab3"
    }
   },
   "source": [
    "<img style=\"float: right;\" src=\"http://www2.le.ac.uk/liscb1.jpg\">\n",
    "Scipy is the first place to look for general-purpose scientific functionality.  The Scipy library is enormous and varied, so covering all the features is a course unto itself.  Having said that, one of the most commonly used features in scipy is it's fitting routines, which we will now explore."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbpresent": {
     "id": "92136e88-e29e-4111-8750-c88ebb58bb51"
    }
   },
   "source": [
    "## Curve fitting\n",
    "\n",
    "Scipy comes with a number of fitting routines.  One that can be extremely useful is `curve_fit`, which fits a function to a given set of data using a least-squares minimization.  Here, we'll fit some (made up) data to the Michaelis-Menten equation:  \n",
    "\n",
    "$$\\large V = \\frac{V_{max}[S]}{K_m+[S]}$$  \n",
    "\n",
    "First, we define the Michaelis-Menten equation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "nbpresent": {
     "id": "e7a24fae-02ea-4f7e-855f-96196cf04ad1"
    }
   },
   "outputs": [],
   "source": [
    "def michaelis_menten(s, km, vmax):\n",
    "    return (vmax*s) / (km + s)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbpresent": {
     "id": "384fdd93-16fc-4462-af13-f4161ead77bb"
    }
   },
   "source": [
    "Now we pick some substrate concentrations we (would) do our measurements at:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "nbpresent": {
     "id": "adb5385b-c9ed-41c0-945d-4041ac358d52"
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "nbpresent": {
     "id": "059ccecf-4c19-4e9b-89c0-de7b3d9310ad"
    }
   },
   "outputs": [],
   "source": [
    "substrate_concentrations = np.array([0.01, 0.05, 0.1, 0.5, 1, 2, 3])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbpresent": {
     "id": "af264ed4-3682-48f5-baf0-6a75b865eca9"
    }
   },
   "source": [
    "For convenience when plotting, we'll create a numpy array of 1000 values between 0 and 3 for our x-axis:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "nbpresent": {
     "id": "277a693f-c70b-42e0-b9bf-104c78772129"
    }
   },
   "outputs": [],
   "source": [
    "substrate_concentration_range = np.linspace(0, 3, 1000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbpresent": {
     "id": "9794ff9d-74b9-43f9-ba34-96ed4ca1486d"
    }
   },
   "source": [
    "Ok, lets make sure this looks right:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "nbpresent": {
     "id": "c70edfd7-6502-4b2d-a01f-909c07343d5a"
    }
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "nbpresent": {
     "id": "15fb9b87-fd8b-49be-8434-ceca9504b7d8"
    }
   },
   "outputs": [],
   "source": [
    "km = 0.1\n",
    "vmax = 0.5\n",
    "\n",
    "initial_velocities = michaelis_menten(s=substrate_concentrations, km=km, vmax=vmax)\n",
    "mm_curve = michaelis_menten(substrate_concentration_range, km, vmax)\n",
    "\n",
    "\n",
    "# Plot all the things!\n",
    "plt.scatter(substrate_concentrations, initial_velocities)\n",
    "plt.plot(substrate_concentration_range, mm_curve)\n",
    "\n",
    "plt.ylabel('Velocity')\n",
    "plt.xlabel('initial [S]')\n",
    "\n",
    "# In M-M kinetics, the Km is the substrate concentration where you've reached half-max rate\n",
    "plt.vlines(km, ymin=0, ymax=vmax/2, linestyle='dashed')\n",
    "plt.hlines(vmax/2, xmin=0, xmax=km, linestyle='dashed')\n",
    "\n",
    "plt.xlim(xmin=0, xmax=3)\n",
    "plt.ylim(ymin=0, ymax=0.6)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbpresent": {
     "id": "91659714-0323-4e6b-a690-b8bdb8fa66f3"
    }
   },
   "source": [
    "Ok, so lets now simulate some noisy data (with Normally distributed noise):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "nbpresent": {
     "id": "fe914a33-f391-4466-a63b-c1a75a242002"
    }
   },
   "outputs": [],
   "source": [
    "number_of_concentrations = len(substrate_concentrations)\n",
    "noise = (np.random.normal(loc=1, scale=0.1, size=number_of_concentrations))\n",
    "simulated_data = michaelis_menten(substrate_concentrations, km, vmax,) * noise"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbpresent": {
     "id": "8b149f2f-ffd7-4b91-8510-cc5227c01078"
    }
   },
   "source": [
    "ok, so how does our simulation look?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "nbpresent": {
     "id": "bd13d505-6702-44ea-9be3-619bbccffa2a"
    }
   },
   "outputs": [],
   "source": [
    "plt.plot(substrate_concentration_range, mm_curve, label='ground truth')\n",
    "plt.scatter(substrate_concentrations, simulated_data, label='simulated data', color='green')\n",
    "plt.legend(loc='lower right');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbpresent": {
     "id": "bd5046a4-4a31-4c56-a058-53d7dcb50be2"
    }
   },
   "source": [
    "Now we can try to fit our original equation.  In order to do this, we need to provide a guess for the parameters that the algorithm can start from.  For something as simple as the M-M equation, even quite bad guesses will do."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "nbpresent": {
     "id": "18bc67a5-1745-42e7-963b-ade414db379b"
    }
   },
   "outputs": [],
   "source": [
    "initial_guess = (100, simulated_data[-1])  # Km and Vmax.  Note the Km is a truly horrible guess, given our data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbpresent": {
     "id": "c571e13e-4227-4ffc-84f9-1c95b8eef6a9"
    }
   },
   "source": [
    "In addition, we can provide the fitting algorithm with bounds - regions of allowed values for the parameters.  In this case, it's not necissary, but it's always good to have a sanity check (in this case, both parameters must be positive.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "nbpresent": {
     "id": "d9ccad30-f3be-47bb-accf-91704a1987d9"
    }
   },
   "outputs": [],
   "source": [
    "lower_bounds = (0, 0)  # Km, Vmax\n",
    "upper_bounds = (np.inf, np.inf)  # Km, Vmax"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "nbpresent": {
     "id": "c05ec43a-8a08-4daa-a925-5a5cfac6327f"
    }
   },
   "outputs": [],
   "source": [
    "from scipy.optimize import curve_fit\n",
    "\n",
    "fitted, covariance = curve_fit(f=michaelis_menten,\n",
    "                               xdata=substrate_concentrations,\n",
    "                               ydata=simulated_data,\n",
    "                               p0=initial_guess,\n",
    "                               bounds=(lower_bounds, upper_bounds)\n",
    "                              )\n",
    "print('Km:', fitted[0])\n",
    "print('Vmax:',fitted[1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbpresent": {
     "id": "c5f2063a-188a-434c-a8de-9cdd40798f33"
    }
   },
   "source": [
    "Remember, always look at your data as much as possible!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "nbpresent": {
     "id": "75495aa9-aca6-46ec-bbd0-1621dbd2494b"
    }
   },
   "outputs": [],
   "source": [
    "calculated_curve = michaelis_menten(substrate_concentration_range, fitted[0], fitted[1])\n",
    "\n",
    "plt.plot(substrate_concentration_range, mm_curve, label='ground truth')\n",
    "plt.scatter(substrate_concentrations, simulated_data, label='simulated data', color='green')\n",
    "plt.plot(substrate_concentration_range, calculated_curve, label='fitted')\n",
    "\n",
    "plt.legend(loc='lower right');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 1: Did the fitting actually work?\n",
    "The easiest way to see if the fitting function returned sensible values is to randomize the starting point. If it keeps giving the same numbers from various starting positions, it's more likely you can trust the results.\n",
    "1. Generate a random starting value for $K_m$ from 0.01 to 10, and for $V_{max}$ from 0.01 to 10.  *Hint: use `np.random.uniform()`\n",
    "2. Fit using these starting values and print both the starting value and resulting parameters.\n",
    "3. Repeat for 10 total fits.  How does it look?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbpresent": {
     "id": "9b7729a3-b3ea-4d49-9c46-a34c87695884"
    }
   },
   "source": [
    "### Exercise 2:  Compute the fitting error\n",
    "Compute the error of the fit.  \n",
    "*Hint: Read the docstring of curve_fit using `.?`, and look at what the function returns*\n",
    "\n",
    "__Bonus:__\n",
    "Remake the plot, but shade the error bounds.  \n",
    "1. The maximum positive error is when you subtract the standard deviation from the Km and add it to the Vmax\n",
    "2. Calculate the curve for the maximum and minimum errors as we did above\n",
    "3. Ask google how to fill between curves in matplotlib\n",
    "4. Try setting `alpha=0.2` in the plotting function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "nbpresent": {
     "id": "b4cde2a2-cea2-4c0e-9017-93436d9bb3bd"
    }
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "nbpresent": {
     "id": "2055d2e7-cc65-47c6-9b6d-e6cca15ffc8a"
    }
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true,
    "nbpresent": {
     "id": "bc357582-a3fc-49e2-9024-8cf4f2e94f65"
    }
   },
   "source": [
    "## Uncertainty"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbpresent": {
     "id": "63f91388-d1aa-4ebd-921b-b64ffc6e02ef"
    }
   },
   "source": [
    "In the fit above, we used the covariance matrix to determine the error bounds.  What's not so apparent from this is that this only works for some fitting routines, is only valid if the errors are normally distributed, and is then only valid if the data is homoscedastic (the measurement error is the same for every point - which isn't really true in this case!)  \n",
    "\n",
    "So what are we to do?\n",
    "1.  Work out the analytical error for every fit, assuming the characteristics of the error (hard - lots of algebra)\n",
    "2.  Use resampling to empirically determine the distribution of parameter values (hard - lots of computation)\n",
    "\n",
    "Except #2 isn't hard - we have computers to do the computation for us.  The type of resampling we'll be doing is call bootstrapping.  The mathematical justification for why this works is (way) beyond this course.  If you follow a few rules, you can approximate the actual underlying distribution on almost any parameter, and then extract what you want directly.  Here are the rules:\n",
    "1. Generate a new set of data by sampling, with replacement, from your current data set.\n",
    "2. The new data set must have the same number of observations as the original (some data will get repeated, some left out)\n",
    "3. Fit your curve with the new dataset\n",
    "4. Repeat a few hundred to tens of thousands of times until the distribution stabilizes\n",
    "\n",
    "An example in code should make this clear:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lower_bounds = (0, 0)  # Km, Vmax\n",
    "upper_bounds = (np.inf, np.inf)  # Km, Vmax\n",
    "initial_guess = (1, simulated_data[-1])  # Km and Vmax.\n",
    "\n",
    "bootstrap_kms =[]\n",
    "bootstrap_vmaxs = []\n",
    "\n",
    "iterations = 100\n",
    "\n",
    "for _ in range(iterations):\n",
    "    # Generate new sample indices with replacement \n",
    "    indices = np.arange(len(substrate_concentrations)-1)\n",
    "    resample_indices = np.random.choice(indices, size=len(indices), replace=True)\n",
    "    resampled_substrate_concentrations = substrate_concentrations[resample_indices]\n",
    "    resampled_simulated_data = simulated_data[resample_indices]\n",
    "    \n",
    "    # Fit the curve with the resampled dataset\n",
    "    bootstrap_fitted, bootstrap_covariance = curve_fit(f=michaelis_menten,\n",
    "                                                       xdata=resampled_substrate_concentrations,\n",
    "                                                       ydata=resampled_simulated_data,\n",
    "                                                       p0=initial_guess,\n",
    "                                                       bounds=(lower_bounds, upper_bounds)\n",
    "                                                      )\n",
    "    bootstrap_kms.append(bootstrap_fitted[0])\n",
    "    bootstrap_vmaxs.append(bootstrap_fitted[1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now look at the results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(14, 4))\n",
    "\n",
    "plt.subplot(121)\n",
    "plt.hist(bootstrap_kms, bins=100)\n",
    "# plt.xlim(0,1)\n",
    "plt.title('Km')\n",
    "\n",
    "plt.subplot(122)\n",
    "plt.hist(bootstrap_vmaxs, bins=100)\n",
    "# plt.xlim(0,2)\n",
    "plt.title('Vmax');\n",
    "\n",
    "len(bootstrap_kms)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 3: Confidence interval\n",
    "The 90% confidence interval is the range of values that contain the central 90% of the distribution.  Extract the 90%  CI for the data above.\n",
    "1. Make your life easier now, change the number of iterations above to 10000 and start it running.\n",
    "2. Make sorted listes of bootstrap_kms and bootstrap_vmaxs\n",
    "3. If you have 10000 sorted data points, the central 9000 data points make up your 90% CI - extract the lowest and hightest values from these 9000\n",
    "4. Use `plt.scatter` to plot all 10000 points.  What do you notice?\n",
    "5. Make two plots comparing the fitting error to the 90% CI.  First, plot the $K_m$ CI using the originally fitted $V_{max}$, then plot the $V_{max}$ CI using the originally fitted $K_m$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One thing to note here:  The confidence intervals *are not centered on the prediction*, this means that, given our data, we're more certain about one bound than the other.  Because we've used random noise, your intervals may appear centered, but if you run the notebook again, you may be able to see that they're not."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 4: Experimental replicates\n",
    "Most experiments will have replicate data points, i.e. the researcher measures the initial rate several times at each substrate concentration.\n",
    "1. Create a new list of substrate concentrations called `substrate_concentrations_with_replicates`, using each concentration three times (order doesn't matter)\n",
    "2. Create new simulated data called `simulated_data_with_replicates` from above, adding the same amount of random noise as for `simulated_data`\n",
    "3. Fit the data and calculate the 90% CI with the replicates.  Make plots of the distributions.\n",
    "What do you notice?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbpresent": {
     "id": "d3b60900-9f96-4f2b-b7d1-962102b13410"
    }
   },
   "source": [
    "## Model choice"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbpresent": {
     "id": "e64b70b6-ff93-4a2b-88fa-12da31ddc7c1"
    }
   },
   "source": [
    "One common difficulty in analysing data is choosing the correct model.  The most tempting choice is the model that fits our data best - *this is extremely dangerous.*  Here's the equation for homotropic alosterism (substrate activation or inhibition of the enzyme):  \n",
    "\n",
    "$$\\large V = \\frac{V_{max}\\frac{[S]}{K_{D1}} + \\beta V_{max}\\frac{[S]^2}{\\alpha K_{D1}K_{D12}}}{1 + \\frac{[S]}{K_{D1}} + \\frac{[S]}{K_{D2}} + \\frac{[S]^2}{\\alpha K_{D1} K_{D2}}}$$  \n",
    "\n",
    "Now, **we know that our system doesn't have any substrate inhibition or activation** - we created the data afterall - but lets see how this does..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "nbpresent": {
     "id": "c3a995f5-78ba-49c3-b74c-ac4e15c513d8"
    }
   },
   "outputs": [],
   "source": [
    "def michaelis_menten_ha(s, kd1, kd2, alpha, beta, vmax):\n",
    "    '''\n",
    "    The Michaelis-Menton function including homotropic allosterism. \n",
    "    '''\n",
    "    numerator = vmax * (s/kd1) + (beta * vmax * (s**2)/(alpha * kd1 *kd2))\n",
    "    denominator = 1 + (s/kd1) + (s/kd2) + (s**2)/(alpha * kd1 * kd2)\n",
    "    \n",
    "    return numerator/denominator"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbpresent": {
     "id": "eac752f9-c35b-473f-b3fb-c192f0ea22c6"
    }
   },
   "source": [
    "If we set $\\alpha$ and $\\beta$ to be 1, and set $K_{D1}$ and $K_{D2}$ to the same value, the equation above reduces to the standard Michaelis-Menton equation, so we'll set our starting parameters there."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "nbpresent": {
     "id": "323208b9-f96b-4fb9-9214-7791ac5ecbab"
    }
   },
   "outputs": [],
   "source": [
    "initial_guess = (fitted[0], fitted[0], 1, 1, fitted[1])  # KD1, KD2, alpha, beta and Vmax\n",
    "lower_bounds_ha = (0, 0, 0, 0, 0)\n",
    "upper_bounds_ha = (np.inf, np.inf, np.inf, np.inf, np.inf)\n",
    "\n",
    "fitted_ha, covariance_ha = curve_fit(f=michaelis_menten_ha,\n",
    "                               xdata=substrate_concentrations,\n",
    "                               ydata=simulated_data,\n",
    "                               p0=initial_guess,\n",
    "                               bounds=(lower_bounds_ha, upper_bounds_ha)\n",
    "                              )\n",
    "\n",
    "print('Original Michaelis-Menton')\n",
    "print('Km:', fitted[0])\n",
    "print('Vmax:',fitted[1])\n",
    "\n",
    "print()\n",
    "\n",
    "print('Michaelis-Menton with aditional parameters')\n",
    "print('KD1', fitted_ha[0])\n",
    "print('KD2:',fitted_ha[1])\n",
    "print('alpha:', fitted_ha[2])\n",
    "print('beta:',fitted_ha[3])\n",
    "print('Vmax:', fitted_ha[4])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbpresent": {
     "id": "e049ca4f-824b-4819-89c6-9db72af79472"
    }
   },
   "source": [
    "Ok, so now which one fits our data better?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "nbpresent": {
     "id": "fbb7a234-1852-4b76-8e66-5069135f765e"
    }
   },
   "outputs": [],
   "source": [
    "calculated_curve_ha = michaelis_menten_ha(substrate_concentration_range,\n",
    "                                          *fitted_ha)  # *fitted_ha => fitted_ha[0], fitted_ha[1]\n",
    "\n",
    "plt.figure(figsize=(14,4))\n",
    "\n",
    "plt.subplot(121)\n",
    "plt.scatter(substrate_concentrations, simulated_data, label='simulated data', color='green')\n",
    "plt.plot(substrate_concentration_range, calculated_curve, label='Fitted Michaelis-Menton', color='red')\n",
    "plt.legend(loc='lower right')\n",
    "\n",
    "plt.subplot(122)\n",
    "plt.scatter(substrate_concentrations, simulated_data, label='simulated data', color='green')\n",
    "plt.plot(substrate_concentration_range, calculated_curve_ha, label='Fitted Michaelis-Menton-HA', color='purple')\n",
    "plt.legend(loc='lower right');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbpresent": {
     "id": "1988469e-94f4-4fc8-9246-da6e4cd12ecb"
    }
   },
   "source": [
    "Which one do you think looks better?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So what does this mean?  Because the  `michaelis_menten_ha` function is the `michaelis_menten` function with three extra parameters, it will **always** fit noisy data better.  Because we used the `michaelis_menten` function to generate the data, we know that's the one we should actually be using, though.  If you're faced with a real-world problem, and the function with more parameters will always fit noisy data better, how can you tell if you really should use the function with more parameters or not?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "nbpresent": {
     "id": "dba9a945-fc3c-4f08-84bb-d90e073d3a4f"
    }
   },
   "source": [
    "## Model selection by cross validation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "How do you know which model to choose?  One of the best ways is cross-validation.  The key to undersanding why this is so powerful is that cross-validation essentially tests how much of the random noise you're modelling.  The function that models the least noise and most signal is almost always the correct one. \n",
    "\n",
    "Cross-validation leaves a bit of data out, then fits (trains) the curve without that data, then tests how well the left out data is predicted by the fit.\n",
    "\n",
    "We'll start with leave-one-out cross validation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "initial_guess = (100, simulated_data[-1])  # Km and Vmax.  Note the Km is a truly horrible guess, given our data\n",
    "lower_bounds = (0, 0)  # Km, Vmax\n",
    "upper_bounds = (np.inf, np.inf)  # Km, Vmax\n",
    "predictions = []\n",
    "\n",
    "# Step through each substrate concentration\n",
    "for i in range(len(substrate_concentrations)):\n",
    "    \n",
    "    # Build a mask excluding the test measurement\n",
    "    mask = np.ones(len(substrate_concentrations)).astype('bool')\n",
    "    mask[i] = False\n",
    "    \n",
    "    # Remove the test measurement from our training data\n",
    "    training_substrate_concentrations = substrate_concentrations[mask]\n",
    "    training_simulated_data = simulated_data[mask]\n",
    "    \n",
    "    # Fit the curve\n",
    "    bootstrap_fitted, bootstrap_covariance = curve_fit(f=michaelis_menten,\n",
    "                                                       xdata=training_substrate_concentrations,\n",
    "                                                       ydata=training_simulated_data,\n",
    "                                                       p0=initial_guess,\n",
    "                                                       bounds=(lower_bounds, upper_bounds)\n",
    "                                                      )\n",
    "    predictions.append(michaelis_menten(substrate_concentrations[i], *fitted))\n",
    "\n",
    "abs_errors = abs(predictions - simulated_data)\n",
    "xval_error = sum(abs_errors)\n",
    "print(\"Cross validation error for 'michaelis_menten':\",xval_error)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 5: Model selection by cross validation\n",
    "Run the cross validation on `michaelis_menten_ha`.  Which function gives lower cross-validation error?  \n",
    "\n",
    "*Hint: If you get `Optimal parameters not found` error, add the parameter `max_nfev=10000` to the `curve_fit` call.*\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So which model do we choose?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "xval_ratio = xval_error / xval_ha_error\n",
    "\n",
    "if xval_ratio < 0.95:\n",
    "    print('Michaelis-Menton is the prefered model.')\n",
    "elif xval_ratio > 1.05:\n",
    "    print('Michaelis-Menton with Homotropic Alostery is the prefered model.')\n",
    "else:\n",
    "    print(\"The cross validation scores are too close, we can't descriminiate between the models\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As always, there are lots of choices for how exactly we do cross-validation, but here are the two most import issues:\n",
    "\n",
    "1. If you only have a few more observations than parameters, do leave-one-out validation.\n",
    "2. If you have correlated data (i.e. several blocks of data,) try to keep the correlated data either in the test or training sets  \n",
    "\n",
    "In general, if you have many more data points than parameters, you can do k-fold cross validation.  For example, with 5-fold cross-validation, you would split your data into 5 chunks, using 4 of the chunks for training and one for error determination.  As you have 5 chunks, you do it 5 times (using a different chunk each time,) and sum up the errors.  The model with the lowest error is the one over-fitting the least."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "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.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}