AustinRochford/gist:92b06d174a7f84fded6e

## gistfile1.txt
{"nbformat_minor": 0, "cells": [{"execution_count": 1, "cell_type": "code", "source": "%matplotlib inline", "outputs": [], "metadata": {"collapsed": true, "trusted": true}}, {"source": "---\ntitle: Maximum Likelihood Estimation of Custom Models in Python with StatsModels\ntags: Statistics, Python\n---", "cell_type": "markdown", "metadata": {}}, {"source": "Maximum likelihood estimation is a common method for fitting statistical models.  In Python, it is quite possible to fit maximum likelihood models using just [`scipy.optimize`](http://docs.scipy.org/doc/scipy/reference/optimize.html). Over time, however, I have come to prefer the convenience provided by [`statsmodels`'](http://statsmodels.sourceforge.net/) [`GenericLikelihoodModel`](http://statsmodels.sourceforge.net/devel/dev/generated/statsmodels.base.model.GenericLikelihoodModel.html).  In this post, I will show how easy it is to subclass `GenericLikelihoodModel` and take advantage of much of `statsmodels`' well-developed machinery for maximum likelihood estimation of custom models.", "cell_type": "markdown", "metadata": {}}, {"source": "#### Zero-inflated Poisson models", "cell_type": "markdown", "metadata": {}}, {"source": "The model we use for this demonstration is a [zero-inflated Poisson model](http://en.wikipedia.org/wiki/Zero-inflated_model#Zero-inflated_Poisson).  This is a model for count data that generalizes the Poisson model by allowing for an overabundance of zero observations.\n\nThe model has two parameters, $\\pi$, the proportion of excess zero observations, and $\\lambda$, the mean of the Poisson distribution.  We assume that observations from this model are generated as follows.  First, a weighted coin with probability $\\pi$ of landing on heads is flipped.  If the result is heads, the observation is zero.  If the result is tails, the observation is generated from a Poisson distribution with mean $\\lambda$.  Note that there are two ways for an observation to be zero under this model:\n\n1. the coin is heads, and\n2. the coin is tails, and the sample from the Poisson distribution is zero.\n\nIf $X$ has a zero-inflated Poisson distribution with parameters $\\pi$ and $\\lambda$, its probability mass function is given by\n\n$$\\begin{align*}\nP(X = 0)\n    & = \\pi + (1 - \\pi)\\ e^{-\\lambda} \\\\\nP(X = x)\n    & = (1 - \\pi)\\ e^{-\\lambda}\\ \\frac{\\lambda^x}{x!} \\textrm{ for } x > 0.\n\\end{align*}$$\n\nIn this post, we will use the parameter values $\\pi = 0.3$ and $\\lambda = 2$.  The probability mass function of the zero-inflated Poisson distribution is shown below, next to a normal Poisson distribution, for comparison.", "cell_type": "markdown", "metadata": {}}, {"execution_count": 2, "cell_type": "code", "source": "from __future__ import division\n\nfrom matplotlib import  pyplot as plt\nimport numpy as np\nfrom scipy import stats\nimport seaborn as sns\nfrom statsmodels.base.model import GenericLikelihoodModel", "outputs": [], "metadata": {"collapsed": true, "trusted": true}}, {"execution_count": 3, "cell_type": "code", "source": "np.random.seed(123456789)", "outputs": [], "metadata": {"collapsed": true, "trusted": true}}, {"execution_count": 4, "cell_type": "code", "source": "pi = 0.3\nlambda_ = 2.", "outputs": [], "metadata": {"collapsed": true, "trusted": true}}, {"execution_count": 5, "cell_type": "code", "source": "def zip_pmf(x, pi=pi, lambda_=lambda_):\n    if pi < 0 or pi > 1 or lambda_ <= 0:\n        return np.zeros_like(x)\n    else:\n        return (x == 0) * pi + (1 - pi) * stats.poisson.pmf(x, lambda_)", "outputs": [], "metadata": {"collapsed": false, "trusted": true}}, {"execution_count": 6, "cell_type": "code", "source": "fig, ax = plt.subplots(figsize=(8, 6))\n\nxs = np.arange(0, 10);\n\npalette = sns.color_palette()\n\nax.bar(2.5 * xs, stats.poisson.pmf(xs, lambda_), width=0.9, color=palette[0], label='Poisson');\nax.bar(2.5 * xs + 1, zip_pmf(xs), width=0.9, color=palette[1], label='Zero-inflated Poisson');\n\nax.set_xticks(2.5 * xs + 1);\nax.set_xticklabels(xs);\nax.set_xlabel('$x$');\n\nax.set_ylabel('$P(X = x)$');\n\nax.legend();", "outputs": [{"output_type": "display_data", "data": {"image/png": "iVBORw0KGgoAAAANSUhEUgAAAf0AAAGCCAYAAADquBqcAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3X+cXHV97/HXkhjIJjEksEEIK9iQ+xGQH+FasOWnCBi8\nGFqtxVwR+U2rKNdeQa1AC4o04A/S0lLkl2ixcKuAwZJCUVRaFU1FFImfmHCTbhIoe8kSQhLMr71/\nzCRulk12sjs7Zyfn9Xw88tiZM+d7znuHkPecc+ac09Ld3Y0kSdr57VJ0AEmS1BiWviRJJWHpS5JU\nEpa+JEklYelLklQSlr4kSSUxsoiVRsR04AZgBHBrZs7axny/C/wQOCMzv7EjYyVJ0tYavqUfESOA\nG4HpwEHAzIg4cBvzzQL+ZUfHSpKkVyti9/6RwMLMXJyZ64G7gdP7mO/DwNeBzgGMlSRJvRRR+pOB\njh7Pl1anbRERk6mU+U3VSZsvG9jvWEmS1LciSr+W6/7eAHwiM7uBluqfWsdKkqQ+FPFFvmVAe4/n\n7VS22Hv678DdEQGwJ3BqRKyvcexWNmzY2D1y5IjBZpYkqZm09DWxiNKfB0yNiP2B5cAZwMyeM2Tm\n72x+HBF3AA9k5pyIGNnf2N66utbUJXRb2zg6O1fVZVmN0myZmy0vmLkRmi0vmLkRmi0vNDZzW9u4\nPqc3fPd+Zm4ALgYeAp4G7snM+RFxUURcNJCxQ51ZkqSdQSHn6WfmXGBur2k3b2Pec/obO9TWrVvH\nggULWLHi5QEvo719P0aNGlXHVJIk7ZhCSr/ZdHQs4bI5VzJmG7tL+rO6cxXXzbiaKVOm1jmZJEm1\ns/RrNKZtHOP22b3oGJIkDZjX3pckqSQsfUmSSsLSlySpJDymL0nD1Lp16+joWFLXZfZ3JtFxxx3J\nlCkHsHHjRvbb7w1cfvlfsuuuu/U577/92/dZvPgZzjzz7Lpm1NCx9CVpmOroWMIl18+hdfykuixv\nzcrnmX3pjO2eSbTrrrtxxx1fA+Dqq6/g/vu/wRlnvK/PeY855jiOOea4umRTY1j6kjSMtY6fxNgJ\nxdxX7NBDD2PRokW89NJLXHvtVSxfvpzddtuNyy77FFOmHMCDDz5A5nw++tHL+M53HuHLX76FXXYZ\nwdixY7nxxi/xzDOLuPbaq9mwYT2bNnXz2c9ez+TJ+3L33f/Agw8+AMBpp/0Bf/zHM3n22eV87GMf\n4dBDp/HUU0/S1jaJa6/9PLvuumshv/vOytKXJL3Khg0bePzxH3LUUb/Pbbf9PREHcu21n+enP53H\nZz5z5Za9AS0tlUu833nnrXzhC3/LnnvuyerVlQuZzZlzL+95z0xOOWU6GzZsYOPGjfzqV/OZO/db\n3HLLnWza1M2FF36AadOOYOzYcSxd2sFVV13Lxz/+Ka688pN873vf4ZRTTi3sPdgZ+UU+SdIW69b9\nhnPO+Z9ccMFZ7LXX3px22un84hdP8va3vwOAI454MytXrmTNmtUAdHdXbn56yCGHcc01f8EDD9zP\nxo0bATj44EP46ldv56677uS5555l11135ec//xnHHfdWdt11N0aPHs3xx5/Ik08+QUtLC3vvPZkD\nDqgceoh4I88+u7yAd2Dn5pa+JGmLUaN23bIV39Pmcv+trW/i9rGPfZKnn36KH/7w3znvvPdz221f\n5eSTp3PwwYfwgx88xsc+dgmXXfbnW/YM9Fzu5mmjRr1my/RddhnBxo2/qc8vpS3c0pckbdehh07j\n4Ycrtzz56U/nsfvuE2htbd1qnmXLlnLQQW/ivPMuYvfdd+f5559n+fJl7L33PvzRH72XY489nkWL\nFnLYYYfz/e9/l9/85hXWrl3LY499l0MPndbHhwoNBbf0JWkYW7Py+YYuq/eWOMC5517ItddezQc+\nMJPRo0dz+eV/uWXezfP/3d/NZunSDrq7u3nzm4/kgAOm8g//8GUeeuhBRo4cyR577MlZZ53LuHHj\neMc7TuOCCz4AwDvf+YdMnfrfePbZ5a9ad19ZNDgtO/unq87OVYP+BRct+jVX/fD6AV97f9XyF/mL\n37u04Tfcabb7TTdbXjBzIzRbXqhf5kaep99s73Oz5YXGZm5rG9fnJya39CVpmBo1apR351RdeUxf\nkqSSsPQlSSoJS1+SpJKw9CVJKglLX5KkkvDb+5I0TBVxa93vfe9RvvzlW7aatmjRQq6/fjZHHfV7\ndc1S6615//ZvZ/OjH/07v/d7R7PbbqMZPbqVj3zkT7c5/2OPfZf29v3Yf/837FCek08+ln/918de\nNX1nut2wpS9Jw1RHxxIum3MlY9rG1WV5qztXcd2Mq7d7GuDxx7+V449/65bn3/zmvTzyyEM1F/7m\na7/UcmGdWm/N+8AD9zF37qO0tLRw++1f6nf+73//uxx99LE7XPq9Ly282c50u2FLX5KGsTFt4wZ8\nYbDB+s//XMKdd97G3//97Vumfe1rX+HRRx9h3br1HHfcCZx33kU8++xy/uzPLubggw8hcz7XX//X\nfOMb9/D44z+gpaWFs846j7e97eRXLb/nrXmvueYvGTNmLJlP88ILL/DBD36EE054Gx//+EdZu3Yt\n5577Ps4885ytxs+Zcx8PPHAf69dvYN999+WKK65mwYLk3//9MX72sye4887buOaa6+nu7uYLX7iO\nF1/sYrfdduPjH/8Ur3/9/ixfvoyrrrqcV15Zy9FH11bUg7nd8D33/GPhtxu29CVJr7Jhwwauuupy\nLr74o0yatBcAP/7xj1i6tINbbvkKmzZt4hOf+N88+eQTTJq0F8uWLeWKK67moIPexHe/+20WLlzA\nnXfezYsvdnH++Wdx+OHT2GOPPbdaR++9AStWvMBNN93O4sX/l0984s844YS3MWvWFzn55OO2bGn3\n3NI/4YQTmTHjDwG45Zab+Na3vsm7330GxxxzHEcffSzHH38iAJdc8qdceumfs+++7fzyl0/x+c/P\nYvbsm5g9+3O8613v4e1vfwf33vtPNb0nzX67YUtfkvQqt9xyE1OmHMCJJ560ZdqPf/wjfvKTxznn\nnP8JwNq1r7B0aQeTJu3FXnvtzUEHvQmAX/ziSU4+eTotLS1MmDCRww8/gvnzn97ubu+WlhaOPfZ4\nAPbf/w2sWLGi34yLFi3klltuYvXql1mzZu1WhyA2H2ZYs2YNTz31c6644uNbXlu/fgMATz31cz77\n2c8B8Pa3n8pNN/1Nn+vZfLthgMMOO4LTTjudCy/8ANdccz3Q/+2GTzzx5C2HTA4++BC+8pXb6ez8\nL44//kT23bd9q9sNA1tuN3zMMcfX/XbDlr4kaSs//ek8vv/9R7n99rte9dqZZ57N6ae/a6tpzz67\nnNGjt/5iW1/3dbn33n/igQfuo6Wlheuum/2q11/zmt/eWnd794XZvIPgs5+9ir/6qy8wZcoBzJ37\nLZ544j96zNNSXc4mxo4d1+ftgmtVr9sN33//fYXfbthT9iRJW1SOVV/N5ZdfzejRo7d67aij3sI/\n//Mc1q5dC0Bn5/N0dXW9ahmHHjqNb3/7X9m0aRNdXV08+eQTHHzwm3jXu97DHXd8jdtvv4s999xz\nwLfT3Txs7do1TJy4Bxs2bOChhx7c8nprayurV1e2useMGcs+++zDo48+Uh3bzcKFvwYqW+Lf/vbD\nADz88L/sUIaB3G74ueeeK/x2w27pS9IwtrqOd2WrZVnf/OY3ePHFLj73uc9uNf397z+XE088icWL\nF/Mnf1L5Ql1raytXXPHprW6xC5UzAH75y59z9tkzaWlp4YMfvIQJEya+al29x9XyuPK88vP88/+E\nCy88m913352DD34Ta9asAeBtbzuFWbOu4etfv4fPfGYWV175GT73ub/izjtvZ8OGDZx00ikccMBU\nLrnkY1x11eXcddedHHPM8ds846Betxt+4xvfyBe/+DeF3m7YW+vWwFvrNkaz5QUzN0Kz5QVvrdsI\nzZYXvLWuJGk7vLWu6s1j+pIklYSlL0lSSVj6kiSVRCHH9CNiOnADMAK4NTNn9Xr9dOBqYFP1z6WZ\n+Z3qa4uBl4CNwPrMPLJxySVJal4NL/2IGAHcCJwELAN+EhFzMnN+j9keycxvVuc/BLgPOKD6Wjdw\nQmb2f7kmSZK0RRG7948EFmbm4sxcD9wNnN5zhsxc3ePpWOD/9VrG4E5UlCSphIrYvT8Z6OjxfClw\nVO+ZIuIPgGuBvYFTerzUDTwSERuBmzPzlt5jJUnSqxWxpV/TxXIy8/7MPBB4J/DVHi8dnZnTgFOB\nD0XEsUOQUZKknU4RW/rLgPYez9upbO33KTMfi4iREbFHZr6Qmc9Wp3dGxH1UDhc8tq3xEya0MnLk\niEEF7uoaO6jxABMnjqWtbdygl7OjiljnYDRbXjBzIzRbXjBzIzRbXig+cxGlPw+YGhH7A8uBM4CZ\nPWeIiCnAM5nZHRFHAGTmCxHRCozIzFURMYbKbv+rtreyrq41gw68YsXLdVlGoy8Z2WyXqWy2vGDm\nRmi2vGDmRmi2vNDwy/D2Ob3hpZ+ZGyLiYuAhKqfs3ZaZ8yPiourrNwPvBs6KiPXAy8B7q8NfB9wb\nEZuz35WZDzf6d5AkqRkVcp5+Zs4F5vaadnOPx9cB1/Ux7hng8CEPKEnSTsgr8kmSVBKWviRJJWHp\nS5JUEpa+JEklYelLklQSlr4kSSVh6UuSVBKWviRJJWHpS5JUEpa+JEklYelLklQSlr4kSSVh6UuS\nVBKWviRJJWHpS5JUEpa+JEklYelLklQSlr4kSSVh6UuSVBKWviRJJWHpS5JUEpa+JEklYelLklQS\nlr4kSSVh6UuSVBKWviRJJWHpS5JUEpa+JEklYelLklQSlr4kSSVh6UuSVBKWviRJJWHpS5JUEiOL\nWGlETAduAEYAt2bmrF6vnw5cDWyq/rk0M79Ty1hJktS3hm/pR8QI4EZgOnAQMDMiDuw12yOZeVhm\nTgPOBr60A2MlSVIfiti9fySwMDMXZ+Z64G7g9J4zZObqHk/HAv+v1rGSJKlvRezenwx09Hi+FDiq\n90wR8QfAtcDewCk7MlaSJL1aEVv63bXMlJn3Z+aBwDuBr0ZEy9DGkiRp51bElv4yoL3H83YqW+x9\nyszHImIkMLE6X81jASZMaGXkyBEDTwt0dY0d1HiAiRPH0tY2btDL2VFFrHMwmi0vmLkRmi0vmLkR\nmi0vFJ+5iNKfB0yNiP2B5cAZwMyeM0TEFOCZzOyOiCMAMvOFiFjZ39jeurrWDDrwihUv12UZnZ2r\nBr2cHdHWNq7h6xyMZssLZm6EZssLZm6EZssLjc28rQ8XDd+9n5kbgIuBh4CngXsyc35EXBQRF1Vn\nezfwi4h4ApgNvHd7Yxv9O0iS1IwKOU8/M+cCc3tNu7nH4+uA62odK0mS+ucV+SRJKglLX5KkkrD0\nJUkqCUtfkqSSsPQlSSoJS1+SpJKw9CVJKglLX5KkkrD0JUkqCUtfkqSSsPQlSSoJS1+SpJKw9CVJ\nKglLX5KkkrD0JUkqCUtfkqSSsPQlSSoJS1+SpJKw9CVJKomRRQfQ1tatW0dHx5JBLaO9fT9GjRpV\np0TbV4+8UMksSRpalv4w09GxhEuun0Pr+EkDGr9m5fPMvnQGU6ZMrXOyvg02L/w28+TJe9QxmSSp\nN0t/GGodP4mxEyYXHaNmzZZXksrKY/qSJJWEpS9JUklY+pIklYSlL0lSSVj6kiSVhKUvSVJJWPqS\nJJWEpS9JUklY+pIklYSlL0lSSVj6kiSVRCHX3o+I6cANwAjg1syc1ev19wGXAS3AKuBPM/Pn1dcW\nAy8BG4H1mXlk45JLktS8Gr6lHxEjgBuB6cBBwMyIOLDXbM8Ax2XmocCngS/1eK0bOCEzp1n4kiTV\nrogt/SOBhZm5GCAi7gZOB+ZvniEzf9hj/seBfXsto2WIM0qStNMp4pj+ZKCjx/Ol1Wnbch7wYI/n\n3cAjETEvIi4YgnySJO2Uiij97lpnjIi3AucCH+8x+ejMnAacCnwoIo6tcz5JknZKRezeXwa093je\nTmVrfysRcShwCzA9M7s2T8/MZ6s/OyPiPiqHCx7b1somTGhl5MgRgwrc1TV2UOMBJk4cS1vbuIav\nq5Z1DkY98kIlMwx93qFg5qHXbHnBzI3QbHmh+MxFlP48YGpE7A8sB84AZvacISJeD9wLnJmZC3tM\nbwVGZOaqiBgDnAJctb2VdXWtGXTgFSterssyOjtXNXRdbW3jalrnYNdVz+UMdd56a8R7XG/NlrnZ\n8oKZG6HZ8kJjM2/rw0XDSz8zN0TExcBDVE7Zuy0z50fERdXXbwauBCYAN0UE/PbUvNcB91anjQTu\nysyHG/07SJLUjAo5Tz8z5wJze027ucfj84Hz+xj3DHD4kAeUJGkn5BX5JEkqCUtfkqSSsPQlSSoJ\nS1+SpJKw9CVJKglLX5KkkrD0JUkqCUtfkqSSsPQlSSoJS1+SpJKw9CVJKglLX5KkkrD0JUkqCUtf\nkqSSsPQlSSoJS1+SpJKw9CVJKomRAxkUEa8B9gBezMxX6htJkiQNhZpLPyIOAd4PtALrgNXA+IgA\nWAHcnJnPDkVISZI0eDWVfkS8H3gF+ERmburj9d2AmRGxODMfrXNGSZJUB/2WfkSMBh7Z3lZ8dRf/\nHRHxO/UMJ0mS6qff0s/MtcDaWhaWmc8MOpEkSRoSO/zt/Yi4PiKOrj4+pvqlPkmSNMwN5JS9XwLz\nq48fB/64fnEkSdJQGcgpe/sBX4mIR4AfAhPrG0mSJA2FgWzpLwbOBZYB5wC71jOQJEkaGgMp/TXA\npsz8J2AWNX7JT5IkFWuHS79a9q+rPn0tlYv1SJKkYW5Al+HNzKeqP58EnqxrIkmSNCQGdMOdiDi7\n+vOcuqaRJElDZqB32Rtf/fnaegWRJElDy1vrSpJUEpa+JEklMaAv8g1WREwHbgBGALdm5qxer78P\nuAxoAVYBf5qZP69lrCRJ6lvDt/QjYgRwIzAdOIjKLXkP7DXbM8BxmXko8GngSzswVpIk9aGILf0j\ngYWZuRggIu4GTue31/MnM3/YY/7HgX1rHStJkvo20C39B6o/vzWAsZOBjh7Pl1anbct5wIMDHCtJ\nkqoGenGeZ6o/Fw1geHetM0bEW6lc5//oHR272YQJrYwcOWJHh22lq2vsoMYDTJw4lra2cQ1fVy3r\nHIx65IVKZhj6vEPBzEOv2fKCmRuh2fJC8ZlrKv2ImDLAgu/LMqC9x/N2Klvsvdd5KHALMD0zu3Zk\nbE9dXWsGFRZgxYqX67KMzs5VDV1XW9u4mtY52HXVczlDnbfeGvEe11uzZW62vGDmRmi2vNDYzNv6\ncFHrlv7/Aj5cpyzzgKkRsT+wHDgDmNlzhoh4PXAvcGZmLtyRsZIkqW+1HtM/NSKO7uuFiPjojqww\nMzcAFwMPAU8D92Tm/Ii4KCIuqs52JTABuCkinoiIH29v7I6sX5Kksqp1S/8E4PURcVpmfgu2HG8/\nH/hD4Is7stLMnAvM7TXt5h6Pz68uu6axkiSpf7WW/m8y8wcR8baI+DKVL9a9DNwKLB6ibJIkqY5q\nLf2vRMQ64FjgG1SOt1+emeuHLJkkSaqrWo/pTwG+Cbw+My8ArgI+GBHjtz9MkiQNF7Vu6V+Wmfdv\nfpKZayLib4A/iYjfz8wzhyaeBmrdunUsWLBgUKfUtbfvx6hRo+qYatvWr1/fVHklqRnVVPo9C7/H\ntE0RcROVy+BqmOnoWMJlc65kzAAvBLG6cxXXzbiaKVOm1jlZ3559djmfmvvppskrSc1oUNfez8zu\niJhdrzCqrzFt4xi3z+5Fx6hZs+WVpGbT7zH9iNgtIt6yrdcz88Ee855Yr2CSJKm++t3Sz8xXImJ9\nRHwMeDAzn+75ekTsAhxF5Vz+u4ckpSRJGrRaj+n/R0SsB/4oItqA3YARVG6A8xLwncy8duhiSpKk\nwar1hjufAN5N5b7291dP25MkSU2k1vP0d8nM383MvYFHI2LGUIaSJEn1V+u391/Y/CAz/09EnDtE\neaQht27dOjo6lgxqGV4TQFIzqrX0IyLGZObq6vPV251bGsY6OpZwyfVzaB0/aUDj16x8ntmXzvCa\nAJKaTq2lPwM4IyJeAH4MdEfEU5n5y4g4MTO/M3QRpfprHT+JsRMmFx1Dkhqq1mP6H8nMycAfAT+g\n8s39eyPiOeDGoQonSZLqp9ZT9h6s/lwALABuB4iIScBnhiydJEmqm1q39PuUmc8Df1enLJIkaQgN\nqvQBMvNn9QgiSZKG1qBLX5IkNQdLX5KkkrD0JUkqCUtfkqSSsPQlSSoJS1+SpJKw9CVJKglLX5Kk\nkrD0JUkqCUtfkqSSsPQlSSoJS1+SpJKw9CVJKglLX5KkkhhZxEojYjpwAzACuDUzZ/V6/Y3AHcA0\n4FOZ+fkery0GXgI2Ausz88gGxZYkqak1vPQjYgRwI3ASsAz4SUTMycz5PWZ7Afgw8Ad9LKIbOCEz\nVwx5WEmSdiJF7N4/EliYmYszcz1wN3B6zxkyszMz5wHrt7GMliHOKEnSTqeI0p8MdPR4vrQ6rVbd\nwCMRMS8iLqhrMkmSdmJFlH73IMcfnZnTgFOBD0XEsXXIJEnSTq+IL/ItA9p7PG+nsrVfk8x8tvqz\nMyLuo3K44LFtzT9hQisjR44YYNSKrq6xgxoPMHHiWNraxjVsXfVQS+Z65AUYP7618jdjEBr9Hm9e\nVy3rHG6aLXOz5QUzN0Kz5YXiMxdR+vOAqRGxP7AcOAOYuY15tzp2HxGtwIjMXBURY4BTgKu2t7Ku\nrjWDDrxixct1WUZn56qGraseaslcr3WtXFmf/06NfI87O1fR1jaupnUOJ82WudnygpkbodnyQmMz\nb+vDRcNLPzM3RMTFwENUTtm7LTPnR8RF1ddvjojXAT8BXgtsiohLgIOAScC9EbE5+12Z+XCjfwdJ\nkppRIefpZ+ZcYG6vaTf3ePwcWx8C2Oxl4PChTSdJ0s7JK/JJklQSlr4kSSVh6UuSVBKWviRJJWHp\nS5JUEpa+JEklYelLklQSlr4kSSVRyMV5GmnRol8Panx7+351SiJJUrF2+tK/5Po5tI6fNKCxa1Y+\nz+xLZ9Q5kSRJxdjpS791/CTGTphcdAxJkgrnMX1JkkrC0pckqSQsfUmSSsLSlySpJHb6L/JJQ2Xd\nunUsWLCAFSteHvAy2tv3Y9SoUXVMJUnbZulLA9TRsYTL5lzJmLZxAxq/unMV1824milTptY5mST1\nzdKXBmFM2zjG7bN70TEkqSYe05ckqSQsfUmSSsLSlySpJCx9SZJKwtKXJKkkLH1JkkrC0pckqSQs\nfUmSSsLSlySpJCx9SZJKwtKXJKkkLH1JkkrC0pckqSQsfUmSSsLSlySpJEYWsdKImA7cAIwAbs3M\nWb1efyNwBzAN+FRmfr7WsZIkqW8N39KPiBHAjcB04CBgZkQc2Gu2F4APA58bwFhJktSHInbvHwks\nzMzFmbkeuBs4vecMmdmZmfOA9Ts6VpIk9a2I0p8MdPR4vrQ6bajHSpJUakWUfndBYyVJKrUivsi3\nDGjv8bydyhb7UI8dkIkTx9ZtOW1t4/qdr6tr8OtrZOZ65AUYP7618l93EHbW93goFLHOwWi2vGDm\nRmi2vFB85iJKfx4wNSL2B5YDZwAztzFvyyDG1sWKFS/XbTmdnasasr5GZq7XulauXDPoZeys73G9\ntbWNa/g6B6PZ8oKZG6HZ8kJjM2/rw0XDSz8zN0TExcBDVE67uy0z50fERdXXb46I1wE/AV4LbIqI\nS4CDMvPlvsY2+neQJKkZFXKefmbOBeb2mnZzj8fPsfVu/O2OlSRJ/fOKfJIklYSlL0lSSVj6kiSV\nRCHH9CXtmHXr1tHRsWRQy2hv349Ro0bVKZGkZmTpS02go2MJl1w/h9bxkwY0fs3K55l96QymTJla\n52SSmomlLzWJ1vGTGDvBq05LGjiP6UuSVBKWviRJJWHpS5JUEpa+JEklYelLklQSlr4kSSVh6UuS\nVBKWviRJJWHpS5JUEpa+JEklYelLklQSlr4kSSVh6UuSVBKWviRJJWHpS5JUEpa+JEklYelLklQS\nlr4kSSVh6UuSVBKWviRJJTGy6ACSGmfdunUsWLCAFSteHvAy2tv3Y9SoUXVMJalRLH2pRDo6lnDZ\nnCsZ0zZuQONXd67iuhlXM2XK1Donk9QIlr5UMmPaxjFun92LjiGpAB7TlySpJCx9SZJKwtKXJKkk\nLH1JkkqikC/yRcR04AZgBHBrZs7qY56/Bk4F1gBnZ+YT1emLgZeAjcD6zDyyQbElSWpqDd/Sj4gR\nwI3AdOAgYGZEHNhrnncAB2TmVOBC4KYeL3cDJ2TmNAtfkqTaFbF7/0hgYWYuzsz1wN3A6b3mmQHc\nCZCZjwO7R8RePV5vaUhSSZJ2IkWU/mSgo8fzpdVptc7TDTwSEfMi4oIhSylJ0k6miNLvrnG+bW3N\nH5OZ06gc7/9QRBxbn1iSJO3civgi3zKgvcfzdipb8tubZ9/qNDJzefVnZ0TcR+VwwWNDFXbixLF1\nW05bDZc+7eoa/PoambkeeQHGj2+t/hceuJ31PYbmzFxPjV5fPZh56DVbXig+cxGlPw+YGhH7A8uB\nM4CZveaZA1wM3B0RbwFezMz/iohWYERmroqIMcApwFVDGXYwNybpvZzOzlUNWV8jM9drXStXrhn0\nMnbW97he62t05nppaxvX0PXVg5mHXrPlhcZm3taHi4bv3s/MDVQK/SHgaeCezJwfERdFxEXVeR4E\nnomIhcDNwAerw18HPBYRPwMeB76VmQ83+neQJKkZFXKefmbOBeb2mnZzr+cX9zHuGeDwoU0nSdLO\nySvySZJUEpa+JEklYelLklQSlr4kSSVh6UuSVBKWviRJJWHpS5JUEoWcpy9p57Zu3To6OpYMahnt\n7fsxatSoOiWSBJa+pCHQ0bGES66fQ+v4SQMav2bl88y+dAZTpkytczKp3Cx9SUOidfwkxk7ofdds\nSUXymL4kSSVh6UuSVBKWviRJJWHpS5JUEpa+JEklYelLklQSnrInadhat24dCxYsYMWKlwe8DC/y\nI/2WpS97RzoKAAAHUUlEQVRp2OroWMJlc65kTNu4AY1f3bmK62Zc7UV+pCpLX9KwNqZtHOP22b3o\nGNJOwWP6kiSVhKUvSVJJWPqSJJWEpS9JUklY+pIklYSlL0lSSVj6kiSVhKUvSVJJeHEeSaJyyd+O\njiWDWoaX/NVwZ+lLEpVL/l5y/Rxax08a0Pg1K59n9qUzvOSvhjVLX5KqWsdPYuyEyUXHkIaMpS9J\ndeSdATWcWfqSVEfeGVDDmaUvSXXmnQE1XBVS+hExHbgBGAHcmpmz+pjnr4FTgTXA2Zn5RK1jJUnS\nqzW89CNiBHAjcBKwDPhJRMzJzPk95nkHcEBmTo2Io4CbgLfUMlaSysLTDLWjitjSPxJYmJmLASLi\nbuB0oGdxzwDuBMjMxyNi94h4HfCGGsZKUinU6zTD9vb9/PJhSRRR+pOBjh7PlwJH1TDPZGCfGsZK\nUmnU4zTDRn35sF57JgA/pAxQEaXfXeN8LfVY2ZqVz9dl7OrOVQNezo6ObbbMg8n72/F7+x7vwHp3\nVLP9vfA9Hth6d9Rg/9/dUR0dS7jwilvZbezEAY1/5eUVfOnT5wPwods+yugJYwa0nLVdq/nb875Y\n8xkSixb9ekDr2Wzzegb7QaUeZ3S0dHfX2sH1ERFvAf4yM6dXn38S2NTzC3kR8ffAdzPz7urzXwHH\nU9m9v92xkiSpb0XccGceMDUi9o+IUcAZwJxe88wBzoItHxJezMz/qnGsJEnqQ8NLPzM3ABcDDwFP\nA/dk5vyIuCgiLqrO8yDwTEQsBG4GPri9sY3+HSRJakYN370vSZKKUcTufUmSVABLX5KkkrD0JUkq\nCW+4U4Nmu95/RNwO/A/g+cw8pOg8/YmIduArwCQq13H4Umb+dbGpti8idgO+B+wKjAK+mZmfLDZV\n/6qXsp4HLM3Mdxadpz8RsRh4CdgIrM/MIwsN1I+I2B24FTiYyt/lczPzR8Wm2raICODuHpN+B7ii\nCf7/+yRwJrAJ+AVwTmb+pthU2xYRlwDnU7n+zC2ZObuoLG7p96PH9f6nAwcBMyPiwGJT9esOKnmb\nxXrgo5l5MPAW4EPD/T3OzFeAt2bm4cChwFsj4piCY9XiEipnvjTLN3i7gRMyc9pwL/yq2cCDmXkg\nlb8Xw/rsoqyYlpnTgP9O5QZn9xUca7siYn/gAuCI6kbNCOC9hYbajoh4E5XC/13gMOC0iJhSVB5L\nv39b7hWQmeupfCo+veBM25WZjwFdReeoVWY+l5k/qz5+mco/lPsUm6p/mbmm+nAUlX94VhQYp18R\nsS/wDipbonW54mWDNEXWiBgPHJuZt0PlFOPMXFlwrB1xErAoMzv6nbNYL1HZUGiNiJFAK5UbsA1X\nbwQez8xXMnMjlT2E7yoqjLv3+1fLvQJUJ9VP8dOAxwuO0q+I2AX4KTAFuCkzny44Un++CFwKvLbo\nIDugG3gkIjYCN2fmLUUH2o43AJ0RcQeVLbr/AC7p8eFwuHsv8LWiQ/QnM1dExOeB/wTWAg9l5iMF\nx9qep4BrImIi8AqVQ68/LiqMW/r9a5bdoE0vIsYCX6fyD+XAL1DdIJm5qbp7f1/guIg4oeBI2xQR\np1H5jscTNMmWc9XR1V3Pp1I57HNs0YG2YyRwBPB3mXkEsBr4RLGRalO9wuk7gX8qOkt/qrvG/xew\nP5U9gmMj4n2FhtqOzPwVMAt4GJgLPEHluwiFsPT7twxo7/G8ncrWvuooIl4DfAP4h8y8v+g8O6K6\nC/efgTcXnWU7fh+YERH/F/hH4MSI+ErBmfqVmc9Wf3ZSOdY8nI/rL6XyBcmfVJ9/ncqHgGZwKvAf\n1fd5uHsz8IPMfKF6ldZ7qfz9HrYy8/bMfHNmHg+8CGRRWSz9/nm9/yEWES3AbcDTmXlD0XlqERF7\nVr+pTUSMBk6m8gl+WMrMP8/M9sx8A5XduN/JzLOKzrU9EdEaEeOqj8cAp1D5pvawlJnPAR0R8d+q\nk04CfllgpB0xk8qHwWbwK+AtETG6+m/HSVS+nDpsRcSk6s/XA39IgYdRPKbfj8zcEBGbr/c/Arht\nuF/vPyL+kcpdCfeIiA7gysy8o+BY23M0ldNvfh4Rm4vzk5n5LwVm6s/ewJ3V4/q7AF/NzG8XnGlH\nNMNhq72A+ypnlTESuCszHy42Ur8+DNxV3UBYBJxTcJ5+VT9QnUTlG/HDXmY+Wd1LNY/KbvKfAl8q\nNlW/vh4Re1D5AuIHM/OlooJ47X1JkkrC3fuSJJWEpS9JUklY+pIklYSlL0lSSVj6kiSVhKUvSVJJ\nWPqSJJWEpS9JUklY+pIklYSX4ZVUdxHxHir3Od8fWAIcnJmXFhpKkqUvqb4i4k3Ao1T2JN4O/A3w\n60JDSQK89r6kIVLd2p/cLHdOlMrA0pdUVxFxGPAS8HHgDip3QTsqM/+t0GCS3L0vqe5OAdYAC4Ej\ngQOA/1NoIkmAW/qSJJWGp+xJklQSlr4kSSVh6UuSVBKWviRJJWHpS5JUEpa+JEklYelLklQSlr4k\nSSXx/wEKf0WcLgB7iwAAAABJRU5ErkJggg==\n", "text/plain": "<matplotlib.figure.Figure at 0x7faf75d57910>"}, "metadata": {}}], "metadata": {"collapsed": false, "trusted": true}}, {"source": "#### Maximum likelihood estimation", "cell_type": "markdown", "metadata": {}}, {"source": "First we generate 1,000 observations from the zero-inflated model.", "cell_type": "markdown", "metadata": {}}, {"execution_count": 7, "cell_type": "code", "source": "N = 1000", "outputs": [], "metadata": {"collapsed": true, "trusted": true}}, {"execution_count": 8, "cell_type": "code", "source": "inflated_zero = stats.bernoulli.rvs(pi, size=N)\nx = (1 - inflated_zero) * stats.poisson.rvs(lambda_, size=N)", "outputs": [], "metadata": {"collapsed": false, "trusted": true}}, {"execution_count": 9, "cell_type": "code", "source": "fig, ax = plt.subplots(figsize=(8, 6))\n\nax.hist(x, width=0.8, bins=np.arange(x.max() + 1), normed=True);\n\nax.set_xticks(np.arange(x.max() + 1) + 0.4);\nax.set_xticklabels(np.arange(x.max() + 1));\nax.set_xlabel('$x$');\n\nax.set_ylabel('Proportion of samples');", "outputs": [{"output_type": "display_data", "data": {"image/png": "iVBORw0KGgoAAAANSUhEUgAAAfoAAAGCCAYAAAAIZAHlAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3X+0XWdd5/F3cq9B0oTSwq2FcCVY4teWnykSGEuxZWon\ntdK6xuWEDIgglIAtVEdbZRwdQWc0INpCFxja0hFEUsUWAzb2B3RY4YeVDC2gDd8xdlJv0kKvTUyT\nVkzuj/ljn9TTy7337Jt7ds49T96vtbJy9j7Ps8/3Wf3xOc/+8ZxFk5OTSJKkMi3udQGSJKk5Br0k\nSQUz6CVJKphBL0lSwQx6SZIKZtBLklSwwSYPHhFrgauAAeC6zNw4Q7uXAl8G1mXmn7f27QIeAcaB\nw5m5pslaJUkqUWNBHxEDwDXAecAe4CsRsSUzd0zTbiPwV1MOMQmck5l7m6pRkqTSNXnqfg2wMzN3\nZeZhYDNw8TTt3g58Ehid5r1FDdYnSVLxmgz6FcBI2/bu1r7HRcQKqvD/UGtX+zJ9k8AdEbE9Ii5p\nsE5JkorVZNDXWVv3KuBXM3OSavbePoM/KzNXAxcAl0bE2Q3UKElS0Zq8GW8PMNy2PUw1q2/3EmBz\nRAA8HbggIg5n5pbMfBAgM0cj4maqSwHbZvqwsbHxycHBgW7WL0nSQtfxEneTQb8dWBURK4EHgHXA\n+vYGmfkDR15HxA3ApzNzS0QsBQYy80BEnACcD7xrtg/bt++xLpcPQ0PLGR090PXj9lqJ43JM/cEx\n9QfH1D+GhpZ3bNPYqfvMHAMuA24F7gVuzMwdEbEhIjZ06H4qsC0i7gHuAj6Tmbc1VaskSaVq9Dn6\nzNwKbJ2yb9MMbd/Y9vo+4MVN1iZJ0vHAlfEkSSqYQS9JUsEMekmSCmbQS5JUMINekqSCGfSSJBXM\noJckqWAGvSRJBTPoJUkqmEEvSVLBDHpJkgpm0EuSVDCDXpKkghn0kiQVzKCXJKlgBr0kSQUz6CVJ\nKphBL0lSwQx6SZIKZtBLklQwg16SpIIZ9JIkFcyglySpYAa9JEkFM+glSSqYQS9JUsEMekmSCmbQ\nS5JUsMEmDx4Ra4GrgAHguszcOEO7lwJfBtZl5p/Ppa8kSZpZYzP6iBgArgHWAmcA6yPi9BnabQT+\naq59JUnS7Jqc0a8BdmbmLoCI2AxcDOyY0u7twCeBlx5F38e97m1X8qSnPKtbtQMwMLiY8bGJ7h7z\n0D/x4T94d1ePKUnSTJoM+hXASNv2buBl7Q0iYgVVgL+KKugn6/adanDZCiZPev48S36isa4erTLw\n6IzfVSRJ6romb8ab7NyEq4BfzcxJYFHrT92+kiSpgyZn9HuA4bbtYaqZebuXAJsjAuDpwAURcbhm\n3ydYPNAfDxAMDg4wNLS812UsiBq6zTH1B8fUHxxTOZoM+u3AqohYCTwArAPWtzfIzB848joibgA+\nnZlbImKwU9+pJsYn+uJZwbGxcUZHD/S0hqGh5T2vodscU39wTP3BMfWPOl9eGsvGzBwDLgNuBe4F\nbszMHRGxISI2HE3fpmqVJKlUjT5Hn5lbga1T9m2aoe0bO/WVJElz0w9nuyVJ0lEy6CVJKphBL0lS\nwQx6SZIKZtBLklQwg16SpIIZ9JIkFcyglySpYAa9JEkFM+glSSqYQS9JUsEMekmSCmbQS5JUMINe\nkqSCGfSSJBXMoJckqWAGvSRJBTPoJUkqmEEvSVLBDHpJkgpm0EuSVDCDXpKkghn0kiQVzKCXJKlg\nBr0kSQUz6CVJKphBL0lSwQx6SZIKZtBLklSwwSYPHhFrgauAAeC6zNw45f2LgXcDE60/V2Tm51rv\n7QIeAcaBw5m5pslaJUkqUWNBHxEDwDXAecAe4CsRsSUzd7Q1uyMz/6LV/gXAzcBzW+9NAudk5t6m\napQkqXRNnrpfA+zMzF2ZeRjYDFzc3iAzH23bXAb805RjLGqwPkmSitfkqfsVwEjb9m7gZVMbRcRP\nAr8DPAM4v+2tSeCOiBgHNmXmtQ3WKklSkZoM+sk6jTLzU8CnIuJs4GNAtN46KzMfjIgh4PaI+GZm\nbpvpOIsH+uO+wsHBAYaGlve6jAVRQ7c5pv7gmPqDYypHk0G/Bxhu2x6mmtVPKzO3RcRgRDwtMx/O\nzAdb+0cj4maqSwEzBv3E+ERfPEIwNjbO6OiBntYwNLS85zV0m2PqD46pPzim/lHny0uT2bgdWBUR\nKyNiCbAO2NLeICJOi4hFrddnAmTmwxGxNCKWt/afQHVK/xsN1ipJUpEam9Fn5lhEXAbcSvV43fWZ\nuSMiNrTe3wT8FPD6iDgMHARe0+p+KnBTRByp8eOZeVtTtUqSVKpGn6PPzK3A1in7NrW9fg/wnmn6\n3Qe8uMnaJEk6HvTDZW1JknSUDHpJkgpm0EuSVDCDXpKkghn0kiQVzKCXJKlgBr0kSQUz6CVJKphB\nL0lSwQx6SZIKZtBLklQwg16SpIIZ9JIkFcyglySpYAa9JEkFM+glSSqYQS9JUsEMekmSCmbQS5JU\nMINekqSCGfSSJBXMoJckqWAGvSRJBTPoJUkqmEEvSVLBDHpJkgpm0EuSVDCDXpKkgg02efCIWAtc\nBQwA12XmxinvXwy8G5ho/bkiMz9Xp68kSepsTjP6iPi+iHh5RCyq0XYAuAZYC5wBrI+I06c0uyMz\nX5SZq4E3AB+eQ19JktRBxxl9RHwBuBBYBHwV2A/cAvxyh65rgJ2Zuat1nM3AxcCOIw0y89G29suA\nf6rbV5IkdVZnRr8sM/cDPwF8HHgB1Uy7kxXASNv27ta+J4iIn4yIHcBW4B1z6StJkmZXJ+if1Pr7\nVVSn2seBsRr9JusUkJmfyszTgVcDH6tzWUCSJNVT52a8OyPiXuB7gLdGxEnUC/o9wHDb9jDVzHxa\nmbktIgaBk1vtavcFWDzQHw8QDA4OMDS0vNdlLIgaus0x9QfH1B8cUznqBP1lwIuAf8jMQxHxFOCS\nGv22A6siYiXwALAOWN/eICJOA+7LzMmIOBMgMx+OiP2d+k41MT7RF88Kjo2NMzp6oKc1DA0t73kN\n3eaY+oNj6g+OqX/U+fLSMRszcwJ4GvD61q4B4NGZezzeb4zqS8KtwL3AjZm5IyI2RMSGVrOfAr4R\nEXcDVwOvma1vx9FIkqQnqHPX/TuBC4BnUD3ytgT4CPCKTn0zcyvVTXbt+za1vX4P8J66fSVJ0tzU\nOdu9HjgPOAiQmSPAU5osSpIkdUedoP+XzDzUeCWSJKnr6tyM948RcTY8vmLdO4G/bbQqSZLUFXWC\n/h3AR4HnA48B24DXNlmUJEnqjo5Bn5kPAj8WEScAizOzvOcTJEkq1IxBHxFnzLAfgMy8t6GaJElS\nl8w2o7+F2ZexfU6Xa5EkSV02Y9Bn5spjWIckSWpAnZvxiIjnA+dSzfDvzMy/a7QqSZLUFR2fo4+I\nS6mWon0B1Zr3t0bEzzddmCRJmr86M/pfAFZn5kMAETEEfAn4YJOFSZKk+auzMt7+IyEPkJmjwP7m\nSpIkSd1SZ0Z/e0RcB1wPLALeANx25PE7H7OTJGnhqhP066luwvv30+wHH7OTJGnBqrMy3spjUIck\nSWpA3cfrlgLPam/vKXtJkha+jkEfEe8A/gewDxhve8tT9pIkLXB1ZvS/CERmPtB0MZIkqbvqPF63\n25CXJKk/1ZnRvzsirgf+EvgO1SN2k5l5S6OVSZKkeasT9BcCPwGs4onX6A16SZIWuDpB/5PAysz8\nl6aLkSRJ3VXnGv0/AGNNFyJJkrqvzoz+74HPRsSngH9t7ZvMTH/URpKkBa5O0H8vcB/Vz9RKkqQ+\nUmcJ3DccgzokSVID6i6BG8CLqGb3AGTmR5sqSpIkdUedJXAvB94CPBP4G+Bs4POAQS9J0gJXZ0b/\nFuBlwBcy8z9ExPOB/17n4BGxFrgKGACuy8yNU95/LXAl1SI8B4C3ZebXW+/tAh6henb/cGauqfOZ\nkiTp39R5vO47mXkQWBwRizPzb4Ef7NQpIgaAa4C1wBnA+og4fUqz+4BXZuYLgd8CPtz23iRwTmau\nNuQlSTo6dWb0j0bEEuDrwO9GxG7qfUFYA+zMzF0AEbEZuBjYcaRBZn65rf1dVD+F225Rjc+RJEkz\nqBPYlwJLgF8Cnga8EviZGv1WACNt27tb+2byJp64rO4kcEdEbI+IS2p8niRJmqLO43XfaL08CLwp\nIhZl5mSNY9dpA0BEnAv8HHBW2+6zMvPBiBgCbo+Ib2bmtrrHlCRJ9e66/33gXcCjwJ3AmRHx1sz8\nWIeue4Dhtu1hqln91OO/ELgWWJuZ+47sz8wHW3+PRsTNVJcCZgz6xQN1Tk703uDgAENDy3tdxoKo\nodscU39wTP3BMZWjzjX68zLzv0TEhVTh/RqqU+ydgn47sCoiVgIPAOuA9e0NIuL7gZuA12Xmzrb9\nS4GBzDwQEScA51N92ZjRxPhEresQvTY2Ns7o6IGe1jA0tLznNXSbY+oPjqk/OKb+UefLy1yy8UeB\nmzNzDzDRqXFmjgGXAbcC9wI3ZuaOiNgQERtazX4DOAn4UETcHRF/09p/KrAtIu6huknvM5l52xxq\nlSRJ1JvRPxQRfwhcAPxORHwP1XPxHWXmVmDrlH2b2l6/GXjzNP3uA15c5zMkSdLM6szo/zOQwLrW\nNfQVwO83WpUkSeqKOnfdPwT8Qdv2LuB/NVeSJEnqln64f02SJB0lg16SpILNGPQR8Uutv19x7MqR\nJEndNNuM/nWtvz9wLAqRJEndN9vNeP8SEZ8BVkbEn015bzIz/1ODdUmSpC6YLehfDZwHvAD4DE/8\nJbna69hLkqTemTHoM/Nh4MaIeCgz7zyGNUmSpC6pszLe5yPirVSze4DbgGtr/oKdJEnqoTpBvxFY\nDdxAdfr+Z4FVwBUN1iVJkrqgTtCvBc7MzMMAEXEj8FUMekmSFry6C+ZMzvBakiQtYHVm9LcCWyOi\n/dT9rY1WJUmSuqJO0P8K8BbgP7a2bwI+3FhFkiSpa+r8et048KHWH0mS1Ef8URtJkgpm0EuSVDCD\nXpKkgtW5GY+IWAo8q719Zt7bVFGSJKk7OgZ9RFwK/C6wDxhve+s5TRUlSZK6o86M/peB52fm/U0X\nI0mSuqvONfoHDXlJkvpTnRn97RHxHmAz8J0jO71GL0nSwlcn6H+Wan37n56y32v0kiQtcHVWxlt5\nDOqQJEkNqPt43RnAq6hm9p/LzB2NViVJkrqi4814EfEzwO3Ai4DVwB0R8bqmC5MkSfNXZ0Z/BfCS\nzPwWQEScCtwG/HGnjhGxFrgKGACuy8yNU95/LXAl1c/fHgDelplfr9NXkiR1VufxuskjIQ/Qej3Z\nqVNEDADXAGuBM4D1EXH6lGb3Aa/MzBcCv0Xr529r9pUkSR3UmdHfFxHvAjZRzbwvoQroTtYAOzNz\nF0BEbAYuBh6/vp+ZX25rfxfVMru1+kqSpM7qzOjfCvwQ8HXga63XG2r0WwGMtG3vbu2byZuAW46y\nryRJmkadx+u+Daw7imN3PL1/REScC/wccNZc+0qSpJnNGPQRcVZmfjEiLmSa4M3MW6bp1m4PMNy2\nPUw1M5/6OS8ErgXWZua+ufRtt3igP35xd3BwgKGh5b0uY0HU0G2OqT84pv7gmMox24z+DcAXqe66\nn26G3SnotwOrImIl8ADVWYH17Q0i4vuBm4DXZebOufSdamJ8otZ1iF4bGxtndPRAT2sYGlre8xq6\nzTH1B8fUHxxT/6jz5WXGoM/MS1p/n3M0H56ZYxFxGXAr1SNy12fmjojY0Hp/E/AbwEnAhyIC4HBm\nrpmp79HUIUnS8azO79F/ITNf0WnfdDJzK7B1yr5Nba/fDLy5bl9JkjQ3dc52n9C+0XrG/eRmypEk\nSd002814V1Jdn39qRIy2vbUU+HjThUmSpPmb7dT9JuBPqVaou5RqsRyARzJzb9OFSZKk+ZvtZrz9\nEXEQWJSZ9x/DmiRJUpfMeo0+M8eBkyOiH55ckyRJU9RZ6/6vgZsi4k+Ag0d21lgwR5Ik9VidoF9N\ntWDO26bsN+glSVrg6qx1f84xqEOSJDWgzoyeiFgLnEc1s789M29rtCpx6NAhRkaauQdy375l7N17\nsHPDORgefjZLliyZtU2JY5Kkha7OynhXAD8LfILqEbv3RcRHM/O9TRd3PBsZuZ/L37uFpSee0utS\nOnps/0NcfcVFnHbaqlnblTgmSVro6szoXw/8u8w8ABARVwNfAgz6hi098RSWnbSi12V0VYljkqSF\nrNZjc0dCfuprSZK0sNWZ0W+PiBuofjN+EfAmqp+RlSRJC1ydGf3bgYeA9wNXt15f1mRRkiSpO+o8\nXncQ+JVjUIskSeqyOnfdPwX4deBVrV2fBX7La/WSJC18dU7df4Tq9+ffDryj9fqGJouSJEndUedm\nvOdl5ult21+MiB1NFSRJkrqnzoz+gYgYOrIREU8H9jRXkiRJ6pY6M/qHga9FxKepHq+7ENgWEe8F\nJjPzyiYLlCRJR69O0N/b+nPEtVRr3i9q/S1JkhaoOo/X/eYxqEOSJDWgzuN1J1A9Xndea9dtwG9n\n5mNNFiZJkuavzs14HwCeAVwO/ALwTOCaJouSJEndUeca/Usz8wVHNiLii8DXmytJkiR1S61fr4uI\nZW2bJzRUiyRJ6rI6M/o/Br4cEZ+gutN+HfCxRquSJEld0XFGn5kbgSuBpwEnAVdm5nuaLkySJM3f\nrDP6iBgE/iYzzwS2zvXgEbEWuAoYAK5rfWlof/+HqNbNXw38Wma+r+29XcAjwDhwODPXzPXzJUk6\n3s06o8/MMeBgRDx5rgeOiAGqu/PXAmcA6yPi9CnNHqb6sZzfm+YQk8A5mbnakJck6ejUuUb/f4HP\nR8QngUdb+yYz84Md+q0BdmbmLoCI2AxcDDz+gziZOQqMRsSFMxxjUY36JEnSDOoE/SDVErhTZ+Od\nrABG2rZ3Ay+bQ/9J4I6IGAc2Zea1c/x8SZKOe52u0Z9Mdfr97zNz/xyPPd918M/KzAdbv5x3e0R8\nMzO3zfOYkiQdV2YM+ohYR3Wj3AHgSRHxU5n52Tkcew8w3LY9TDWrryUzH2z9PRoRN1NdCpgx6BcP\n1FoSoOcGBwcYGlresd2+fcs6tllITj55WcdxlTimpvX685vgmPqDYyrHbDP6/wb8SGbeExHnAr8J\nzCXotwOrImIl8ADV8/frZ2j7hGvxEbEUGMjMA6219s8H3jXbh02MT9Rb/afHxsbGGR090LHd3r0H\nj0E13bN378GO4ypxTE0aGlre089vgmPqD46pf9T58jJb0I9n5j0AmXlnRPz+XD48M8ci4jLgVqrH\n667PzB0RsaH1/qaIOBX4CvAUYCIiLqe6Q/8U4KaIOFLjxzPztrl8viRJmj3onxQRZ7ReLwK+t22b\nzLx3+m7/JjO3MuX5+8zc1Pb6Wzzx9P4RB4EXdzq+JEma3WxB/2TgL9u2F03Zfk4jFUmSpK6ZMegz\nc+UxrEOSJDWgH+5fkyRJR8mglySpYAa9JEkFM+glSSqYQS9JUsEMekmSCmbQS5JUMINekqSCGfSS\nJBXMoJckqWAGvSRJBTPoJUkqmEEvSVLBZvuZWkk1HDp0iJGR+7t+3H37lrF378GuHnN4+NksWbKk\nq8eUtLAZ9NI8jYzcz+Xv3cLSE0/pdSmzemz/Q1x9xUWcdtqqXpci6Rgy6KUuWHriKSw7aUWvy5Ck\n7+I1ekmSCmbQS5JUMINekqSCGfSSJBXMoJckqWAGvSRJBTPoJUkqmEEvSVLBDHpJkgpm0EuSVDCD\nXpKkgjW61n1ErAWuAgaA6zJz45T3fwi4AVgN/Fpmvq9uX0mS1FljM/qIGACuAdYCZwDrI+L0Kc0e\nBt4O/N5R9JUkSR00eep+DbAzM3dl5mFgM3Bxe4PMHM3M7cDhufaVJEmdNRn0K4CRtu3drX1N95Uk\nSS1NXqOfPJZ9Fw/0x32Fg4MDDA0t79hu375lx6Ca7jn55GUdx1XimKC/xlV3TE3q9ec3wTH1hxLH\nVEeTQb8HGG7bHqaamTfSd2J8oi8eIRgbG2d09EDHdnv3HjwG1XTP3r0HO46rxDEdadcv6o6pKUND\ny3v6+U1wTP2hxDFBvS8vTQb9dmBVRKwEHgDWAetnaLtoHn0lSdIMGgv6zByLiMuAW6kekbs+M3dE\nxIbW+5si4lTgK8BTgImIuBw4IzMPTte3qVolSSpVo8/RZ+ZWYOuUfZvaXn+LJ56in7WvJEmam364\nrC1Jko6SQS9JUsEMekmSCmbQS5JUMINekqSCGfSSJBXMoJckqWAGvSRJBTPoJUkqmEEvSVLBDHpJ\nkgpm0EuSVLBGf9RGUn86dOgQIyP3d/24+/YtY+/eg1095vDws1myZElXjymVxKCX9F1GRu7n8vdu\nYemJp/S6lFk9tv8hrr7iIk47bVWvS5EWLINe0rSWnngKy05a0esyJM2T1+glSSqYQS9JUsEMekmS\nCmbQS5JUMINekqSCGfSSJBXMoJckqWAGvSRJBTPoJUkqmEEvSVLBDHpJkgpm0EuSVDCDXpKkgjX6\n63URsRa4ChgArsvMjdO0eT9wAfAY8IbMvLu1fxfwCDAOHM7MNU3WKklSiRqb0UfEAHANsBY4A1gf\nEadPafPjwHMzcxXwFuBDbW9PAudk5mpDXpKko9Pkqfs1wM7M3JWZh4HNwMVT2lwE/BFAZt4FPDUi\nvq/t/UUN1idJUvGaDPoVwEjb9u7WvrptJoE7ImJ7RFzSWJWSJBWsyaCfrNlupln7KzJzNdX1+0sj\n4uzulCVJ0vGjyZvx9gDDbdvDVDP22do8q7WPzHyg9fdoRNxMdSlg20wftnigPx4gGBwcYGhoecd2\n+/YtOwbVdM/JJy/rOK4SxwT9Na7jeUxN6vXnN8ExlaPJoN8OrIqIlcADwDpg/ZQ2W4DLgM0R8XLg\nnzPz2xGxFBjIzAMRcQJwPvCu2T5sYnyiL54VHBsbZ3T0QMd2e/cePAbVdM/evQc7jqvEMR1p1y+O\n5zE1ZWhoeU8/vwmOqX/U+fLSWDZm5hhViN8K3AvcmJk7ImJDRGxotbkFuC8idgKbgJ9vdT8V2BYR\n9wB3AZ/JzNuaqlWSpFI1+hx9Zm4Ftk7Zt2nK9mXT9LsPeHGTtUmSdDzoh7PdkiTpKBn0kiQVzKCX\nJKlgBr0kSQUz6CVJKphBL0lSwQx6SZIKZtBLklQwg16SpIIZ9JIkFcyglySpYAa9JEkFM+glSSqY\nQS9JUsEMekmSCmbQS5JUMINekqSCDfa6AEk6Fg4dOsTIyP1dP+6+fcvYu/dgV485PPxslixZ0tVj\n6vhl0Es6LoyM3M/l793C0hNP6XUps3ps/0NcfcVFnHbaql6XokIY9JKOG0tPPIVlJ63odRnSMeU1\nekmSCmbQS5JUMINekqSCGfSSJBXMoJckqWAGvSRJBTPoJUkqmM/RS1Kf6qfV/sAV/3ql0aCPiLXA\nVcAAcF1mbpymzfuBC4DHgDdk5t11+0rS8axfVvuD+iv++eWl+xoL+ogYAK4BzgP2AF+JiC2ZuaOt\nzY8Dz83MVRHxMuBDwMvr9JUklbfaX4lfXnqtyRn9GmBnZu4CiIjNwMVAe1hfBPwRQGbeFRFPjYhT\ngefU6CtJKlBpX156rcmb8VYAI23bu1v76rR5Zo2+kiSpgyZn9JM12y3qxof96/4RntSNA7UZGFzM\n+NhEV4956NBo7baP7X+oq5/dlLnUWeKYjqZ9Lzgmx9RL/n+idxZNTtbN47mJiJcDv5mZa1vb7wQm\n2m+qi4g/BP53Zm5ubX8T+FGqU/ez9pUkSZ01eep+O7AqIlZGxBJgHbBlSpstwOvh8S8G/5yZ367Z\nV5IkddBY0GfmGHAZcCtwL3BjZu6IiA0RsaHV5hbgvojYCWwCfn62vk3VKklSqRo7dS9JknrPJXAl\nSSqYQS9JUsEMekmSCuaP2sygtLX2I+IjwIXAQ5n5gl7X0w0RMQx8FDiFat2GD2fm+3tb1fxFxPcC\nnweeBCwB/iIz39nbquavtbT1dmB3Zr661/V0Q0TsAh4BxoHDmbmmpwV1QUQ8FbgOeB7Vf1c/l5l/\n3duqjl5EBLC5bdcPAL/e7/+vaD12/jpgAvgG8MbM/Nfp2jqjn0bbWvtrgTOA9RFxem+rmrcbqMZT\nksPAL2bm84CXA5cW8M+JzPwOcG5mvhh4IXBuRLyix2V1w+VUT9GUdAfwJHBOZq4uIeRbrgZuyczT\nqf796+snnrKyOjNXAy+h+gG1m3tc1rxExErgEuDM1sRtAHjNTO0N+uk9vk5/Zh6m+jZ4cY9rmpfM\n3Abs63Ud3ZSZ38rMe1qvD1L9D+mZva2qOzLzsdbLJVT/Ee/tYTnzFhHPAn6caqbYldUwF5BixhMR\nJwJnZ+ZHoHrUOTP397isbjoP+IfMHOnYcmF7hGqiszQiBoGlVD8ANy2Dfnp11unXAtL6hrsauKvH\npXRFRCyOiHuAbwN3Zua9va5pnv4AuILqNGNJJoE7ImJ7RFzS62K64DnAaETcEBFfjYhrI2Jpr4vq\notcAf9LrIuYrM/cC7wP+EXiAarG5O2Zqb9BPr6RTi8WLiGXAJ4HLWzP7vpeZE61T988CXhkR5/S4\npKMWET9BdW/I3RQ0+205q3VK+AKqS0dn97qgeRoEzgQ+mJlnAo8Cv9rbkrqjtcrqq4E/63Ut8xUR\npwG/AKykOou5LCJeO1N7g356e4Dhtu1hqlm9FpiI+B7gz4E/zsxP9bqebmudNv1L4Id7Xcs8/Ahw\nUUT8P+ATwKsi4qM9rqkrMvPB1t+jVNd9+/06/W6qmyW/0tr+JFXwl+AC4P+0/ln1ux8GvpSZD7dW\nkr2J6r+zaRn003Ot/T4QEYuA64F7M/OqXtfTLRHx9Nadz0TEk4EfA+7ubVVHLzP/a2YOZ+ZzqE6d\nfi4zX9/ruuYrIpZGxPLW6xOA86nufu5bmfktYCQifrC16zzg73pYUjetp/qiWYJvAi+PiCe3/j94\nHtWNrtMy6KdR4lr7EfEJ4EvAD0bESES8sdc1dcFZVI+XnBsRd7f+lPBkwTOAz7Wu0d8FfDozP9vj\nmrqplEsAsXFpAAABFklEQVRj3wdsa/vn9JnMvK3HNXXD24GPR8TXqO66/589rmfeWl/EzqOa+fa9\nzPwa1aPF24Gvt3Z/eKb2rnUvSVLBnNFLklQwg16SpIIZ9JIkFcyglySpYAa9JEkFM+glSSqYQS9J\nUsEMekmSCmbQS5JUsMFeFyCpf0XET1P9FvZK4H7geZl5RU+LkvQEBr2koxIRzwfupDoz+BHgA8Df\n97QoSd/Fte4lzUtrVr+ipF8QlEpi0Es6KhHxIuAR4FeAG4CvAi/LzC/0tDBJT+Cpe0lH63zgMWAn\nsAZ4LvCnPa1I0ndxRi9JUsF8vE6SpIIZ9JIkFcyglySpYAa9JEkFM+glSSqYQS9JUsEMekmSCmbQ\nS5JUsP8PSe13gWnlsaUAAAAASUVORK5CYII=\n", "text/plain": "<matplotlib.figure.Figure at 0x7faf80e98a50>"}, "metadata": {}}], "metadata": {"collapsed": false, "trusted": true}}, {"source": "We are now ready to estimate $\\pi$ and $\\lambda$ by maximum likelihood.  To do so, we define a class that inherits from `statsmodels`' `GenericLikelihoodModel` as follows.", "cell_type": "markdown", "metadata": {}}, {"execution_count": 10, "cell_type": "code", "source": "class ZeroInflatedPoisson(GenericLikelihoodModel):\n    def __init__(self, endog, exog=None, **kwds):\n        if exog is None:\n            exog = np.zeros_like(endog)\n            \n        super(ZeroInflatedPoisson, self).__init__(endog, exog, **kwds)\n    \n    def nloglikeobs(self, params):\n        pi = params[0]\n        lambda_ = params[1]\n\n        return -np.log(zip_pmf(self.endog, pi=pi, lambda_=lambda_))\n    \n    def fit(self, start_params=None, maxiter=10000, maxfun=5000, **kwds):\n        if start_params is None:\n            lambda_start = self.endog.mean()\n            excess_zeros = (self.endog == 0).mean() - stats.poisson.pmf(0, lambda_start)\n            \n            start_params = np.array([excess_zeros, lambda_start])\n            \n        return super(ZeroInflatedPoisson, self).fit(start_params=start_params,\n                                                    maxiter=maxiter, maxfun=maxfun, **kwds)", "outputs": [], "metadata": {"collapsed": false, "trusted": true}}, {"source": "The key component of this class is the method `nloglikeobs`, which returns the negative log likelihood of each observed value in `endog`.  Secondarily, we must also supply reasonable initial guesses of the parameters in `fit`.  Obtaining the maximum likelihood estimate is now simple.", "cell_type": "markdown", "metadata": {}}, {"execution_count": 11, "cell_type": "code", "source": "model = ZeroInflatedPoisson(x)\nresults = model.fit()", "outputs": [{"output_type": "stream", "name": "stdout", "text": "Optimization terminated successfully.\n         Current function value: 1.586641\n         Iterations: 37\n         Function evaluations: 70\n"}], "metadata": {"collapsed": false, "trusted": true}}, {"source": "We see that we have estimated the parameters fairly well.", "cell_type": "markdown", "metadata": {}}, {"execution_count": 12, "cell_type": "code", "source": "pi_mle, lambda_mle = results.params\n\npi_mle, lambda_mle", "outputs": [{"execution_count": 12, "output_type": "execute_result", "data": {"text/plain": "(0.31542487710071976, 2.0451304204850853)"}, "metadata": {}}], "metadata": {"collapsed": false, "trusted": true}}, {"source": "There are many advantages to buying into the `statsmodels` ecosystem and subclassing `GenericLikelihoodModel`.  The already-written `statsmodels` code handles storing the observations and the interaction with `scipy.optimize` for us.  (It is possible to control the use of `scipy.optimize` through keyword arguments to `fit`.)\n\nWe also gain access to many of `statsmodels`' built in model analysis tools.  For example, we can use bootstrap resampling to estimate the variation in our parameter estimates.", "cell_type": "markdown", "metadata": {}}, {"execution_count": 13, "cell_type": "code", "source": "boot_mean, boot_std, boot_samples = results.bootstrap(nrep=500, store=True)\nboot_pis = boot_samples[:, 0]\nboot_lambdas = boot_samples[:, 1]", "outputs": [{"output_type": "stream", "name": "stdout", "text": "<class '__main__.ZeroInflatedPoisson'>\n"}], "metadata": {"collapsed": false, "trusted": true}}, {"execution_count": 14, "cell_type": "code", "source": "fig, (pi_ax, lambda_ax) = plt.subplots(ncols=2, figsize=(16, 6))\n\nsns.boxplot(boot_pis, ax=pi_ax, names=['$\\pi$'], color=palette[0]);\nsns.boxplot(boot_lambdas, ax=lambda_ax, names=['$\\lambda$'], color=palette[1]);\n\nfig.suptitle('Boostrap distribution of maximum likelihood estimates');", "outputs": [{"output_type": "display_data", "data": {"image/png": "iVBORw0KGgoAAAANSUhEUgAAA6cAAAGSCAYAAAAfNLdHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XucXWWd5/tPWVUqRUKITXmDwihN/wQaMKiBY7g6DCfi\nBUeYgSDSiDIMLQ7tdKPd9EFtOa+2g+Mx0MxBIkg3gmDLbXCaNMo4QriIRMJFib+eCIEQRAtMQkIJ\nVJKaP9Yq3BRVtStJJU921ef9euVVez3rWWv/9qWy67uftZ7VNjAwgCRJkiRJJb2qdAGSJEmSJBlO\nJUmSJEnFGU4lSZIkScUZTiVJkiRJxRlOJUmSJEnFGU4lSZIkScUZTiVpjCJiQ0QsiYj7I+KnEfF/\njfP+zxnP/TW5r+UR8br69p1N+o5aV0T8c0TsFBEzIuKhTazjsMbnMSJOj4iPbco+NlVEfCUifhYR\n87bm/TTc3zciYq9tcV/D3PdLr0lEvCsiLqhvfzEi/nwT9vOK90tEHB4R39sKNY/Lfoe+b5u9zzdh\nvy97z0qSxk9H6QIkqYX0ZeZMgIg4CvgycPg47v+vgL8d2hgRbQCZOZ4Xpn5pX5k5ewvren+9/LrN\nqOMIYC1wd72vSzZjH5vqNGD6OD+fI8rM07bF/TSTmYuBxfXipj72TXm/bC9e9r4dx7pf9p6VJI0f\nw6kkbZ5pwG/hpZB2PjCH6o/4/zcz/2mU9jcB3wGmUv0/fAbwAWCHiFgC/Az4f4DvAz8G3gkcHRF/\nCbwb2AG4NjO/WN//8np/7wN+B5yYmb9sLDYi/gC4Gngz1R/VbQ3r1mXmlC2o63bggHp3HRFxZb38\nc+DkzPxdXeMBmfnbiHgX8BXgFOB0YENEnAR8GjgSWJuZX42IdwBfrx/vL4FTM3N1RPyovv8jgJ2B\nT2TmHUNfoIj4yjDP/U3AFOC+iPhyZv5TQ/8vAm+t/+0O/BfgPcBRwErgg5m5PiLOBT5Y13VXZp4e\nER3AXcDZmXlbRHwZWJ+Z59b1/pfMvC8i1gH/P3A08CvgXGAesBvwZ5n5vYg4BXhnZn66rut/AOdn\n5u1j2X7o89Dw+A4H/jwzP1g3DdTtpwH/DvgIcFz9OrwauAf408zcOGQ/6zJzSr04JSK+C/wx8NPM\nPKnu82+oXuMO4F7gjMx8cZT2OcDXgD7gFa9lvc924O+Aw4DXAP8tMxeM5X2bmR9reJ8fDvwNsArY\nF/gu1Xv108BrgQ9n5iMR8UHgr+vn4hngo0AXL3/Pngn8K3Ax1XuG+nW4KyIOA+Y3PNeHZua6EV4e\nSRIe1itJm2KH+rDepcA3gPPq9o8A+wP7UYWrr0TEG0dpPxH4l3oUdn/g/sz8S+B3mTkzMz9GFR7/\nkOoP8D/OzMeBv87Md9fbHBYRf1zf/wCwOjP3Ay7i938QN/oCcHtm/jFwA7//Q3pwe7agrsZRuKjX\n7Q08C/zpkPt4SWY+RhU+/796/3fU/Qb7XkEV9vYHHqofw+C+2jPzQODPGtp/X0TEsbzyuX9DZn6o\n4fH809DtqILpEcCHgCuBH9TP6++A99d9LsrMWZm5L9V74gOZuZ4qbF8cEUcC/zdVABr62LuA/1m/\nDmuBLwHvpQqHXxqmnvHYfiRtEXEmVdA9pn7s/wF4T/0e2EgVyEarZyZwFrA38LaIeE9EvBa4HPgP\n9XPXAZzRpH0B8IHMfCfwRoYf2f0E1ft8FjALOC0iZgBzaf6+HVr3flQhcy/gY8Ae9X4vpQqpAIsy\n86DMPIAq/H42M5fz8vfsncAFwNfq7Y+r9wHw51ThfiZwMNV7SJI0CsOpJI3d4B+7e1GNyH2rbj8Y\n+HZmDmTmb4DbqEY4Z4/Q/hPg4xHxBWDfUUZTHsvMnzQsHx8RPwXuA/ahCgSDrq5/XgMMdz7cIVRh\ni8y8mWrUaKjNravRiswcPNzxSqrnppm2oQ0RsRMwLTMX1U3/CBza0OX6+ud9wIxh9jnScz+aAWBh\nZm6gGiV+VWbeUq97qOF+3hsRP46IB6mC4T4Amfkw1WP+HvDxOrAO9eKQff6vhvsb7nGM9/aD2oCT\nqd7Hx2VmP/BvqEbDF9cjju+lCqyj+UlmPlkfIn1/3T+ARzNzWd1n8LX7oxHaB/sPjvZfyTDvCaoR\n7JPr2n4MvI7qi5J7Gdv7ttG9mfnrzHwRWAYMPqeNz2NPRHy/fp3/gpf/vjXWdyRwUV3XfwemRsSO\nwJ3A1yLi01SHkW8YQ12SNKkZTiVpM2Tmj4FdIqKbKtQM98c0w7QP1IHrEKpDRf9hlAmAnhu8ERFv\npRqJeW89kvjPVIcgDmek8wlHqhGAzamryX23NSyv5/efOSPVPZqhtb9Q/9zAyKeotI1wezQvAtSH\nsvY3tG8E2iPiNcB/A46tR/++wcsfz75Uwf8NI+x/6D4b72/wcTQ+VwzZ/1i2H4sBqnD7FqCnof0f\n6y9gZmbm2zOz2WjsCw23B1+Loe+/sf5uNGsHOLOhvj0y89ZNeN+OVPfGhuXG5/HvgQvr1/l0qsO4\nR6r3wIa6ejLzucycRzXauwNwZ0TEGOqSpEnNcCpJmyEi3k71f+jTwCKqUc1X1WH1UKrz9YZr/0lE\n7A70ZualwGVUh0YC9NfnLg5nJ6pQ+GxEvIHq/NJGxzf8vGuY7W+nOmyXiHgfMH2Yx7Q5dQ21e0Qc\nVN8+keo5AFgOvKu+fWxD/7VU5wo2asvMZ4FVETE48vox4EdjrAFe+dwfQjUyvCXa+H1QfCYipgD/\nnt+fu/kRqnNgDwP+PiKmbeb9LAfeERFtEdFDdQjreGsDlgD/CbipPm/zfwLH1c8XEfG6+j2xKQaA\nBGZExB512+BrN1L7L+r2t9Xtc0fY9y3Anw6+FyPijyKia5zet8PZCXiyvn1KQ/vQ9+z3gf88uFCf\nK01E7JGZP8/M86lGdw2nktSE4VSSxm7wnNMlVIfP/kl92OgNwIPAA1R/4J+dmb8ZqZ1qht/7I+I+\nqnBzQb3/BcCDEfEtXn7uJZn5AFWY+AVwFa+cNGZ6RDxAdb7cZ4ap/W+AQyPiZ1TnJz7WsG7wfo7Y\n1LqGbA9VAPlURDxMNWnUxQ33f0FE3Es1Mji4zfeAfxcR9zUE0cF1f0J1rugDVOcIjuWczKqIkZ/7\nYfuPsK9XPM7MXEM1Wvoz4F+ovoQYnHDqy8AnM/N/M/K5v6M9dwN17XcAjwIPU70GP92U7Zvc50DD\nz4H6nMm/oBqJ/w31hFf1c/59qvM/x7K/l2TmC8DHge/Wh8SuB77epP0/Av9cH7b+6xEey6VUz8l9\nUV0e52KqUc7Daf6+bVp34/NS3/5iXetioJdXvmeXRMRsqmD6roh4ICJ+Xj8WgLMi4qH6uXwRWDjC\nfUqSam0DA9tkJn1J0lYSEY9Sze7629K1SJIkbS5HTiWp9fktoyRJanmOnEqSJEmSinPkVJIkSZJU\nnOFUkiRJklSc4VSSJEmSVJzhVJIkSZJUnOFUkiRJklSc4VSSJEmSVJzhVJIkSZJUnOFUkiRJklSc\n4VSSJEmSVJzhVJIkSZJUnOFUkiRJklSc4VSSJEmSVJzhVJIkSZJUnOFUkiRJklSc4VSSJEmSVJzh\nVJIkSZJUnOFUkiRJklSc4VSSJEmSVJzhVJIkSZJUXEezDhExB5gPtAOXZua8IeuPAb4EbKz/nZ2Z\nP6zX/RVwUt3+EPDxzHxhXB+BJEmTTET0AFcArwcGgAWZeeGQPh8FPgu0AWuBMzLzwXrdcuBZYAPQ\nn5mztlnxkiSNYNSR04hoBy4C5gB7A3MjYq8h3W7NzP0zcyZwCrCg3nYGcBpwQGbuSxVuTxjX6iVJ\nmpz6gc9k5j7AQcCnhvl8fgQ4NDP3A86j/nyuDQCHZ+ZMg6kkaXvRbOR0FrAsM5cDRMQ1wDHA0sEO\nmflcQ/8pwNP17WepPjy7ImID0AWsHJ+yJUmavDLzKeCp+va6iFgKvJmXfz7f3bDJPcBuQ3bTtrXr\nlCRpUzQ753RXYEXD8hN128tExIfrD8aFwH8GyMzfAl8FHgeeBFZn5q3jUbQkSarURyrNpAqgI/kE\ncHPD8gBwa0QsjojTtmJ5kiSNWbOR04Gx7CQzbwRujIhDgG8BERF7AH8GzADWAN+NiI9m5lUj7Wf9\n+g0DHR3tYypckqQxmNCjgxExBbgWOCsz143Q5wjgVGB2Q/PszPxVRHQDP4iIX2TmouG297NZGt26\ndes444wzAPj617/OjjvuWLgiabs34mdzs3C6EuhpWO6hGj0dVmYuioiOiNgFeBdwV2Y+AxAR1wPv\nAUYMp6tW9TUpR5Kksevunlq6hK0mIjqB64Ar6y+Jh+uzH/ANYE5mrhpsz8xf1T97I+IGqtN4hg2n\nfjZLo+vr62Pjxo0APP30Ovr6NhauSNq+jfbZ3CycLgb2rA8ZehI4Hpjb2KEeIX0kMwci4gCAzHw6\nIhI4NyJ2AJ4HjgR+srkPQpIkVSKiDbgMeDgz54/QZ3fgeuCkzFzW0N4FtGfm2ojYETgK+JttULY0\nIXV1dXHccXNpa2ujq6urdDlSSxs1nGbm+og4E7iFarbdyzJzaUScXq+/BDgWODki+oF11DPyZub9\nEXEFVcDdCNzHy2cKlCRJm2c21aXaHoyIJXXbOcDu8NLn8+eB6cDFEQG/v2TMG4Hr67YO4KrM/P62\nLV+aWGbPPrR0CdKE0DYwMKbTSreJ3t61208xkqSW1909dUKfc7ot+NksSRpPo302N5utV5IkSZKk\nrc5wKkmSJEkqznAqSZIkSSrOcCpJkiRJKs5wKkmSJEkqznAqSZIkSSrOcCpJkiRJKs5wKkmSJEkq\nznAqSZIkSSrOcCpJkiRJKs5wKkmSJEkqznAqSZIkSSrOcCpJkiRJKs5wKkmSJEkqznAqSZIkSSrO\ncCpJkiRJKs5wKkmSJEkqrqN0AZIkSSrnzjtv5447bitdRktbs2Y1ANOm7Vy4ktZ28MGHMXv2oaXL\nUEGOnEqSJElbYM2aNaxZs6Z0GVLLaxsYGChdw0t6e9duP8VIklped/fUttI1tDo/m6Xm5s07D4DP\nfe7cwpVI27/RPpsdOZUkSZIkFWc4lSRJkiQVZziVJEmSJBVnOJUkSZIkFWc4lSRJkiQVZziVWszA\nwADb0yzbkiRJ0ngwnEot5s47b+euuxaVLkOSJEkaVx2lC5A0dn19z3HttVcDMHPmu+jq6ipckSRJ\nkjQ+HDmVWsqI1yyWJEmSWpojp1IL6erq4thjT6Ctrc1RU0mSJE0ojpxKLcgJkSRJkjTROHIqtZC+\nvuf49revAOCd75zl6KkkSZImDEdOpRbyu9/18cILz/PCC8/zu9/1lS5HkiRJGjeOnEotZGDACZGk\nyS4ieoArgNcDA8CCzLxwSJ+PAp+lmkVtLXBGZj5Yr5sDzAfagUszc942LF+SpBE5ciq1kK6uLl7z\nmtfymte8lh128JBeaZLqBz6TmfsABwGfioi9hvR5BDg0M/cDzgMWAEREO3ARMAfYG5g7zLaSJBXh\nyKnUQrq6unj3uw90tl5pEsvMp4Cn6tvrImIp8GZgaUOfuxs2uQfYrb49C1iWmcsBIuIa4JjGbSVJ\nKsVwKrWQvr7nePDBJfXtPgOqNMlFxAxgJlUAHckngJvr27sCKxrWPQEcuFWKkyRpExlOpZbiOaeS\nKhExBbgWOCsz143Q5wjgVGB23bTJ16GaPr2Ljo72za5Tmgw6O6vfke7uqYUrkVqb4VTb1J133s4d\nd9xWuowW10ZbG/z933+1dCEt7eCDD2P27ENLlyFtlojoBK4DrszMG0fosx/wDWBOZq6qm1cCPQ3d\neqhGT0e0apUzg0vN9PdvAKC3d23hSqTt32hf4jQNp81m9YuIY4AvARvrf2dn5g/rdTsDlwL7UH1b\ne2pm/njzHoYkgBdffLF0CZIKiog24DLg4cycP0Kf3YHrgZMyc1nDqsXAnvXhwE8CxwNzt27FkiSN\nzajhtGFWvyOpvm29NyJuyszGiRNuzcz/XvffF7gB+MN63QXAzZl5XER0ADuO9wNQa5k9+1BHq7bQ\nvHnnAfC5z51buBJJhcwGTgIejIgldds5wO4AmXkJ8HlgOnBxRAD0Z+aszFwfEWcCt1B96XzZkM90\nSZKKaTZy2nRWv8x8rqH/FODpuu804JDM/JO633pgzbhVLknSJJSZd9DkUnCZ+UngkyOsWwgs3Aql\nSZK0RZqF0zHN6hcRHwa+DLwJOKpufivQGxGXA/sDP6WatMGTVyRJkiRJL9MsnI5pVr96MoYbI+IQ\n4FtA1Ps+ADgzM++NiPnAX1IdajQsZwSUmnNGQEmSJE1EzcLpJs3ql5mLIqIjIv6g7vdEZt5br76W\nKpyOyBkBpeacEVAaO7/EkSSpdYx6zgoNs/pFxKupZvW7qbFDROxRzxxIRBwAkJnPZOZTwIqI+KO6\n65HAz8e1ekmSJEnShDDqyOlIs/pFxOn1+kuAY4GTI6IfWAec0LCLTwNX1cH2l8DHt8JjkCRJkiS1\nuKbXOR1uVr86lA7ePh84f4RtHwDevYU1SpIkSZImuGaH9UqSJEmStNUZTiVJkiRJxRlOJUmSJEnF\nGU4lSZIkScUZTiVJkiRJxRlOJUmSJEnFGU4lSZIkScUZTiVJkiRJxRlOJUmSJEnFGU4lSZIkScUZ\nTiVJkiRJxRlOJUmSJEnFGU4lSZIkScUZTiVJkiRJxRlOJUmSJEnFGU4lSZIkScUZTiVJkiRJxRlO\nJUmSJEnFGU4lSZIkScUZTiVJkiRJxRlOJUmSJEnFGU4lSZIkScUZTiVJkiRJxRlOJUmSJEnFGU4l\nSZIkScUZTiVJkiRJxXWULkCSJG2aiOgBrgBeDwwACzLzwiF93g5cDswE/jozv9qwbjnwLLAB6M/M\nWdumckmSRmY4lSSp9fQDn8nM+yNiCvDTiPhBZi5t6PMM8Gngw8NsPwAcnpm/3Qa1SpI0Jh7WK0lS\ni8nMpzLz/vr2OmAp8OYhfXozczFVkB1O29atUpKkTWM4lSSphUXEDKpDd+/ZhM0GgFsjYnFEnLZV\nCpMkaRMZTiVJalH1Ib3XAmfVI6hjNTszZwLvAz4VEYdslQIlSdoEnnMqSVILiohO4Drgysy8cVO2\nzcxf1T97I+IGYBawaLi+06d30dHRvqXlShNaZ2f1O9LdPbVwJVJrM5xKktRiIqINuAx4ODPnN+n+\nsnNLI6ILaM/MtRGxI3AU8DcjbbxqVd+WlitNeP39GwDo7V1buBJp+zfalziGU0mSWs9s4CTgwYhY\nUredA+wOkJmXRMQbgXuBnYCNEXEWsDfV5Weujwio/g64KjO/v43rlyTpFQynkiS1mMy8gybzRmTm\nU0DPMKvWAe/YGnVJkrQlnBBJkiRJklSc4VSSJEmSVJzhVJIkSZJUXNNzTiNiDjAfaAcuzcx5Q9Yf\nA3wJ2Fj/Ozszf9iwvh1YDDyRmR8cx9olSZIkSRPEqCOndbC8CJhDNcPf3IjYa0i3WzNz//pi3qcA\nC4asPwt4GBgYl4olSZIkSRNOs8N6ZwHLMnN5ZvYD1wDHNHbIzOcaFqcATw8uRMRuwNHApQy5zpok\nSZIkSYOaHda7K7CiYfkJ4MChnSLiw8CXgTdRXcx70NeAs6musSZJkiRJ0rCajZyO6VDczLwxM/cC\nPgh8KyLaIuIDwG8ycwmOmkqSJEmSRtFs5HQlL7+Adw/V6OmwMnNRRHQAfwC8B/hQRBwNvBbYKSKu\nyMyTR9p++vQuOjrax1y8NBl1dla/I93dUwtXIkmSJI2fZuF0MbBnRMwAngSOB+Y2doiIPYBHMnMg\nIg4AyMyngXPqf0TEYcBfjBZMAVat6tucxyBNKv39GwDo7V1buBJp++eXOJIktY5Rw2lmro+IM4Fb\nqC4lc1lmLo2I0+v1lwDHAidHRD+wDjhhhN05W68kSZIkaVhNr3OamQuBhUPaLmm4fT5wfpN93Abc\ntpk1SpIkSZImuGYTIkmSJEmStNUZTiVJkiRJxRlOJUmSJEnFGU4lSZIkScUZTiVJkiRJxRlOJUmS\nJEnFGU4lSZIkScUZTiVJkiRJxRlOJUmSJEnFGU4lSZIkScUZTiVJkiRJxRlOJUmSJEnFGU4lSZIk\nScUZTiVJkiRJxRlOJUmSJEnFGU4lSZIkScUZTiVJkiRJxRlOJUmSJEnFGU4lSZIkScUZTiVJkiRJ\nxRlOJUmSJEnFdZQuQJIkjV1E9ABXAK8HBoAFmXnhkD5vBy4HZgJ/nZlfbVg3B5gPtAOXZua8bVW7\nJEmjceRUkqTW0g98JjP3AQ4CPhURew3p8wzwaeC/NjZGRDtwETAH2BuYO8y2kiQVYTiVJKmFZOZT\nmXl/fXsdsBR485A+vZm5mCrINpoFLMvM5ZnZD1wDHLMNypYkqSnDqSRJLSoiZlAdunvPGDfZFVjR\nsPxE3SZJUnGGU0mSWlBETAGuBc6qR1DHYmArliRJ0hZxQiRJklpMRHQC1wFXZuaNm7DpSqCnYbmH\navR0RNOnd9HR0b7pRUqTSGdn9TvS3T21cCVSazOcSpLUQiKiDbgMeDgz5zfp3jZkeTGwZ3048JPA\n8cDc0XawalXfZlYqTR79/RsA6O1dW7gSafs32pc4htNN8O1vX8GKFY+VLkOT3OOPV+/BefPOK1yJ\nBD09b+HEE08uXcZkMxs4CXgwIpbUbecAuwNk5iUR8UbgXmAnYGNEnAXsnZnrIuJM4BaqS8lclplL\nt/kjkCRpGIbTTbBixWPk/15G+2t3Ll2KJrGNG6pDh5ateLpwJZrsNjy/unQJk1Jm3kGTOSMy8yle\nfvhu47qFwMKtUJokSVvEcLqJ2l+7M1NmHF66DEkqbt3yH5UuQfKoJm0XPKpJ25NWPqrJcCpJklrW\nihWP8a+/TDqmvbp0KZrENrZX55w+8vSjhSvRZLd+zYulS9gihlNJktTSOqa9mp0P8XKtkrR60crS\nJWwRr3MqSZIkSSrOcCpJkiRJKs5wKkmSJEkqznAqSZIkSSrOcCpJkiRJKs5wKkmSJEkqbkyXkomI\nOcB8oB24NDPnDVl/DPAlYGP97+zM/GFE9ABXAK8HBoAFmXnhONYvSZIkSZoAmo6cRkQ7cBEwB9gb\nmBsRew3pdmtm7p+ZM4FTgAV1ez/wmczcBzgI+NQw20qSJEmSJrmxHNY7C1iWmcszsx+4BjimsUNm\nPtewOAV4um5/KjPvr2+vA5YCbx6PwiVJkiRJE8dYDuvdFVjRsPwEcODQThHxYeDLwJuAo4ZZPwOY\nCdyzOYVKkiRJkiausYTTgbHsKDNvBG6MiEOAbwExuC4ipgDXAmfVI6jDmj69i46O9rHcXRGdndtv\nbZJUQmdnO93dU0uXIUmSJoCxhNOVQE/Dcg/V6OmwMnNRRHRExB9k5jMR0QlcB1xZB9gRrVrVN5aa\ni+nv31C6BEnarvT3b6C3d23pMkZkcJYkqXWMJZwuBvasD8t9EjgemNvYISL2AB7JzIGIOACgDqZt\nwGXAw5k5f1wrlyRJkiRNGE3DaWauj4gzgVuoLiVzWWYujYjT6/WXAMcCJ0dEP7AOOKHefDZwEvBg\nRCyp2/4qM/9lnB+HJEmSJKmFjek6p5m5EFg4pO2ShtvnA+cPs90djG1GYEmSJEnSJGZwlCRJkiQV\nZziVJEmSJBVnOJUkSZIkFWc4lSRJkiQVZziVJEmSJBVnOJUkSZIkFWc4lSRJkiQVZziVJEmSJBVn\nOJUkSZIkFddRuoBWsmbNajY8v5p1y39UuhRJKm7D86tZs8aPEUmSND4cOZUkSZIkFedX3ptg2rSd\n6X12PVNmHF66FEkqbt3yHzFt2s6ly5AkSROEI6eSJEmSpOIMp5IkSZKk4gynkiRJkqTiPOdUkqQW\nExE9wBXA64EBYEFmXjhMvwuB9wF9wCmZuaRuXw48C2wA+jNz1rapXJKkkTlyKklS6+kHPpOZ+wAH\nAZ+KiL0aO0TE0cAfZuaewH8ELm5YPQAcnpkzDaaSpO2F4VSSpBaTmU9l5v317XXAUuDNQ7p9CPjH\nus89wM4R8YaG9W3bolZJksbKcCpJUguLiBnATOCeIat2BVY0LD9Rt0E1cnprRCyOiNO2epGSJI2B\n55xKktSiImIKcC1wVj2COtRIo6MHZ+aTEdEN/CAifpGZi4brOH16Fx0d7eNU8fjr7Nx+a5OkEjo7\n2+nunlq6jM1iOJUkqQVFRCdwHXBlZt44TJeVQE/D8m51G5n5ZP2zNyJuAGYBw4bTVav6xrPscdff\nv6F0CZK0Xenv30Bv79rSZYxotODsYb2SJLWYiGgDLgMezsz5I3S7CTi57n8QsDozfx0RXRExtW7f\nETgKeGgblC1J0qgcOZUkqfXMBk4CHoyIJXXbOcDuAJl5SWbeHBFHR8Qy4Dng43W/NwLXRwRUfwdc\nlZnf36bVj6M1a1azfs0LrF60snQpklTc+jUvsKZzdekyNpvhVJKkFpOZdzCGo58y88xh2h4B3rE1\n6pIkaUsYTiVJUsuaNm1nnulfxc6H7Nq8syRNcKsXrWTatJ1Ll7HZPOdUkiRJklSc4VSSJEmSVJzh\nVJIkSZJUnOFUkiRJklSc4VSSJEmSVJzhVJIkSZJUnOFUkiRJklSc4VSSJEmSVJzhVJIkSZJUnOFU\nkiRJklSc4VSSJEmSVJzhVJIkSZJUXEfpAlrNhudXs275j0qXoUls4/rnAXhVx2sLV6LJbsPzq4Fd\nSpchSZImCMPpJujpeUvpEiQef/wxAHbvMRSotF38f1GSJI2bpuE0IuYA84F24NLMnDdk/THAl4CN\n9b+zM/OHY9m21Zx44smlS5CYN+88AD73uXMLVyJJkiSNn1HPOY2IduAiYA6wNzA3IvYa0u3WzNw/\nM2cCpwALNmFbSZIkSZKaTog0C1iWmcszsx+4BjimsUNmPtewOAV4eqzbSpIkSZIEzQ/r3RVY0bD8\nBHDg0E4R8WHgy8CbgKM2ZVtJkiRJkpqF04Gx7CQzbwRujIhDgG9FxNs3p5jp07vo6GjfnE2lSaOz\ns/od6e6eWrgSSZIkafw0C6crgZ6G5R6qEdBhZeaiiOgAXlf3G/O2AKtW9TUpR1J//wYAenvXFq5E\n2v75JY6In3+8AAAQiklEQVQkSa2jWThdDOwZETOAJ4HjgbmNHSJiD+CRzByIiAMAMvOZiFjTbFtJ\nkiRJkqDJhEiZuR44E7gFeBj4TmYujYjTI+L0utuxwEMRsQS4ADhhtG23zsOQJEmSJLWyptc5zcyF\nwMIhbZc03D4fOH+s20qSJEmSNFSzS8lIkiRJkrTVGU4lSZIkScUZTiVJkiRJxRlOJUmSJEnFGU4l\nSZIkScUZTiVJkiRJxRlOJUmSJEnFGU4lSZIkScUZTiVJkiRJxXWULkCSJI1dRPQAVwCvBwaABZl5\n4TD9LgTeB/QBp2Tmkrp9DjAfaAcuzcx526p2SZJGYziVJKm19AOfycz7I2IK8NOI+EFmLh3sEBFH\nA3+YmXtGxIHAxcBBEdEOXAQcCawE7o2Imxq3bUXr17zI6kUrS5ehSWzjCxsAeNVr2gtXoslu/ZoX\nYZfSVWw+w6kkSS0kM58Cnqpvr4uIpcCbgcaA+SHgH+s+90TEzhHxRuCtwLLMXA4QEdcAxwzZtqX0\n9LyldAkSjz/+GAC77+L7UYXt0tr/LxpOJUlqURExA5gJ3DNk1a7AioblJ+q2Nw/TfuBWLHGrO/HE\nk0uXIDFv3nkAfO5z5xauRGpthlNJklpQfUjvtcBZmblumC5t43E/06d30dHhoYrSaDo7q9+R7u6p\nhSuRWpvhVJKkFhMRncB1wJWZeeMwXVYCPQ3Lu1GNknYOae+p20e0alXflhUrTQL9/dU5p729awtX\nIm3/RvsSx0vJSJLUQiKiDbgMeDgz54/Q7Sbg5Lr/QcDqzPw1sBjYMyJmRMSrgePrvpIkFefIqSRJ\nrWU2cBLwYEQsqdvOAXYHyMxLMvPmiDg6IpYBzwEfr9etj4gzgVuoLiVzWavP1CtJmjgMp5IktZDM\nvIMxHPmUmWeO0L4QWDjedUmStKU8rFeSJEmSVJzhVJIkSZJUnOFUkiRJklSc4VSSJEmSVJzhVJIk\nSZJUnOFUkiRJklSc4VSSJEmSVJzhVJIkSZJUnOFUkiRJklSc4VSSJEmSVJzhVJIkSZJUnOFUkiRJ\nklSc4VSSJEmSVJzhVJIkSZJUnOFUkiRJklSc4VSSJEmSVJzhVJIkSZJUnOFUkiRJklSc4VSSJEmS\nVJzhVJIkSZJUnOFUkiRJklRcR7MOETEHmA+0A5dm5rwh6z8KfBZoA9YCZ2Tmg/W6vwJOAjYCDwEf\nz8wXxvURSJIkSZJa3qgjpxHRDlwEzAH2BuZGxF5Duj0CHJqZ+wHnAQvqbWcApwEHZOa+VOH2hHGt\nXpIkSZI0ITQbOZ0FLMvM5QARcQ1wDLB0sENm3t3Q/x5gt/r2s0A/0BURG4AuYOX4lC1JkiRJmkia\nnXO6K7CiYfmJum0knwBuBsjM3wJfBR4HngRWZ+atm1+qJEmSJGmiajZyOjDWHUXEEcCpwOx6eQ/g\nz4AZwBrguxHx0cy8aqR9TJ/eRUdH+1jvUpqUOjur35Hu7qmFK5EkSZLGT7NwuhLoaVjuoRo9fZmI\n2A/4BjAnM1fVze8C7srMZ+o+1wPvAUYMp6tW9Y29cmmS6u/fAEBv79rClUjbP7/EkSSpdTQLp4uB\nPevJjZ4EjgfmNnaIiN2B64GTMnNZw6pfAOdGxA7A88CRwE/GqW5JkiRJ0gQy6jmnmbkeOBO4BXgY\n+E5mLo2I0yPi9Lrb54HpwMURsSQiflJv+wBwBVXAfbDuu2ArPAZJkiRJUotrep3TzFwILBzSdknD\n7U8Cnxxh2/OB87ewRkmSJEnSBNdstl5JkiRJkrY6w6kkSZIkqTjDqSRJkiSpuKbnnEqSpO1LRHwT\neD/wm8zcd5j104FvAm+jmjH/1Mz8eb1uOfAssAHoz8xZ26hsSZJG5cipJEmt53JgzijrzwHuy8z9\ngZOBCxrWDQCHZ+ZMg6kkaXtiOJUkqcVk5iJg1Shd9gL+V903gRkR0d2wvm0rlidJ0mYxnEqSNPE8\nAHwEICJmAW8BdqvXDQC3RsTiiDitUH2SJL2C55xKkjTx/B1wQUQsAR4CllCdYwpwcGY+WY+k/iAi\nflGPxA5r+vQuOjrat37FUgvr7Kx+R7q7pxauRGpthlNJkiaYzFwLnDq4HBGPAo/U656sf/ZGxA3A\nLGDEcLpqVd/WLVaaAPr7q+9+envXFq5E2v6N9iWOh/VKkjTBRMS0iHh1ffs04LbMXBcRXRExtW7f\nETiKamRVkqTiHDmVJKnFRMTVwGHALhGxAvgC0AmQmZcAewP/EBEDwM+AT9SbvgG4ISKg+hvgqsz8\n/jYuX5KkYRlOJUlqMZk5t8n6u4EYpv1R4B1bqy5JkraEh/VKkiRJkooznEqSJEmSijOcSpIkSZKK\nM5xKkiRJkooznEqSJEmSijOcSpIkSZKKM5xKkiRJkooznEqSJEmSijOcSpIkSZKKM5xKkiRJkooz\nnEqSJEmSijOcSpIkSZKKM5xKkiRJkooznEqSJEmSiusoXYAmlzvvvJ077ritdBkt7fHHHwNg3rzz\nClfS2g4++DBmzz60dBmSJEmqGU6lFjNt2rTSJUiSJEnjznCqbWr27EMdrZIkSZL0Cp5zKkmSJEkq\nznAqSZIkSSrOcCpJkiRJKs5wKkmSJEkqznAqSZIkSSrOcCpJkiRJKs5wKkmSJEkqznAqSZIkSSrO\ncCpJkiRJKs5wKkmSJEkqrqNZh4iYA8wH2oFLM3PekPUfBT4LtAFrgTMy88F63c7ApcA+wABwamb+\neFwfgSRJkiSp5Y06choR7cBFwBxgb2BuROw1pNsjwKGZuR9wHrCgYd0FwM2ZuRewH7B0vAqXJEmS\nJE0czUZOZwHLMnM5QERcAxxDQ8jMzLsb+t8D7Fb3nQYckpl/UvdbD6wZt8qlSWpgYACAtra2wpVI\nkiRJ46dZON0VWNGw/ARw4Cj9PwHcXN9+K9AbEZcD+wM/Bc7KzL7NrFUScOedt9PW1sbs2YeWLkWS\nJEkaN83C6cBYdxQRRwCnArMb9n0AcGZm3hsR84G/BD4/0j6mT++io6N9rHcpTTrr1q3j+uu/A8C/\n/beHs+OOOxauSJIkSRofzcLpSqCnYbmHavT0ZSJiP+AbwJzMXFU3PwE8kZn31svXUoXTEa1a5aCq\nNJq+vj42btwIwNNPr6Ovb2PhiqTtW3f31NIlSJKkMWoWThcDe0bEDOBJ4HhgbmOHiNgduB44KTOX\nDbZn5lMRsSIi/igz/xU4Evj5eBYvTTZdXV0ce+wJtLW10dXVVbocSYVExDeB9wO/ycx9h1k/Hfgm\n8DbgearZ8n9erxt1Fn5JkkoZdbbeehKjM4FbgIeB72Tm0og4PSJOr7t9HpgOXBwRSyLiJw27+DRw\nVUQ8QDVb79+O+yOQJGnyuZxqJv2RnAPcl5n7AydTzZ4/1ln4JUkqoul1TjNzIbBwSNslDbc/CXxy\nhG0fAN69hTVKqvX1Pcd1110DwAEHvNvRU2mSysxF9VFNI9kL+Lu6b0bEjIh4PbAHTWbhlySplFFH\nTiVtb7x8jKQxeQD4CEBEzALeQnWpt+Fm4d91m1cnSdIwmo6cStp+dHV1cdxxcz3nVFIzfwdcEBFL\ngIeAJcAGNmEW/kHOpC8119lZ/Y44CZu0ZQynUovx+qaSmsnMtVSXdwMgIh4FfgnswBhm4W/kTPpS\nc/39GwDo7V1buBJp+zfalzge1iu1mLa2NtraPLxX0sgiYlpEvLq+fRpwW2auo2EW/nr98cBNBUuV\nJOkljpxKktRiIuJq4DBgl4hYAXwB6ISXJi3cG/iHiBgAfgZ8ol63PiIGZ+FvBy7LTCdDkiRtFwyn\nkiS1mMyc22T93UCMsO4Vs/BLkrQ98LBeSZIkSVJxhlNJkiRJUnGGU0mSJElScYZTSZIkSVJxbQMD\nm3w97q2mt3ft9lOMJKnldXdP9bpLW8jP5onvzjtv5447bitdRkt7/PHHANh997cUrqS1HXzwYV7P\nfRIY7bPZ2XolSZKkLTBt2rTSJUgTgiOnkqQJy5HTLednsyRpPI322ew5p5IkSZKk4gynkiRJkqTi\nDKeSJEmSpOIMp5IkSZKk4gynkiRJkqTiDKeSJEmSpOIMp5IkSZKk4gynkiRJkqTiDKeSJEmSpOIM\np5IkSZKk4gynkiRJkqTiDKeSJEmSpOIMp5IkSZKk4gynkiRJkqTiDKeSJEmSpOIMp5IkSZKk4gyn\nkiRJkqTiDKeSJEmSpOIMp5IkSZKk4gynkiRJkqTiDKeSJEnSFhgYGGBgYKB0GVLLM5xKkiRJW+DO\nO2/nrrsWlS5DankdpQuQJEmSWlVf33Nce+3VAMyc+S66uroKVyS1LkdOJUmSpM3WVroAacJw5FSS\nJEnaTF1dXRx33Fza2tocNZW2kOFUkiRJ2gKzZx9augRpQmgaTiNiDjAfaAcuzcx5Q9Z/FPgs1TEN\na4EzMvPBhvXtwGLgicz84DjWLknSpBQR3wTeD/wmM/cdZv0uwJXAG6k+6/9rZv5DvW458CywAejP\nzFnbpmpp4mpr89BeaTyMes5pHSwvAuYAewNzI2KvId0eAQ7NzP2A84AFQ9afBTwMOL+2JEnj43Kq\nz+aRnAksycx3AIcDX42IwS+kB4DDM3OmwVSStD1pNiHSLGBZZi7PzH7gGuCYxg6ZeXdmrqkX7wF2\nG1wXEbsBRwOX4tnikiSNi8xcBKwapcuvgJ3q2zsBz2Tm+ob1fiZLkrY7zcLprsCKhuUn6raRfAK4\nuWH5a8DZwMbNqk6SJG2ObwD7RMSTwANURzENGgBujYjFEXFakeokSRpGs3NOx3wobkQcAZwKzK6X\nP0B1LsySiDh8LPvo7p7qN7mSJG25c4D7M/PwiNgD+EFE7J+Za4HZmfmriOiu239Rj8QOy89mSdK2\n0mzkdCXQ07DcQzV6+jIRsR/Vt7QfyszBw4zeA3woIh4FrgbeGxFXbHnJkiSpifcA3wXIzF8CjwJR\nL/+q/tkL3EB1Co8kScU1C6eLgT0jYkZEvBo4HripsUNE7A5cD5yUmcsG2zPznMzsycy3AicAP8zM\nk8e3fEmSNIxfAEcCRMQbqILpIxHRFRFT6/YdgaOAh4pVKUlSg1EP683M9RFxJnAL1aVkLsvMpRFx\ner3+EuDzwHTg4oiAkaeld7ZeSZLGQURcDRwG7BIRK4AvAJ3w0mfz3wKXR8QDVF9EfzYzfxsRbwOu\nrz+vO4CrMvP7JR6DJElDtQ0MmBklSZIkSWU1O6xXkiRJkqStznAqSZIkSSrOcCpJkiRJKs5wKkmS\nJEkqznAqSZIkSSrOcCpJkiRtpoiYHhHfjojXla5FanVeSkZqERGxK/AVYE9gA/AMcFN9TUNJklRI\nRHwSeFVmLihdi9TKOkoXIGnM3pKZJ0bEicBAZl5duiBJkgTA94BvAoZTaQt4WK/UIjLzrogIYA3Q\nXboeSZJUycxfAztGxE6la5FameFUai0fBe4C9ogIj3yQJGk7EBGvBdYB7y9di9TKDKdSa+nJzFXA\nb4A9ShcjSdJkFxHtwBeBc4EPl61Gam1OiCRJkiRtpoiYD3wjM38eEXcDh2Xmi6XrklqRI6eSJEnS\nZoiI44CfZubP66b/ARxdsCSppTlyKkmSJEkqzpFTSZIkSVJxhlNJkiRJUnGGU0mSJElScYZTSZIk\nSVJxhlNJkiRJUnGGU0mSJElScYZTSZIkSVJx/wcVS32yEo37SgAAAABJRU5ErkJggg==\n", "text/plain": "<matplotlib.figure.Figure at 0x7faf758dbb50>"}, "metadata": {}}], "metadata": {"collapsed": false, "trusted": true}}, {"source": "The next time you are fitting a model using maximum likelihood, try integrating with `statsmodels` to take advantage of the significant amount of work that has gone into its ecosystem.", "cell_type": "markdown", "metadata": {}}, {"execution_count": null, "cell_type": "code", "source": "", "outputs": [], "metadata": {"collapsed": true, "trusted": true}}], "nbformat": 4, "metadata": {"kernelspec": {"display_name": "Python 2", "name": "python2", "language": "python"}, "language_info": {"mimetype": "text/x-python", "nbconvert_exporter": "python", "version": "2.7.6", "name": "python", "file_extension": ".py", "pygments_lexer": "ipython2", "codemirror_mode": {"version": 2, "name": "ipython"}}}}