{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\nCalculate Confidence Intervals\n==============================\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import matplotlib.pyplot as plt\nfrom numpy import argsort, exp, linspace, pi, random, sign, sin, unique\nfrom scipy.interpolate import interp1d\n\nfrom lmfit import (Minimizer, Parameters, conf_interval, conf_interval2d,\n                   report_ci, report_fit)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Define the residual function, specify \"true\" parameter values, and generate\na synthetic data set with some noise:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "def residual(pars, x, data=None):\n    argu = (x*pars['decay'])**2\n    shift = pars['shift']\n    if abs(shift) > pi/2:\n        shift = shift - sign(shift)*pi\n    model = pars['amp']*sin(shift + x/pars['period']) * exp(-argu)\n    if data is None:\n        return model\n    return model - data\n\n\np_true = Parameters()\np_true.add('amp', value=14.0)\np_true.add('period', value=5.33)\np_true.add('shift', value=0.123)\np_true.add('decay', value=0.010)\n\nx = linspace(0.0, 250.0, 2500)\nnoise = random.normal(scale=0.7215, size=x.size)\ndata = residual(p_true, x) + noise"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Create fitting parameters and set initial values:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "fit_params = Parameters()\nfit_params.add('amp', value=13.0)\nfit_params.add('period', value=2)\nfit_params.add('shift', value=0.0)\nfit_params.add('decay', value=0.02)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Set-up the minimizer and perform the fit using leastsq algorithm, and show\nthe report:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "mini = Minimizer(residual, fit_params, fcn_args=(x,), fcn_kws={'data': data})\nout = mini.leastsq()\n\nfit = residual(out.params, x)\nreport_fit(out)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Calculate the confidence intervals for parameters and display the results:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "ci, tr = conf_interval(mini, out, trace=True)\n\nreport_ci(ci)\n\nnames = out.params.keys()\ni = 0\ngs = plt.GridSpec(4, 4)\nsx = {}\nsy = {}\nfor fixed in names:\n    j = 0\n    for free in names:\n        if j in sx and i in sy:\n            ax = plt.subplot(gs[i, j], sharex=sx[j], sharey=sy[i])\n        elif i in sy:\n            ax = plt.subplot(gs[i, j], sharey=sy[i])\n            sx[j] = ax\n        elif j in sx:\n            ax = plt.subplot(gs[i, j], sharex=sx[j])\n            sy[i] = ax\n        else:\n            ax = plt.subplot(gs[i, j])\n            sy[i] = ax\n            sx[j] = ax\n        if i < 3:\n            plt.setp(ax.get_xticklabels(), visible=False)\n        else:\n            ax.set_xlabel(free)\n\n        if j > 0:\n            plt.setp(ax.get_yticklabels(), visible=False)\n        else:\n            ax.set_ylabel(fixed)\n\n        res = tr[fixed]\n        prob = res['prob']\n        f = prob < 0.96\n\n        x, y = res[free], res[fixed]\n        ax.scatter(x[f], y[f], c=1-prob[f], s=200*(1-prob[f]+0.5))\n        ax.autoscale(1, 1)\n        j += 1\n    i += 1"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "It is also possible to calculate the confidence regions for two fixed\nparameters using the function ``conf_interval2d``:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "names = list(out.params.keys())\n\nplt.figure()\ncm = plt.cm.coolwarm\nfor i in range(4):\n    for j in range(4):\n        plt.subplot(4, 4, 16-j*4-i)\n        if i != j:\n            x, y, m = conf_interval2d(mini, out, names[i], names[j], 20, 20)\n            plt.contourf(x, y, m, linspace(0, 1, 10), cmap=cm)\n            plt.xlabel(names[i])\n            plt.ylabel(names[j])\n\n            x = tr[names[i]][names[i]]\n            y = tr[names[i]][names[j]]\n            pr = tr[names[i]]['prob']\n            s = argsort(x)\n            plt.scatter(x[s], y[s], c=pr[s], s=30, lw=1, cmap=cm)\n        else:\n            x = tr[names[i]][names[i]]\n            y = tr[names[i]]['prob']\n\n            t, s = unique(x, True)\n            f = interp1d(t, y[s], 'slinear')\n            xn = linspace(x.min(), x.max(), 50)\n            plt.plot(xn, f(xn), 'g', lw=1)\n            plt.xlabel(names[i])\n            plt.ylabel('prob')\n\nplt.show()"
      ]
    }
  ],
  "metadata": {
    "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.7.6"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}