usALEX - Corrections - Gamma factor fit - extended.ipynb

{
 "metadata": {
  "name": "",
  "signature": "sha256:fe0a31dc4cb03e999787044bd1e24d18b3a9e42b6f5bdea398d5a058fad72aef"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "# Fit Gamma factor\n",
      "\n",
      "> *This notebook estimates the gamma factor from a set of 5 \u03bcs-ALEX smFRET measurements.*\n",
      "\n",
      "## What this notebook does?\n",
      "\n",
      "According to [Lee 2005](http://dx.doi.org/10.1529/biophysj.104.054114) ([PDF](http://www.chem.ucla.edu/~michalet/papers/BJ2005.pdf), [SI PDF](http://www.chem.ucla.edu/~michalet/papers/BJ2005SI.pdf)), we estimate the $\\gamma$-factor \n",
      "from Proximity Ratio (PR) and S values (with background, leakage and direct excitation correction) \n",
      "for a set of 5 \u03bcs-ALEX measurements.\n",
      "\n",
      "The PR and S values are computed by the template notebook:\n",
      "\n",
      "- [usALEX-5samples-PR-leakage-dir-ex-all-ph](usALEX-5samples-PR-leakage-dir-ex-all-ph.ipynb)\n",
      "\n",
      "which is executed for each sample using [usALEX-Batch](usALEX-Batch.ipynb).\n",
      "\n",
      "From [Lee 2005](http://dx.doi.org/10.1529/biophysj.104.054114) (equation 20), the following linear relation holds:\n",
      "\n",
      "$$\\frac{1}{S} = \\Omega + \\Sigma \\cdot E_{PR}$$\n",
      "\n",
      "Once $\\Omega$ and $\\Sigma$ are fitted, we can compute the $\\gamma$-factor as ([equation 22](http://www.sciencedirect.com/science/article/pii/S0006349505733464#eq22)):\n",
      "\n",
      "$$\\gamma = (\\Omega-1)/(\\Omega + \\Sigma-1)$$"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Import libraries"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from __future__ import division\n",
      "import numpy as np\n",
      "import pandas as pd\n",
      "import lmfit\n",
      "from scipy.stats import linregress"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 3
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Computation"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "This notebook read data from the file:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "data_file = 'results/usALEX-5samples-PR-leakage-dir-ex-all-ph.txt'"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 6
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "data = pd.read_csv(data_file, sep=\"\\s+\").set_index('sample')\n",
      "data"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "html": [
        "<div style=\"max-height:1000px;max-width:1500px;overflow:auto;\">\n",
        "<table border=\"1\" class=\"dataframe\">\n",
        "  <thead>\n",
        "    <tr style=\"text-align: right;\">\n",
        "      <th></th>\n",
        "      <th>n_bursts_all</th>\n",
        "      <th>n_bursts_do</th>\n",
        "      <th>n_bursts_fret</th>\n",
        "      <th>E_pr_fret</th>\n",
        "      <th>E_pr_fret_sig</th>\n",
        "      <th>S_pr_fret</th>\n",
        "      <th>S_pr_fret_sig</th>\n",
        "      <th>E_pr_do_kde</th>\n",
        "      <th>nt_mean</th>\n",
        "      <th>E_pr_fret_kde</th>\n",
        "    </tr>\n",
        "    <tr>\n",
        "      <th>sample</th>\n",
        "      <th></th>\n",
        "      <th></th>\n",
        "      <th></th>\n",
        "      <th></th>\n",
        "      <th></th>\n",
        "      <th></th>\n",
        "      <th></th>\n",
        "      <th></th>\n",
        "      <th></th>\n",
        "      <th></th>\n",
        "    </tr>\n",
        "  </thead>\n",
        "  <tbody>\n",
        "    <tr>\n",
        "      <th>7d</th>\n",
        "      <td> 1169</td>\n",
        "      <td> 582</td>\n",
        "      <td>  544</td>\n",
        "      <td> 0.926413</td>\n",
        "      <td> 0.002498</td>\n",
        "      <td> 0.553918</td>\n",
        "      <td> 0.004235</td>\n",
        "      <td> 0.0014</td>\n",
        "      <td> 22.265796</td>\n",
        "      <td> 0.9302</td>\n",
        "    </tr>\n",
        "    <tr>\n",
        "      <th>12d</th>\n",
        "      <td> 1304</td>\n",
        "      <td> 329</td>\n",
        "      <td>  945</td>\n",
        "      <td> 0.731053</td>\n",
        "      <td> 0.003017</td>\n",
        "      <td> 0.565995</td>\n",
        "      <td> 0.003207</td>\n",
        "      <td> 0.0164</td>\n",
        "      <td> 21.905752</td>\n",
        "      <td> 0.7432</td>\n",
        "    </tr>\n",
        "    <tr>\n",
        "      <th>17d</th>\n",
        "      <td> 2483</td>\n",
        "      <td> 462</td>\n",
        "      <td> 1958</td>\n",
        "      <td> 0.427139</td>\n",
        "      <td> 0.002573</td>\n",
        "      <td> 0.543997</td>\n",
        "      <td> 0.002362</td>\n",
        "      <td> 0.0126</td>\n",
        "      <td> 21.094356</td>\n",
        "      <td> 0.4318</td>\n",
        "    </tr>\n",
        "    <tr>\n",
        "      <th>22d</th>\n",
        "      <td> 2051</td>\n",
        "      <td> 319</td>\n",
        "      <td> 1670</td>\n",
        "      <td> 0.183379</td>\n",
        "      <td> 0.001969</td>\n",
        "      <td> 0.549327</td>\n",
        "      <td> 0.002511</td>\n",
        "      <td> 0.0008</td>\n",
        "      <td> 22.987313</td>\n",
        "      <td> 0.1796</td>\n",
        "    </tr>\n",
        "    <tr>\n",
        "      <th>27d</th>\n",
        "      <td>  787</td>\n",
        "      <td> 160</td>\n",
        "      <td>  586</td>\n",
        "      <td> 0.083907</td>\n",
        "      <td> 0.002890</td>\n",
        "      <td> 0.562717</td>\n",
        "      <td> 0.004169</td>\n",
        "      <td>-0.0068</td>\n",
        "      <td> 16.941794</td>\n",
        "      <td> 0.0844</td>\n",
        "    </tr>\n",
        "  </tbody>\n",
        "</table>\n",
        "</div>"
       ],
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 7,
       "text": [
        "        n_bursts_all  n_bursts_do  n_bursts_fret  E_pr_fret  E_pr_fret_sig  \\\n",
        "sample                                                                       \n",
        "7d              1169          582            544   0.926413       0.002498   \n",
        "12d             1304          329            945   0.731053       0.003017   \n",
        "17d             2483          462           1958   0.427139       0.002573   \n",
        "22d             2051          319           1670   0.183379       0.001969   \n",
        "27d              787          160            586   0.083907       0.002890   \n",
        "\n",
        "        S_pr_fret  S_pr_fret_sig  E_pr_do_kde    nt_mean  E_pr_fret_kde  \n",
        "sample                                                                   \n",
        "7d       0.553918       0.004235       0.0014  22.265796         0.9302  \n",
        "12d      0.565995       0.003207       0.0164  21.905752         0.7432  \n",
        "17d      0.543997       0.002362       0.0126  21.094356         0.4318  \n",
        "22d      0.549327       0.002511       0.0008  22.987313         0.1796  \n",
        "27d      0.562717       0.004169      -0.0068  16.941794         0.0844  "
       ]
      }
     ],
     "prompt_number": 7
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "data[['E_pr_fret', 'E_pr_fret_kde']]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "html": [
        "<div style=\"max-height:1000px;max-width:1500px;overflow:auto;\">\n",
        "<table border=\"1\" class=\"dataframe\">\n",
        "  <thead>\n",
        "    <tr style=\"text-align: right;\">\n",
        "      <th></th>\n",
        "      <th>E_pr_fret</th>\n",
        "      <th>E_pr_fret_kde</th>\n",
        "    </tr>\n",
        "    <tr>\n",
        "      <th>sample</th>\n",
        "      <th></th>\n",
        "      <th></th>\n",
        "    </tr>\n",
        "  </thead>\n",
        "  <tbody>\n",
        "    <tr>\n",
        "      <th>7d</th>\n",
        "      <td> 0.926413</td>\n",
        "      <td> 0.9302</td>\n",
        "    </tr>\n",
        "    <tr>\n",
        "      <th>12d</th>\n",
        "      <td> 0.731053</td>\n",
        "      <td> 0.7432</td>\n",
        "    </tr>\n",
        "    <tr>\n",
        "      <th>17d</th>\n",
        "      <td> 0.427139</td>\n",
        "      <td> 0.4318</td>\n",
        "    </tr>\n",
        "    <tr>\n",
        "      <th>22d</th>\n",
        "      <td> 0.183379</td>\n",
        "      <td> 0.1796</td>\n",
        "    </tr>\n",
        "    <tr>\n",
        "      <th>27d</th>\n",
        "      <td> 0.083907</td>\n",
        "      <td> 0.0844</td>\n",
        "    </tr>\n",
        "  </tbody>\n",
        "</table>\n",
        "</div>"
       ],
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 8,
       "text": [
        "        E_pr_fret  E_pr_fret_kde\n",
        "sample                          \n",
        "7d       0.926413         0.9302\n",
        "12d      0.731053         0.7432\n",
        "17d      0.427139         0.4318\n",
        "22d      0.183379         0.1796\n",
        "27d      0.083907         0.0844"
       ]
      }
     ],
     "prompt_number": 8
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "res = linregress(data.E_pr_fret_kde, 1/data.S_pr_fret)\n",
      "slope, intercept, r_val, p_val, stderr = res"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 9
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "For more info see [`scipy.stats.linearregress`](http://docs.scipy.org/doc/scipy-0.13.0/reference/generated/scipy.stats.linregress.html)."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "Sigma = slope \n",
      "Sigma"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 10,
       "text": [
        "-0.009993289402355424"
       ]
      }
     ],
     "prompt_number": 10
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "Omega = intercept\n",
      "Omega"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 11,
       "text": [
        "1.8063091259676307"
       ]
      }
     ],
     "prompt_number": 11
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "[Pearson correlation coefficient](http://en.wikipedia.org/wiki/Pearson_product-moment_correlation_coefficient):"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "r_val"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 10,
       "text": [
        "0.030997465223121211"
       ]
      }
     ],
     "prompt_number": 10
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "[Coefficient of determination](http://en.wikipedia.org/wiki/Coefficient_of_determination) $R^2$:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "r_val**2"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 11,
       "text": [
        "0.00096084285025860891"
       ]
      }
     ],
     "prompt_number": 11
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "P-value (to test the null hypothesis that the slope is zero):"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "p_val"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 12,
       "text": [
        "0.96053912269852137"
       ]
      }
     ],
     "prompt_number": 12
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Gamma computed from the previous fitted values:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "gamma = (Omega - 1)/(Omega + Sigma - 1)\n",
      "gamma"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 12,
       "text": [
        "1.012549404323615"
       ]
      }
     ],
     "prompt_number": 12
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "with open('results/usALEX - gamma factor - all-ph.txt', 'w') as f:\n",
      "    f.write(str(gamma))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 13
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Corrected FRET values"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "With the knoledge of $\\gamma$ we can compute the corrected FRET efficiencies.\n",
      "\n",
      "Let define the function to apply gamma correction (from `burstlib.py`):"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def gamma_correct_E(Er, gamma):\n",
      "    \"\"\"Apply gamma correction to the uncorrected FRET `Er`.\"\"\"\n",
      "    Er = np.asarray(Er)\n",
      "    return Er/(gamma-gamma*Er+Er)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 14
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The corrected values are:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "E_alex = gamma_correct_E(data.E_pr_fret, gamma)\n",
      "E_alex"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 15,
       "text": [
        "array([ 0.92456139,  0.73109077,  0.42336969,  0.18263484,  0.07914672])"
       ]
      }
     ],
     "prompt_number": 15
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "data.E_pr_fret"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 16,
       "text": [
        "sample\n",
        "7d        0.924344\n",
        "12d       0.730478\n",
        "17d       0.422609\n",
        "22d       0.182170\n",
        "27d       0.078920\n",
        "Name: E_pr_fret, dtype: float64"
       ]
      }
     ],
     "prompt_number": 16
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "# Non-linear least squares, using [Lmfit](http://cars9.uchicago.edu/software/python/lmfit/)\n",
      "\n",
      "With a non-linear least-squares minimization we are not limited to fit a linear relation.\n",
      "\n",
      "Let's rewrite the relation without inverting $S_R$:\n",
      "\n",
      "$$S_R = \\frac{1}{\\Omega + \\Sigma \\cdot E_{PR}}$$"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import lmfit\n",
      "from lmfit import minimize, Parameters, report_fit, report_ci, report_errors"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 17
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def residual(params, Er, Sr, eps_data=None):\n",
      "    slope = params['slope'].value\n",
      "    intercept = params['intercept'].value\n",
      "\n",
      "    model = 1./(intercept + slope*Er)\n",
      "    if eps_data is None:\n",
      "        eps_data = np.ones_like(Er)\n",
      "    return (Sr - model)/eps_data\n",
      "\n",
      "params = Parameters()\n",
      "params.add('intercept', value=2)\n",
      "params.add('slope', value=0.)\n",
      "\n",
      "out = minimize(residual, params, args=(data.E_pr_fret, data.S_pr_fret))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 18
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "out.values"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 19,
       "text": [
        "{'intercept': 1.7406596630381004, 'slope': 0.0030796355446684408}"
       ]
      }
     ],
     "prompt_number": 19
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "report_fit(params)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "[[Variables]]\n",
        "     intercept:     1.74066 +/- 0.02472581 (1.42%) initial =  2\n",
        "     slope:         0.003079636 +/- 0.0436766 (1418.24%) initial =  0\n",
        "[[Correlations]] (unreported correlations are <  0.100)\n"
       ]
      }
     ],
     "prompt_number": 20
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "Omega = out.values['intercept']\n",
      "Sigma = out.values['slope']\n",
      "\n",
      "gamma = (Omega - 1)/(Omega + Sigma - 1)\n",
      "gamma"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 21,
       "text": [
        "0.99585925397442776"
       ]
      }
     ],
     "prompt_number": 21
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We can compute the confidence intervals ([method details](http://cars9.uchicago.edu/software/python/lmfit/confidence.html)) and derive the effect on $\\gamma$:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "ci = lmfit.conf_interval(out)\n",
      "report_ci(ci)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stderr",
       "text": [
        "/Users/anto/anaconda/lib/python2.7/site-packages/lmfit/confidence.py:289: UserWarning: Warning, maxiter=200 reachedand prob(slope=-0.197920364455) = 0.980791305555 < max(sigmas).\n",
        "  warn(errmsg)\n",
        "/Users/anto/anaconda/lib/python2.7/site-packages/lmfit/confidence.py:289: UserWarning: Warning, maxiter=200 reachedand prob(slope=0.204079635545) = 0.98002107331 < max(sigmas).\n",
        "  warn(errmsg)\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "             99.70%    95.00%    67.40%     0.00%    67.40%    95.00%    99.70%\n",
        "    slope      -inf  -0.13579  -0.04813   0.03940   0.05448   0.14337       inf\n",
        "intercept   1.53797   1.66423   1.71197   1.74066   1.77002   1.82203   1.98211\n"
       ]
      }
     ],
     "prompt_number": 22
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Let's compute gamma in the 95% confidence interval for the 2 parameters using a simple brute-force method:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "gamma_list = []\n",
      "for sig_err in np.linspace(ci['slope'][1][1], ci['slope'][5][1], 20):\n",
      "    for omega_err in  np.linspace(ci['intercept'][1][1], ci['intercept'][5][1], 20):\n",
      "        Omega = out.values['intercept'] + omega_err\n",
      "        Sigma = out.values['slope'] + sig_err\n",
      "\n",
      "        gamma = (Omega - 1)/(Omega + Sigma - 1)\n",
      "        gamma_list.append(gamma)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 23
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "np.min(gamma_list), np.max(gamma_list)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 24,
       "text": [
        "(0.94259831075324385, 1.0584055044039404)"
       ]
      }
     ],
     "prompt_number": 24
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "> **NOTE:** Using the 95% confidence interval for the 2 fit parameters we observe a variation on gamma ~ 12%. This only considering the empirical error of the 5 data points respect to the model. Moreover, the underlying assumption is that the errors are equally distributed for all the data points. So a more realistic estimation could be worst than that."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Just for test we can fit the same linear relation we fitted initially with `linregress`:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def residual2(params, x, data, eps_data=None):\n",
      "    slope = params['slope'].value\n",
      "    intercept = params['intercept'].value\n",
      "\n",
      "    model = (intercept + slope*x)\n",
      "    if eps_data is None:\n",
      "        eps_data = np.ones_like(data)\n",
      "    return (1./data-model)/eps_data\n",
      "\n",
      "params2 = Parameters()\n",
      "params2.add('intercept', value=2)\n",
      "params2.add('slope', value=0.)\n",
      "\n",
      "out2 = minimize(residual2, params2, args=(data.E_pr_fret, data.S_pr_fret))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 25
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "ci = lmfit.conf_interval(out2)\n",
      "report_ci(ci)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stderr",
       "text": [
        "/Users/anto/anaconda/lib/python2.7/site-packages/lmfit/confidence.py:289: UserWarning: Warning, maxiter=200 reachedand prob(slope=-0.198576396493) = 0.9808242408 < max(sigmas).\n",
        "  warn(errmsg)\n",
        "/Users/anto/anaconda/lib/python2.7/site-packages/lmfit/confidence.py:289: UserWarning: Warning, maxiter=200 reachedand prob(slope=0.203423603507) = 0.980824277424 < max(sigmas).\n",
        "  warn(errmsg)\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "             99.70%    95.00%    67.40%     0.00%    67.40%    95.00%    99.70%\n",
        "    slope      -inf  -0.13627  -0.04863   0.03839   0.05348   0.14112       inf\n",
        "intercept   1.52166   1.66272   1.71237   1.74130   1.77024   1.81989   1.96096\n"
       ]
      }
     ],
     "prompt_number": 26
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "Omega = out2.values['intercept']\n",
      "Sigma = out2.values['slope']\n",
      "\n",
      "gamma = (Omega - 1)/(Omega + Sigma - 1)\n",
      "gamma"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 27,
       "text": [
        "0.9967412671491358"
       ]
      }
     ],
     "prompt_number": 27
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Alternative method\n",
      "\n",
      "By definition the stoichiometry is independent from E:\n",
      "\n",
      "$$ S = \\frac{\\gamma F_D + F_A}{\\gamma F_D + F_A + F_{AA}}$$\n",
      "\n",
      "Neglecting $\\gamma$ we obtain a raw stoichiometry (analogous of the proxhimity ratio):\n",
      "\n",
      "$$ S_R  = \\frac{F_D + F_A}{F_D + F_A + F_{AA}}$$\n",
      "\n",
      "Changes of $S_R$ as a function of $E$ indicates $\\gamma \\ne 1$.\n",
      "\n",
      "We can fit $\\gamma$ by minimizing the variance of $S$. However, experimentally, we have $E_R$ and $S_R$ and not $F_D$, $F_A$ and $F_{AA}$. We can express $S_R$ as a function of $E_R$:\n",
      "\n",
      "$$ S_R  = \\frac{1}{1 + \\frac{F_{AA}}{F_D + F_A}} = \\frac{1}{1 + E_R\\frac{F_{AA}}{F_A}} = \\frac{1}{1 + E_R k_2}$$\n",
      "\n",
      "where $k_2 = F_{AA}/F_A$. $k_2$ as a function of $E_R$ and $S_R$ is:\n",
      "\n",
      "$$ k_2 = \\left(\\frac{1}{S_R} - 1\\right)/E_R$$\n",
      "\n",
      "Using this last definition, we can write $S$ as:\n",
      "\n",
      "$$S = \\frac{1}{1 + E k_2} = \\frac{1}{1 + g(E_R, \\gamma) k_2} = f(\\gamma, E_R, S_R)$$\n",
      "\n",
      "where $g(\\cdot)$ function that links $E$ and $E_R$. To find $\\gamma$ we can minimize the variance of $S$ (i.e. $f(\\cdot)$)."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Let start defining $g(\\cdot)$ and $k_2(\\cdot)$:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def correct_E(Er, gamma):\n",
      "    Er = np.asanyarray(Er)\n",
      "    return Er/(gamma-gamma*Er+Er)\n",
      "\n",
      "def k2(Er, Sr):\n",
      "    return (1./Sr - 1)/Er"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 28
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "A direct brute-force search gives:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "gamma_range = np.arange(0.9, 1.1, 0.00002)\n",
      "y_std = np.zeros_like(gamma_range)\n",
      "for i, gamma in enumerate(gamma_range):\n",
      "    y = 1/(1 + correct_E(data.E_pr_fret, gamma)*k2(data.E_pr_fret, data.S_pr_fret))\n",
      "    y_std[i] = y.std()\n",
      "    #plot(y)\n",
      "gamma_range[y_std == y_std.min()]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 29,
       "text": [
        "array([ 0.99696])"
       ]
      }
     ],
     "prompt_number": 29
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Expressing the fit as a minimization problem gives exactly the same result:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def residual(params, Er, Sr):\n",
      "    gamma = params['gamma'].value\n",
      "    S = 1/(1+correct_E(Er, gamma)*k2(Er, Sr))\n",
      "    errors = S - S.mean()\n",
      "    return errors\n",
      "\n",
      "params = Parameters()\n",
      "params.add('gamma', value=1.0)\n",
      "\n",
      "out = minimize(residual, params, args=(data.E_pr_fret, data.S_pr_fret))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 30
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "out.values['gamma']"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 31,
       "text": [
        "0.99696873118951224"
       ]
      }
     ],
     "prompt_number": 31
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "report_fit(params)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "[[Variables]]\n",
        "     gamma:     0.9969687 +/- 0.05088592 (5.10%) initial =  1\n",
        "[[Correlations]] (unreported correlations are <  0.100)\n"
       ]
      }
     ],
     "prompt_number": 32
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "> Here we obtained that the standard deviation of the fitted $\\gamma$ is around 5%, consistent with the previous estimations."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "#Fit plot"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%matplotlib inline\n",
      "from matplotlib.pyplot import *"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 33
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "plot(data.E_pr_fret, 1/data.S_pr_fret, 's')\n",
      "x = np.arange(0, 1, 0.01)\n",
      "plot(x, intercept + slope*x, 'k')\n",
      "ylim(1, 2)\n",
      "xlabel('E (uncorrected)')\n",
      "ylabel('1/S (uncorrected)')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 34,
       "text": [
        "<matplotlib.text.Text at 0x108af1b50>"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEPCAYAAACtCNj2AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAFERJREFUeJzt3X+UXGV9x/H3JoTfQdStwSbAaqBCbQIYlSgog/FgQBRa\nqwZQD/6oHlsRSmul+Ee29VCkRz1qOQfCr6hVEX+gkopQBIcfEShI2GxAMEZTg3CCRBC0YAhM/3ju\nPDM7zMzeuzv3zo99v87Zs/feeWbmu3eT+exzn3ufC5IkSZIkSZIkSZIkSZIkSVLP2hf4EXAPsAH4\naIt2XwA2AmPAYcWUJkkq2j7AocnynsD9wMENbY4Drk6WDwduK6Y0SVK3fRdY1rDtQuCddev3AfMK\nq0iSNMGsgt5nhHBo6PaG7fOBLXXrDwALCqpJktSgiFDYE/gWcDrw+yaPDzWsV3KvSJLU1E45v/4c\n4NvAVwiHjxr9mjAgXbUg2TbBwoULK5s2bcqlQEkaYJuAA7I8Ic+ewhBwKXAv8LkWba4C3pMsLwUe\nA7Y2Ntq0aROVSsWvSoWVK1d2vYZe+XJfuC/cF+2/gIVZP7jz7CkcAbwLWA+sS7adDeyXLK8inHl0\nHPBz4A/Ae3OsR5I0iTxD4RbS9UQ+kmMNkqQMijr7SB1SKpW6XULPcF/UuC9q3BfT03jmT6+qJMfH\nJEkpDQ0NQcbPeXsKkqTIUJAkRYaCJCkyFCRJkaEgSYoMBUlSZChIkiJDQZIUGQqSpMhQkCRFhoIk\nKTIUJEmRoSBJigwFSVJkKEiSIkNBkhQZCpKkyFCQJEWGgiQpMhQkSZGhIEmKDAVJUmQoSJIiQ0GS\nFBkKkqTIUJAkRYaCJCkyFCRJkaEgSYoMBUlSZChIkiJDQZIUGQqSpMhQkCRFeYfCZcBWYLzF48PA\nNcDdwAbg1JzrkSS1kXcorAaWt3n8I8A64FCgBHwG2CnnmiRJLeQdCjcDj7Z5/CFgr2R5L2AbsCPn\nmiRJLXT7r/KLgRuAB4G5wDu6W44kzWzdDoWzCeMJJWAhcB1wCPBEY8PR0dG4XCqVKJVKRdQnSX2j\nXC5TLpen9RpDnSmlrRFgDbCoyWNXA+cAa5P164GPA3c2tKtUKpW86pOkgTQ0NAQZP+e7fUrqfcAb\nk+V5wMuAX3SvHEma2fLuKVwOHEU49XQrsBKYkzy2Ktm+GtiPEFDnAl9r8jozoqdw6qmjbN783O0j\nI/DFL44WXI2kfjeVnkLeYwonTfL4I8Bbcq6hb2zeDDfeONrkkWbbJKnzun34SJLUQwwFSVLU7VNS\npY5xTEaaPkNBA8MxGfWyfvmjxVDoISMj0OwDLGyX1M/65Y8WQ6GH9NJfC5JmJgeaJUlR3/QUxsfD\nLRnqL2Jrttz4vdfa1a9P9px27drVMVPbPfTQXcClQP3zKjz00E+46KKLWr5G43Lj687kdq2e163f\nd+NzerFd489T/f7Tn44DP2OiCvfeO86KFfd15Pd45JFHcuaZZzIdfRMKJ510UvXqvPi91XKvt6tf\nn+w57dq1q2Mmtnv22S3ss8/51Udju0oF7rjjjsJ/P73crtW+7Yd/Z43P6cV2jT/P0NAQd931JA8/\n/GYaLzAeHt7OiSeemPn1mrVbsGAB09U3obBhw4ZulyBJU3bddRuZN2/jc7aPjCxixYoVXaiouSJm\nSe2EGTH3kTSo+uV0zEHTi3MfSVLfnI4pzz6SJNWZ8T0Fu7WSVDPjQ8FurSTVePhIkhTN+J6CpPw5\nr1f/MBQk5c7xuf7h4SNJUjTjewp2ayWpJu2VbnsA+xJmGnsA+ENuFTXnFc2SlFGnr2ieC/wNsAIY\nBrYmLz4P2AZ8FbgY+P0UapUk9aB2YwrfBZ4A3gK8FHgNsBR4CXA8obfwvbwLlCQVxwnxJGlAdfrw\n0Suqr8vEu5ZU3ZXljSRJva9dgpQJYbAbsARYn2xfDNxJOJxUFHsKkpTRVHoK7cYUSsDRwIOEXsOS\n5OuwZJskacCkuXjtIGC8bn0DcHA+5UiSuinNxWvrgUuArxC6IScDY3kWJUnqjjTHmnYDPgy8Llm/\nCbgAeCqvoppwTEGSMprKmELaxrsD+wH3ZaypUwwFScqo0wPNVW8F1gHXJOuHAVdlqkyS1BfShMIo\ncDjwaLK+jnCFsyRpwKQJhaeBxxq2PZtDLZKkLksTCvcApxDOVDoQ+A/gx3kWJUnqjjShcBrwcuCP\nwOXA48AZKZ53GWFm1fE2bUqEw1EbCFdQS5K6KM2o9NuBb6bY1uh1hGm1vwwsavL43sBa4E2EezQM\nA4+0eC3PPpKkjPI6++jslNsa3UxtcLqZk4FvEwIBWgeCJKkg7a5oPhY4DpgPfIFa2swlDD5P14HA\nHOBHyWt+HvjPDryuJGmK2oXCg8BPgBOS79UptJ8A/r4D7z2HMNHeMsLFcbcCtwEbO/DakqQpaBcK\nY8nXlYS7rD2TbJ8N7NKB995COGT0ZPJ1E3AILUJhdHQ0LpdKJUqlUgdKkKTBUS6XKZfL03qNNAMQ\ntwFvpHYv5rnAtcBrUzx3BFhD84Hmg4DzCQPNuwC3A+8E7m3S1oFmScqo03deq9qVWiBAOHy0e4rn\nXQ4cRTiraAuwknDICGAVYR6lawizsD4LXEzzQJAkFSRNgqwFPkoYVwB4JeECNu+8Jkk9LK+ewhmE\naxKqd1t7MeEwjyRpwKRNkDnAy5L299GZU1KzsKcgSRnldfHaHsBZwOmEKStGgOMz1iZJ6gNpQmE1\nsJ3a2UYPAufkVpEkqWvShMJC4DxCMEC4ZkGSNIDShMIfCfdprlqYbJMkDZg0Zx+NEq4nWAB8DTgC\nODW/kiRJ3TJZKMwCng+8DViabDsd+E2eRUmSuiPNqUo/AZbkXcgkPCVVkjKayimpaRp/ijBx3RVM\nHGT+bZY3miZDQZIyyisUNhOmzK5XAV6a5Y2myVCQpIzyCIVZhFtvXjHFmjrFUJCkjPLqKTimIEl9\nyDEFSVLkmIIkKcorFHqBoSBJGeV1P4WdgQ8Dryf0EG4ELqT46bMlSTlLkyCXEsLjS0n7dwM7gA/k\nWFcjewqSlFFeh4/WA4tTbMuToSBJGeV1k50dwAF16wuTbZKkAZNmTOFjwA3AL5P1EeC9eRUkSeqe\ntN2KXQn3aK4APwOeyq2i5jx8JEkZ5XX46COEm+yMEcYSdgP+NmtxkqTelyZBxoBDGrbdDRza+XJa\nsqcgSRnl1VOY1dBuNjAny5tIkvpDmoHma4GvA6sIifMhwu05JUkDJk23YjbwQWBZsn4dcAnwTF5F\nNeHhI0nKyLmPJElRXnMfHQmsJFyfUG1f9CypkqQCpEmQ+4EzgLuYeMjokVwqas6egiRllFdP4THg\nB1MpSJLUX9LeeW02cCXwx7rtd+VSUXP2FCQpo7wGmss8985rAEdneaNpMhQkKSPPPpIkRXmNKawk\n9BSGmNhj+NcsbyRJ6n1pQuEP1MJgN+B44N7cKpIkdc1UDh/tAvw3cNQk7S4D3gw8DCxq0+5VwK3A\nOwiD2c14+EiSMsprQrxGewDzU7RbDSyfpM1s4DzCXEr9Mr4hSQMrzeGj8brlWcCLSDeecDPhKuh2\nTgO+RegtSJK6LE0ovKVueQewFXi6A+89HzgBeAMhFDw+JEld1i4U5gJPAJtTtJmKzwFnUTuzycNH\nktRl7ULhO4R5j74H3An8Ntn+QuCVwInAgcAbp/jeSwj3aQAYBo4l9ECuatZ4dHQ0LpdKJUql0hTf\nVpIGU7lcplwuT+s1Jvvr/A3AycARwJ8m2x4EbgG+SrjauZ0RYA3tzz6CMCi9Bs8+kqSOyePitRuS\nr6m4nHDa6jCwhXARXPU2nqum+JqSpBz1y3F8ewqSlFFR1ylIkgaUoSBJitqNKexBOBtoe7J+EHAc\n4RTVVgPCkqQ+1q6ncA2wf7J8AGF+opcAf0e48Y4kacC0G4AYp3Yq6SeBFxACYWfCXdf+It/SJnCg\nWZIy6vRAc/2n8DLgh8nyduDZTJVJkvpCuzGFceDThIvVFhKmywZ4Ps5TJEkDqV1P4QPANsK4wjGE\nm+0AHEwIC0nSgGl3rOki4AfA9cDjxZTTkmMKkpTRVMYU2jVeSrhJzjLCqanXEs5IGptifdNhKEhS\nRp0OhXrDhENIy4HFwDpCL+IbWd5sGgwFScooz1BofM4S4E3AOVN4/lQYCpKUUVFzH1UI1y8UFQiS\npIJMdZbULcC+nSxkEvYUJCmjTt9PYbzNY/OyvIkkqT+0C4UXEQaWH23y2I/zKUeS1E3tQuH7wJ6E\nM40a3ZhPOZKkbvLOa5I0oLzzmiRpWgwFSVJkKEiSIkNBkhQZCpKkyFCQJEWGgiQpMhQkSZGhIEmK\nDAVJUmQoSJIiQ0GSFBkKkqTIUJAkRYaCJCkyFCRJkaEgSYoMBUlSlHcoXAZsBcZbPH4KMAasB9YC\ni3OuR5LURt6hsBpY3ubxXwCvJ4TBJ4GLcq5HktRG3qFwM/Bom8dvBX6XLN8OLMi5HklSG700pvB+\n4OpuFyFJM9lO3S4gcTTwPuCIVg1GR0fjcqlUolQq5V6UJPWTcrlMuVye1msMdaaUtkaANcCiFo8v\nBq4kjD38vEWbSqVS6XxlkjTAhoaGIOPnfLcPH+1HCIR30ToQJEkFybuncDlwFDBMODV1JTAneWwV\ncAnwl8Cvkm1PA69u8jr2FCQpo6n0FIo4fNQJhoIkZdSPh48kST3EUJAkRYaCJCkyFCRJkaEgSYoM\nBUlSZChIkiJDQZIUGQqSpMhQkCRFhoIkKTIUJEmRoSBJigwFSVJkKEiSIkNBkhQZCpKkyFCQJEWG\ngiQpMhQkSZGhIEmKDAVJUmQoSJIiQ0GSFBkKkqTIUJAkRYaCJCkyFCRJkaEgSYoMBUlSZChIkiJD\nQZIUGQqSpMhQkCRFhoIkKco7FC4DtgLjbdp8AdgIjAGH5VyPJKmNvENhNbC8zePHAQcABwIfBC7I\nuZ6+Vy6Xu11Cz3Bf1LgvatwX05N3KNwMPNrm8bcCX0qWbwf2BublXFNf8x98jfuixn1R476Ynm6P\nKcwHttStPwAs6FItkjTjdTsUAIYa1itdqUKS9JwP5DyMAGuARU0euxAoA19P1u8DjiIMTtf7ObAw\nn/IkaWBtIozb9pQRWp99dBxwdbK8FLitiIIkSd1xOfAgsJ0wdvA+4EPJV9X5hJ7AGPCKoguUJEmS\n1OOWE8YUNgIfb9FmplzoNtm+OIWwD9YDa4HFxZVWuDT/LgBeBewA/qqIorokzb4oAeuADYTxukE1\n2b4YBq4B7ibsi1MLq6xYA3uB8GzCIaQRYA7hF3lwQ5v68YfDGdzxhzT74jXA85Ll5czsfVFtdwPw\nX8DbiiquYGn2xd7APdRO6x4uqriCpdkXo8C5yfIwsA3YqZjyCvU6wgd9mnHbVJ+bvXBKKsCrCb/k\nzcDThLORTmhoM1MudEuzL24Ffpcs387gXtuRZl8AnAZ8C/hNYZUVL82+OBn4NuF6H4BHiiquYGn2\nxUPAXsnyXoRQ2FFQfUXq+AXCvRIKzS5im5+izSB+GKbZF/XeT+0vgUGT9t/FCdSmSBnU61zS7IsD\ngRcAPwLuBN5dTGmFS7MvLgZeTjjRZQw4vZjSek7mz81e6U6l/Y88Ey50y/IzHU04o+uInGrptjT7\n4nPAWUnbIYq59qYb0uyLOYQz+JYBuxN6lLcRjicPkjT74mzCYaUS4Rqn64BDgCfyK6tnZfrc7JVQ\n+DWwb936vtS6wK3aLEi2DZo0+wLC4PLFhDGFdt3HfpZmXyyhdvHjMHAs4ZDCVblXV6w0+2IL4ZDR\nk8nXTYQPwkELhTT74rXAOcnyJuCXwMsIPaiZpG8/N3ci/OJGgJ2ZfKB5kC90S7Mv9iMcU11aaGXF\nS7Mv6q1mcM8+SrMvDgJ+SBiI3Z0w+PjnxZVYmDT74rPAymR5HiE0XlBQfUUbYUAvED4WuJ/wYffP\nybaZeqHbZPviEsLA2brk63+KLrBAaf5dVA1yKEC6ffGPhDOQxoGPFlpdsSbbF8OE6XXGCPvi5KIL\nLIgXCEuSJEmSJEmSJEmSJEmSJEnS1D1D7dqMdcA/tWj3acL0Bv3mecCHp/C8UeAfkuXPEmbOlKSB\nl2bemrl072K+2ZOsT2aE9vPit7KSWigcyOBN8aEC9cosqVKnnECY6qFqM7XpDV5JmEEUwl/XlyXr\nmwjTb1e9h3D1593Al5NtI4R7Nowlr1+dT+aLwIWE6QP+nXBVdXX9PMJkbD8gzLlzE2H+HQhTL3wn\neY+7CffI+FTSfl3yXICPEUJuLKm56hOEK3pvrntNCPMcjRCmSJakgbaDiYeP3t6kzQVMnOril7QO\nhVsIM4u+kDCR3GzCdMv31z2n+uG6htpU1O8lfKBDCIWrqM1Eubph/XrggGT58GQd4Apq01DMIsz5\nvz8TewrHAKvq2qwhHBpaQrjr3q6EntFG4My6532JMA2ElFmvzJIqpfEkk99OcH/CDVYmUwG+T5hR\ndRvwMLAP8AbgG8Bvk3aPJd+XAicmy18h9Aqqr/NNJk5HXF3fk9AD+GbdYzsn348G3pUsPws8znMn\nbDsm+VqXrO9BODw0F7gSeCr5qg8hCHPhjLT8yaU2DAUNovrDojvq1ndtaLe9bvkZwv+H6n0Zmmm1\n/f9arM8ihEqrIEtz74dzgYsatp3e8NzG1xliMO81ogI4pqBB87+Ev/irNhMOG8HE+zc3+0CuEMYN\n3k7tr/bnJ99/DKxIlk8hjA9M5nHC4au/rnvPxcny9dTONJpNOHz0BKEXUHUtYdbLPZL1+cCfJO99\nIrXDR8czMQReTPi5pcwMBfWT3Zg4pvBvTdrcQi0EAP4F+DxwB6HXUP3wrND8r+l7CTdnuZEwAPyZ\nZPtphLGEMUIo1N/esfF16tdPIdwy9W5gA+GeuSTPP5owNnAn4X4A24C1hHGF8wh3C/sa4Q5q6wmH\ntfZMfvYrklqu5rlnWx2WPEeSZrw9CQEwU/0ZnpKqach6HrXU67YTBmN3MDMPoXyCcBOmX3W7EEmS\nJEmSJEmSJEmSJEmSJEkz0v8De16Xk7bRHNUAAAAASUVORK5CYII=\n",
       "text": [
        "<matplotlib.figure.Figure at 0x108ae7910>"
       ]
      }
     ],
     "prompt_number": 34
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 30
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 30
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 30
    }
   ],
   "metadata": {}
  }
 ]
}