Parameter Estimation using SciPy

gistfile1.json

{
 "metadata": {
  "name": ""
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import numpy as np\n",
      "import matplotlib.pyplot as plt\n",
      "from scipy.integrate import odeint\n",
      "from scipy.optimize import leastsq\n",
      "from matplotlib import cm\n",
      "from itertools import combinations\n",
      "from collections import OrderedDict"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 22
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def michelis_menten(y, t, *args):\n",
      "\tVmax = args[0]\n",
      "\tkm = args[1]\n",
      "\tSt = args[2]\n",
      "\tP = y[0]\n",
      "\tS = St - P\n",
      "\n",
      "\tdP = Vmax * (S / (S+km))\n",
      "\treturn dP\n",
      "\n",
      "def hill_michelis_menten(y, t, *args):\n",
      "\tVmax = args[0]\n",
      "\tkm = args[1]\n",
      "\tSt = args[2]\n",
      "\tH = np.max((0, args[3]))\n",
      "\tP = y[0]\n",
      "\tS = St - P\n",
      "\n",
      "\tdP = Vmax * np.power(S,H) / (np.power(S, H)+km)\n",
      "\treturn dP"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 23
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def mm_residuals(params, *args):\n",
      "\n",
      "\tt_exp = args[0]\n",
      "\ty_exp = args[1]\n",
      "\t# Experimental noisy data\n",
      "\n",
      "\tP_0 = 0\n",
      "\t# Initial Condition is 0\n",
      "\n",
      "\tparams = tuple(np.exp(params))\n",
      "\n",
      "\tt_sim = np.linspace(0, 100, 1000)\n",
      "\ty_sim = odeint(michelis_menten, [P_0], t_sim, args = params).flatten()\n",
      "\tmapped_t = np.argmin(np.abs(np.subtract.outer(t_sim, t_exp)), axis = 0)\n",
      "\t# Maps the experimental timepoints to the closest simulated timepoint\n",
      "\n",
      "\treturn (y_sim[mapped_t] - y_exp)\n",
      "\n",
      "def hmm_residuals(params, *args):\n",
      "\n",
      "\tt_exp = args[0]\n",
      "\ty_exp = args[1]\n",
      "\t# Experimental noisy data\n",
      "\n",
      "\tP_0 = 0\n",
      "\t# Initial Condition is 0\n",
      "\n",
      "\tparams = tuple(np.exp(params))\n",
      "\n",
      "\tt_sim = np.linspace(0, 100, 1000)\n",
      "\ty_sim = odeint(hill_michelis_menten, [P_0], t_sim, args = params).flatten()\n",
      "\tmapped_t = np.argmin(np.abs(np.subtract.outer(t_sim, t_exp)), axis = 0)\n",
      "\t# Maps the experimental timepoints to the closest simulated timepoint\n",
      "\n",
      "\treturn (y_sim[mapped_t] - y_exp)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 24
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Parameters MM\n",
      "Vmax = 1\n",
      "km = 3\n",
      "St = 10\n",
      "mm_params = (Vmax, km, St)\n",
      "\n",
      "# Initial Conditions MM\n",
      "P_0 = 0\n",
      "\n",
      "# Timesteps\n",
      "n_steps = 1000\n",
      "t = np.linspace(0, 100, n_steps)\n",
      "\n",
      "# Regular Solution\n",
      "num_P = odeint(michelis_menten, [P_0], t, args = (mm_params))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 25
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Generating noisy experimental data\n",
      "noisy_P = num_P + np.random.randn(n_steps, 1)*0.5\n",
      "noisy_P[noisy_P < 0] = 0\n",
      "noisy_P = noisy_P[::20].flatten()\n",
      "noisy_t = t[::20]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 26
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Michelis-Menten Least Squares\n",
      "mm_init_guess = (0.2, 0.3, 15)\n",
      "mm_initial_P = odeint(michelis_menten, [P_0], t, args = tuple(mm_init_guess)).flatten()\n",
      "log_init_guess = np.log(mm_init_guess) # We log transform the parameters to avoid negative values\n",
      "mm_initial_rss = np.sum(mm_residuals(log_init_guess, *(noisy_t, noisy_P))**2) # Residual Sum of Squares\n",
      "\n",
      "out = leastsq(mm_residuals, log_init_guess, args = (noisy_t, noisy_P), full_output = 1)\n",
      "mm_log_params = out[0].flatten()\n",
      "mm_lsq_params = np.exp(mm_log_params)\n",
      "mm_lsq_P = odeint(michelis_menten, [P_0], t, args = tuple(mm_lsq_params)).flatten()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 27
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Scaling MM Covariance Matrix:\n",
      "mm_cov_P = out[1]\n",
      "n = len(noisy_P)\n",
      "m = len(mm_init_guess) # Number of Parameters\n",
      "rss = np.sum(mm_residuals(mm_log_params, *(noisy_t, noisy_P))**2) # Residual Sum of Squares\n",
      "\n",
      "# Note - this covariance matrix of f(exp(x)) - which is:\n",
      "# e^2x*f''(x) + e^x*f'(x)\n",
      "# At minima, f'(x) ~ 0\n",
      "mm_cov_P = mm_cov_P / (np.exp(2*mm_log_params))\n",
      "\n",
      "# Rescale by residuals\n",
      "mm_cov_P = mm_cov_P * rss / (n - m)\n",
      "mm_aic = n * np.log(rss/n) + 2 * m"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 28
    },
    {
     "cell_type": "heading",
     "level": 1,
     "metadata": {},
     "source": [
      "Trying out more generic model now:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Hill MM Least Squares\n",
      "hmm_init_guess = (0.2, 0.3, 15, 1)\n",
      "log_init_guess = np.log(hmm_init_guess)\n",
      "out = leastsq(hmm_residuals, log_init_guess, args = (noisy_t, noisy_P), \n",
      "\tfull_output = 1, maxfev = 10000, xtol = 0.0001)\n",
      "hmm_log_params = out[0].flatten()\n",
      "hmm_lsq_params = np.exp(hmm_log_params)\n",
      "\n",
      "# Hill MM ODE solution with lsq params\n",
      "hmm_lsq_P = odeint(hill_michelis_menten, [P_0], t, args = tuple(hmm_lsq_params)).flatten()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 29
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Scaling Hill MM Covariance Matrix:\n",
      "hmm_cov_P = out[1]\n",
      "\n",
      "n = len(noisy_P)\n",
      "m = len(hmm_init_guess) # Number of Parameters\n",
      "rss = np.sum(hmm_residuals(hmm_log_params, *(noisy_t, noisy_P))**2) # Residual Sum of Squares\n",
      "\n",
      "# Note - this covariance matrix of f(exp(x)) - which is:\n",
      "# e^2x*f''(x) + e^x*f'(x)\n",
      "# At minima, f'(x) ~ 0\n",
      "hmm_cov_P = hmm_cov_P / (np.exp(2*hmm_log_params))\n",
      "\n",
      "# Rescale by residuals\n",
      "hmm_cov_P = hmm_cov_P * rss / (n - m)\n",
      "hmm_aic = n * np.log(rss/n) + 2 * m"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 30
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Drawing vertical lines between sim with initial solutions and exp points\n",
      "mapped_t = np.argmin(np.abs(np.subtract.outer(t, noisy_t)), axis = 0)\n",
      "paired_measures = np.vstack((mm_initial_P[mapped_t], noisy_P))\n",
      "paired_measures.sort(axis = 0)\n",
      "y_min = paired_measures[0, :]\n",
      "y_max = paired_measures[1, :]\n",
      "\n",
      "\n",
      "plt.figure()\n",
      "plt.plot(noisy_t, noisy_P, 'bo')\n",
      "plt.plot(t, mm_initial_P, 'g-')\n",
      "plt.vlines(noisy_t, y_min, y_max, 'k')\n",
      "plt.xlabel('Time')\n",
      "plt.ylabel('P')\n",
      "plt.legend(['Noisy Data', 'Initial Guess'], loc = 'best')\n",
      "plt.title('Michelis-Menten Initial Solution')\n",
      "plt.text(25, 2, ('RSS = %.3f' %mm_initial_rss), fontsize = 18)\n",
      "plt.savefig('Michelis-Menten_Initial_Solution.jpg', dpi = 500,\n",
      "\tbbox_inches='tight')\n",
      "\n",
      "paired_measures = np.vstack((mm_lsq_P[mapped_t], noisy_P))\n",
      "paired_measures.sort(axis = 0)\n",
      "y_min = paired_measures[0, :]\n",
      "y_max = paired_measures[1, :]\n",
      "\n",
      "plt.figure()\n",
      "plt.plot(t, num_P, 'r-')\n",
      "plt.plot(noisy_t, noisy_P, 'bo')\n",
      "plt.plot(t, mm_lsq_P, 'g-')\n",
      "plt.vlines(noisy_t, y_min, y_max, 'k')\n",
      "plt.xlabel('Time')\n",
      "plt.ylabel('P')\n",
      "plt.legend(['Ideal Parameters', 'Noisy Data', 'Least Square Fit'], loc = 'best')\n",
      "plt.title('Michelis-Menten Least Squares')\n",
      "plt.text(35, 5.5, ('Vmax = %.3f +/- %.3f (%d)'\n",
      "\t\t% (mm_lsq_params[0], np.sqrt(mm_cov_P[0][0]), Vmax)), fontsize = 15)\n",
      "plt.text(35, 4.75, ('Km = %.3f +/- %.3f (%d)' \n",
      "\t\t% (mm_lsq_params[1], np.sqrt(mm_cov_P[1][1]), km)), fontsize = 15)\n",
      "plt.text(35, 4, ('St = %.3f +/- %.3f (%d)' \n",
      "\t\t% (mm_lsq_params[2], np.sqrt(mm_cov_P[2][2]), St)), fontsize = 15)\n",
      "\n",
      "plt.text(25, 2, ('RSS = %.3f' %rss), fontsize = 18)\n",
      "plt.savefig('Michelis-Menten_Least_Squares.jpg', dpi = 500, bbox_inches='tight')\n",
      "\n",
      "plt.text(25, 1, ('AIC = %.3f' %mm_aic), fontsize = 18)\n",
      "\n",
      "plt.savefig('Michelis-Menten_Least_Squares_AIC.jpg', dpi = 500, bbox_inches='tight')\n",
      "\n",
      "# Simple Binding Plots:\n",
      "paired_measures = np.vstack((hmm_lsq_P[mapped_t], noisy_P))\n",
      "paired_measures.sort(axis = 0)\n",
      "y_min = paired_measures[0, :]\n",
      "y_max = paired_measures[1, :]\n",
      "\n",
      "\n",
      "\n",
      "\n",
      "plt.figure()\n",
      "plt.plot(t, num_P, 'r-')\n",
      "plt.plot(noisy_t, noisy_P, 'bo')\n",
      "plt.plot(t, hmm_lsq_P, 'g-')\n",
      "plt.vlines(noisy_t, y_min, y_max, 'k')\n",
      "plt.xlabel('Time')\n",
      "plt.ylabel('P')\n",
      "plt.legend(['Ideal Parameters', 'Noisy Data', 'Least Square Fit'], loc = 'best')\n",
      "plt.title('Hill Function Least Squares')\n",
      "plt.text(35, 5.5, ('Vmax = %.3f +/- %.3f (%d)'\n",
      "\t\t% (mm_lsq_params[0], np.sqrt(hmm_cov_P[0][0]), Vmax)), fontsize = 15)\n",
      "plt.text(35, 4.75, ('Km = %.3f +/- %.3f (%d)' \n",
      "\t\t% (mm_lsq_params[1], np.sqrt(hmm_cov_P[1][1]), km)), fontsize = 15)\n",
      "plt.text(35, 4, ('St = %.3f +/- %.3f  (%d)'\n",
      "\t\t% (hmm_lsq_params[2], np.sqrt(hmm_cov_P[2][2]), St)), fontsize = 15)\n",
      "plt.text(35, 3.25, ('H = %.3f +/- %.3f (%d)' \n",
      "\t\t% (hmm_lsq_params[3], np.sqrt(hmm_cov_P[3][3]), 1)), fontsize = 15)\n",
      "plt.text(25, 2, ('RSS = %.3f' %rss), fontsize = 18)\n",
      "plt.savefig('Hill_Function_Least_Squares.jpg', dpi = 500, bbox_inches='tight')\n",
      "\n",
      "plt.text(25, 1, ('AIC = %.3f' %hmm_aic), fontsize = 18)\n",
      "plt.savefig('Hill_Function_Least_Squares_AIC.jpg', dpi = 500, bbox_inches='tight')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 31
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Parameters MM\n",
      "Vmax = 1.0\n",
      "km = 3.0\n",
      "St = 10.0\n",
      "\n",
      "mm_param_ranges = {}\n",
      "mm_param_ranges['Vmax'] = np.logspace(np.log10(Vmax-0.5), np.log10(Vmax+3), 400)\n",
      "mm_param_ranges['km']= np.logspace(np.log10(km-2), np.log10(km+2), 400)\n",
      "mm_param_ranges['St']= np.logspace(np.log10(St-2), np.log10(St+2), 400)\n",
      "\n",
      "for param_1_name, param_2_name in combinations(mm_param_ranges.keys(), 2):\n",
      "\n",
      "    mm_params_dict = OrderedDict()\n",
      "    mm_params_dict['Vmax'] = 1\n",
      "    mm_params_dict['km'] = 3\n",
      "    mm_params_dict['St'] = 10\n",
      "\n",
      "    param_1_range = mm_param_ranges[param_1_name]\n",
      "    param_2_range = mm_param_ranges[param_2_name]\n",
      "\n",
      "    rss = np.zeros((400, 400))\n",
      "    for i, param_1 in enumerate(param_1_range):\n",
      "        for j, param_2 in enumerate(param_2_range):\n",
      "            mm_params_dict[param_1_name] = param_1\n",
      "            mm_params_dict[param_2_name] = param_2\n",
      "            mm_params = tuple(mm_params_dict.values())\n",
      "            mm_log_params = np.log(mm_params)\n",
      "            rss[i, j] = np.sum(mm_residuals(mm_log_params, *(noisy_t, noisy_P))**2)\n",
      "\n",
      "    rss = np.log(rss)\n",
      "\n",
      "    param_1_range, param_2_range = np.meshgrid(param_1_range, param_2_range)\n",
      "    fig = plt.figure()\n",
      "    ax = fig.add_subplot(111,aspect='equal')\n",
      "    im = ax.pcolor(param_1_range, param_2_range, rss, cmap = plt.cm.RdBu)\n",
      "\n",
      "    height = np.min(param_2_range) - np.max(param_2_range)\n",
      "    length = np.min(param_1_range) - np.max(param_1_range)\n",
      "    ax.set_aspect(float(length)/height)\n",
      "\n",
      "    CS = ax.contour(param_1_range, param_2_range, rss,\n",
      "        cmap = cm.spectral)\n",
      "\n",
      "    plt.clabel(CS, inline=1, fmt='%1.2f', fontsize=14)\n",
      "    plt.xlabel(param_1_name)\n",
      "    plt.ylabel(param_2_name)\n",
      "    plt.title('log RSS')\n",
      "\n",
      "    plt.title('Fitness Space of Parameters')\n",
      "    plt.savefig('%s_%s_parameter_space.jpg' % (param_1_name, param_2_name),\n",
      "        dpi = 500,  bbox_inches='tight')\n",
      "plt.show()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 32
    }
   ],
   "metadata": {}
  }
 ]
}