Notebook: add example notebook on FRET histogram fitting

notebooks/Example - FRET histogram fitting.ipynb

@@ -0,0 +1,319 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Example - FRET histogram fitting\n",
+    "\n",
+    "*This notebook is part of smFRET burst analysis software [FRETBursts](http://tritemio.github.io/FRETBursts/).*"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "> In this notebook shows how to fit a FRET histogram.\n",
+    "> For a complete tutorial on burst analysis see \n",
+    "> [FRETBursts - us-ALEX smFRET burst analysis](FRETBursts - us-ALEX smFRET burst analysis.ipynb)."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from fretbursts import *"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "sns = init_notebook(apionly=True)\n",
+    "import lmfit\n",
+    "print('lmfit version:', lmfit.__version__)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Tweak here matplotlib style\n",
+    "import matplotlib as mpl\n",
+    "mpl.rcParams['font.sans-serif'].insert(0, 'Arial')\n",
+    "mpl.rcParams['font.size'] = 12\n",
+    "%config InlineBackend.figure_format = 'retina'"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Get and process data"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "url = 'http://files.figshare.com/2182601/0023uLRpitc_NTP_20dT_0.5GndCl.hdf5'\n",
+    "download_file(url, save_dir='./data')\n",
+    "full_fname = \"./data/0023uLRpitc_NTP_20dT_0.5GndCl.hdf5\"\n",
+    "\n",
+    "d = loader.photon_hdf5(full_fname)\n",
+    "loader.alex_apply_period(d)\n",
+    "d.calc_bg(bg.exp_fit, time_s=1000, tail_min_us=(800, 4000, 1500, 1000, 3000))\n",
+    "d.burst_search(L=10, m=10, F=6)\n",
+    "ds = d.select_bursts(select_bursts.size, add_naa=True, th1=30)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Fitting the FRET histogram\n",
+    "\n"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "We start defining the model. Here we choose a 3-Gaussian model:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "model = mfit.factory_three_gaussians()\n",
+    "model.print_param_hints()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "The previsou cell prints all the model parameters. \n",
+    "Each parameters has an initial value and bounds (min, max).\n",
+    "The column `vary` tells if a parameter is varied during the fit\n",
+    "(if False the parameter is fixed).\n",
+    "Parameters with an expression (`Expr` column) are not free but\n",
+    "the are computed as a function of other parameters.\n",
+    "\n",
+    "We can modify the paramenters constrains as follows:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "model.set_param_hint('p1_center', value=0.1, min=-0.1, max=0.3)\n",
+    "model.set_param_hint('p2_center', value=0.4, min=0.4, max=0.6)\n",
+    "model.set_param_hint('p2_sigma', value=0.04, min=0.02, max=0.18)\n",
+    "model.set_param_hint('p3_center', value=0.85, min=0.6, max=1.1)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Then, we fit and plot the model:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "E_fitter = bext.bursts_fitter(ds, 'E', binwidth=0.03)\n",
+    "E_fitter.fit_histogram(model=model, pdf=False, method='nelder')\n",
+    "E_fitter.fit_histogram(model=model, pdf=False, method='leastsq')\n",
+    "dplot(ds, hist_fret, show_model=True, pdf=False);"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "The results are in `E_fitter`:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "res = E_fitter.fit_res[0]\n",
+    "res.params.pretty_print()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "To get a dictionary of values:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "res.values"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "print(res.fit_report())"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# ci = res.conf_interval()\n",
+    "# lmfit.report_ci(ci)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Tidy fit results\n",
+    "\n",
+    "A dataframe of fitted parameters is already in `E_fitter`:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "E_fitter.params"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "With [`pybroom`](http://pybroom.readthedocs.io/) we can get a \"tidy\" DataFrame\n",
+    "with more complete fit results:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import pybroom as br"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "df = br.tidy(res)\n",
+    "df"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Now, we can easily select parameters by name:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "df.loc[df.name.str.contains('center')]"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "df.loc[df.name.str.contains('sigma')]"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python [default]",
+   "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.5.2"
+  },
+  "toc": {
+   "colors": {
+    "hover_highlight": "#DAA520",
+    "running_highlight": "#FF0000",
+    "selected_highlight": "#FFD700"
+   },
+   "moveMenuLeft": true,
+   "nav_menu": {
+    "height": "264px",
+    "width": "252px"
+   },
+   "navigate_menu": true,
+   "number_sections": false,
+   "sideBar": true,
+   "threshold": 4,
+   "toc_cell": false,
+   "toc_position": {
+    "height": "673px",
+    "left": "0px",
+    "right": "1139.11px",
+    "top": "107px",
+    "width": "212px"
+   },
+   "toc_section_display": "block",
+   "toc_window_display": true,
+   "widenNotebook": false
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 1
+}