{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "## Physics 503, Lecture 19: Nonlinear Regression\n",
    "### Phys 503, J.R. Gladden, University of Mississippi"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "### Fitting data to a nonlinear model\n",
    "For models which can **not** be linearized, we must take a different approach.  This is always an *iterative* process. A long standing and popular algorithm is the Levenberg-Marquardt scheme. Conceptual details are presented in the slides.  We will not develop \"from scratch\" code to perform this optiomization, but will take advantage of the 'curve_fit' function supplied by scipy.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "### Example: Damped harmonic oscillator\n",
    "The model to fit against the data is: $$y(t) = A_0 e^{-t/\\tau} sin(\\omega t + \\phi) $$\n",
    "so there are 4 free parameters to be determined - these are the 'degrees of freedom'.\n",
    "\n",
    "** Imports and function definitions **"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "%matplotlib wx\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.optimize import curve_fit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "\n",
    "def avg(x):\n",
    "    return sum(x)/float(len(x))\n",
    "\n",
    "\n",
    "def sst(x):\n",
    "    sum=0.0\n",
    "    xbar=avg(x)\n",
    "    for i in range(len(x)):\n",
    "        sum += (x[i]-xbar)**2\n",
    "    return sum\n",
    "\n",
    "def ssr(x,y,model,params):\n",
    "    sum = 0.0\n",
    "    for i in range(len(y)):\n",
    "        sum += (y[i] - model(x[i],*params))**2\n",
    "    return sum\n",
    "\n",
    "def fitstats(fit):\n",
    "\t#Fit statistics\n",
    "\tfitparams=fit[0]\n",
    "\tcov=fit[1]                      #Covariance matrix\n",
    "\tdof=len(xdata)-len(fitparams)   #degrees of freedom\n",
    "\tchisqr=ssr(xdata,ydata,funct,fitparams)     #chi-square is sum of squares of diagonal covariance\n",
    "\terrbars=20\n",
    "\t\n",
    "\t#Computer uncertainties in estimated parameters from covariance matrix and reduced chisqr\n",
    "\tdel_params=[]\n",
    "\tfor i in range(len(fitparams)):\n",
    "\t\tdel_params.append(np.sqrt(cov[i,i]*np.sqrt(chisqr/dof)))\n",
    "\t\n",
    "\treturn chisqr,del_params"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "### Usage"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Fit Results (error)\n",
      "++++++++++++++++++++++++++++++\n",
      "A= 2.136 +/- 0.063 cm\n",
      "omega= 2.979 +/- 0.012 Hz\n",
      "tau= 3.915 +/- 0.179 sec\n",
      "phi= 0.468 +/- 0.031 radians\n",
      "++++++++++++++++++++++++++++++\n",
      "Correlation Coefficient (R^2): 0.916\n"
     ]
    }
   ],
   "source": [
    "def funct(t,a,w,tau,phi):\n",
    "    return a*np.exp(-t/tau)*np.sin(w*t+phi)\n",
    "plt.close('all')\n",
    "\n",
    "#\n",
    "errs=0.2\n",
    "#Generated simulated data\n",
    "par=[2.,3.0,4.0,0.5]\n",
    "\n",
    "xdata=np.arange(0,10,0.2)\n",
    "ydata=funct(xdata,*par)+errs*np.random.randn(len(xdata))\n",
    "\n",
    "#initial guess for parameters\n",
    "par0=(10.5,10.0,1.0,0.1)\n",
    "\n",
    "#Perform fit\n",
    "fit=curve_fit(funct,xdata,ydata,p0=par0)\n",
    "\n",
    "#Parse output and extract some statistics\n",
    "fitparams=fit[0]\n",
    "a_fit=fitparams[0]\n",
    "omega_fit=fitparams[1]\n",
    "tau_fit=fitparams[2]\n",
    "phi_fit=fitparams[3]\n",
    "\n",
    "\n",
    "chisqr,del_params=fitstats(fit)\n",
    "\n",
    "#Generate fit curve\n",
    "xfit=np.linspace(min(xdata),max(xdata),100)\n",
    "yfit=funct(xfit,a_fit,omega_fit,tau_fit,phi_fit)\n",
    "\n",
    "\n",
    "#Plot it all\n",
    "plt.plot(xdata,ydata,'o',label='Data')\n",
    "plt.plot(xfit,yfit,'-',label='Fit')\n",
    "plt.xlabel('Time (s)')\n",
    "plt.ylabel('Amplitude (cm)')\n",
    "plt.legend()\n",
    "plt.show()\n",
    "\n",
    "print 'Fit Results (error)'\n",
    "print '+'*30\n",
    "print 'A= %2.3f +/- %1.3f cm' % (a_fit, del_params[0]) \n",
    "print 'omega= %2.3f +/- %1.3f Hz' % (omega_fit, del_params[1])\n",
    "print 'tau= %2.3f +/- %1.3f sec' % (tau_fit, del_params[2])\n",
    "print 'phi= %2.3f +/- %1.3f radians' % (phi_fit, del_params[3])\n",
    "\n",
    "r=np.sqrt(1.0-ssr(xdata,ydata,funct,fitparams)/sst(ydata))\n",
    "rsq=r**2\n",
    "\n",
    "print '+'*30\n",
    "print 'Correlation Coefficient (R^2): %2.3f' % rsq\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "### Exercise: Gaussian Fit\n",
    "The model for a Gaussian function, commonly used to describe distributions, is $$ y(x) = A \\exp\\left[-\\frac{(x-x_0)^2}{2 \\sigma^2} \\right]$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "def gauss(x,a,x0,sigma): return a*np.exp(-(x-x0)**2/(2.0*sigma**2))\n",
    "    \n",
    "xdata=np.linspace(0,10,100)\n",
    "errs=0.2\n",
    "ydata=gauss(xdata,5.0,3.5,0.6)+errs*np.random.randn(len(xdata))\n",
    "\n",
    "#Initial Guess\n",
    "par0=(2.0,5.0,1.0)\n",
    "\n",
    "fit=curve_fit(gauss,xdata,ydata,p0=par0)\n",
    "\n",
    "fitparams=fit[0]\n",
    "A_fit,x0_fit,sigma_fit=fitparams\n",
    "xfit=np.linspace(min(xdata),max(xdata),100)\n",
    "yfit=gauss(xfit,A_fit,x0_fit,sigma_fit)\n",
    "\n",
    "plt.figure(2)\n",
    "plt.plot(xdata,ydata,'bo',label='Data')\n",
    "plt.plot(xfit,yfit,'g-',label='Gaussian Fit',lw=2)\n",
    "plt.xlabel('x')\n",
    "plt.ylabel('y')\n",
    "plt.legend()\n",
    "plt.show()\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[  6.66066746e-03,   4.91943626e-11,  -5.51747415e-04],\n",
       "       [  4.91943626e-11,   1.37114737e-04,  -5.98494217e-12],\n",
       "       [ -5.51747415e-04,  -5.98494217e-12,   1.37114734e-04]])"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fit[1]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}