Likelihood with autodiff.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Fitting with autograd\n",
    "\n",
    "Let's import some libraries and define a helper function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "import autograd.numpy as np\n",
    "from autograd.scipy.stats import norm\n",
    "from autograd import grad, hessian\n",
    "from scipy.optimize import minimize\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# For mapping values in (0, 1) to (-inf, +inf)\n",
    "def expit(x):\n",
    "    return 1 / (1 + np.exp(-x))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's generate a random dataset:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYEAAAEACAYAAABVtcpZAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl4FFXa9/HvDWFfEgiQAAEhAWQRJwSMqCBBWWRRHB19\nfdARlHFEZBl1HEFBwQ13RBSQeVBcUFDgQUAHQozBERUU0bBFEGUxkgACIUEEQu73jzQYkSV0J326\n0/fnurioPn266pet76pTp6pFVTHGGBOayrkOYIwxxh0rAsYYE8KsCBhjTAizImCMMSHMioAxxoQw\nKwLGGBPCzlgERGS6iGSLSPpJnrtHRApEpHaRtlEisklENohIjyLtCSKSLiIbReT5kvsSjDHGeKs4\nRwKvAj1PbBSRGKA7sLVIWyvgeqAV0AuYLCLieXoKMEhVWwAtROQP6zTGGONfZywCqvoJsPckT00A\n7j2hrR8wS1XzVXULsAlIFJFooIaqfuHp9zpwtdepjTHGlAivzgmIyFXAdlVdc8JTDYHtRR5netoa\nAj8Waf/R02aMMcahsLN9gYhUAe6ncCjIGGNMEDvrIgDEAU2Abzzj/THAVyKSSOGef+MifWM8bZlA\no5O0n5SI2A2NjDHGC6oqZ+71m+IOB4nnH6q6VlWjVTVWVZtSOLTTTlV3AguA/yciFUWkKdAMWKmq\nWUCOiCR6CsfNwHtn+EIC6t9DDz3kPINlKlu5LJNlKul/3ijOFNG3gE8pnNGzTURuOfH9mt8KxHrg\nHWA98AEwRH9LdicwHdgIbFLVxV4lNsYYU2LOOBykqv3P8HzsCY/HA+NP0m8V0PZsAxpjjCk9dsVw\nMSUlJbmO8AeWqfgCMZdlKh7LVLrE23Gk0iQiGoi5jDEmkIkIWkonho0xxpRBVgSMMSaEWREwxpgQ\nZkXAGGNCmBUBY4wJYVYEjDEmhFkRMMaYEGZFwBhjQpgVAWOMCWFWBIwxJoRZETDGmBBmRcAYY0KY\nFQFjjAlhVgSM8ZOsrCxuvvlmhg0bxuHDh13HMQawImCM33zxxResXbuWN954g927d7uOYwxgRcCU\ngkOHDjF16lRmzZrl9eeellUNGzakatWqrmMYc9wZP17SmLO1ZMkSnnjiCXbu3EmXLl2oX7++60jG\nmFOwIwFTKs4//3zCw8NdxzDGnIEVAWOMCWFWBIwxJoRZETA++eabb5g3bx65ubmuoxhjvHDGIiAi\n00UkW0TSi7Q9JSIbRORrEZkrIjWLPDdKRDZ5nu9RpD1BRNJFZKOIPF/yX4pxoW/fvtx999288MIL\nrqMYY7xQnCOBV4GeJ7QlA21UNR7YBIwCEJHWwPVAK6AXMFlExPOaKcAgVW0BtBCRE9dpglB+fj49\nevTgyJEjzJkzh+nTp//uQqgxY8aQkZHhMKEx5nTOWARU9RNg7wltKapa4Hn4ORDjWb4KmKWq+aq6\nhcICkSgi0UANVf3C0+914OoSyG8CxLGrYR955BFGjx6NiPD666+zfft2Zs+e7TpeQBER7rnnHtLT\n08/c2ZhSVhLnBG4FPvAsNwS2F3ku09PWEPixSPuPnjZTRqgqERERLF68mH/+8588+eSTdO/enY4d\nOwKQmprK/PnzOXr0qOOk7s2ePRsRoUuXLiQkJDB+/HgOHDjgOpYJUT5dLCYiDwBHVPXtEspz3Nix\nY48vJyUlkZSUVNKbMF7Kz8/nwIEDVKtW7Q/PtWzZkpYtW/6uLTs7mz59+hAREUGlSpXo1auXv6IG\npE6dOnHhhReyZs0a3n//fV555RXOO+88rrzyStfRTJBJS0sjLS3Np3V4XQREZCDQG7isSHMm0KjI\n4xhP26naT6loETCBZfDgwbz22mv079//+J796fZkv/nmG2rVqkVCQoIdCXhUqFCBhIQEEhISWLly\npes4JkiduIM8bty4s15HcYeDxPOv8IHIFcC9wFWqeqhIvwXADSJSUUSaAs2AlaqaBeSISKLnRPHN\nwHtnndYEhKysLB599FE++ugjoqOjufjii0lNTeXCCy/8Q9/rr7+edu3aefXLGSpEhOeff5758+e7\njmJC0BmPBETkLSAJiBSRbcBDwP1ARWCpZ/LP56o6RFXXi8g7wHrgCDBEf7uD2J3ADKAy8IGqLi7h\nr8X4UZs2bdi2bdvxxwMHDjxlvxdffBGA994LzbqfnJzM/v37CQs7+Z/bM888w0svvcT06dO5+mqb\nL2H864xFQFX7n6T51dP0Hw+MP0n7KqDtWaUzJogtX76cf/7zn3z++efUr1+fq6666qT9WrRoQbdu\n3Zg2bZqfExpjVwybYkpNTaVLly7cc889rqMEhcmTJzNixAji4uJYv3497du3p6Cg4MwvNMbPrAiY\nYlm2bBk1atRg1qxZrqMEhUceeYTrrruO5557jlatWp3Va5OTk+nQoQODBw8upXTG/MaKgCm2Bg0a\nuI4QVG6++Wbq1asHQJMmTfj3v/9NkyZNTtl/3759bNq0iU8//ZRGjRqxYMECPyU1ocyKgDF+MGnS\nJFSVSZMmnfT51q1bk5+fz5/+9Cd2795N3bp1/ZzQhCorAsYEgLi4OD777DMiIiLsQ+iNX1kRMH61\nevVqqlevTkxMTJm7/bSqkpOTw+OPP17mvjZTdlkRMH61fft2unTpQn5+Pnl5ea7jlKihQ4cSERHB\nK6+8wuOPP05UVJRP69u/fz+jR48mJyenhBIa80dWBIxfrVu3jvLly/PbHcbLjm3btrFgwQK+++47\nhg8fTrly3v951a9fn4kTJzJ37lw+/vjjEkxpzO9ZETB+c+eddxIfH8+QIUNcRwlYVatWZf78+VSr\nVo1BgwbRrFkz15FMGefTXUSNORu9evUK+TuInslHH31EZmYm7dq1cx3FhAg7EjDFVrlyZXJyckhN\nTSUiIsJ1nDKpUaNGdOzYkUqVKrmOYkKEHQmYYouMjGTz5s0cOnSIxo0b+7y+f//73/z9738nOjq6\nBNK5sWPHDt566y2+/fZb0tPTy+S5DlO22ZGAOStRUVElUgAeffRR3n//fV5++eUSSOXOyy+/zOzZ\ns6lZsyYPPvggPXr0cB3JmLNiRwLGiUGDBv3uVtTBrHfv3qX6IUiqyjfffEPVqlVp3rx5qW3HhCYr\nAsYEsNq1a9OvXz8AqlSpwrZt26hTp47jVKYsseEgYwLYyy+/zNatWzl48CD169enV69epKamuo5l\nyhArAsYEsMqVK9O4cWMqV65Mamoqbdu2tYvHTImyImCcW7RoETNnziQ/P991lIB2zjnnlMhJeWOK\nsiJgnFq1ahXXXXcd9913H8nJya7jBLw6derw7LPP2mcRmxJjRcA407VrV/Ly8hg6dCjx8fH28YvF\nMGTIED7++GNWrFjhOoopI2x2kHEmKSmJpKQkAPr27es2jBeysrLIysry68Vu5cqVC+qL60zgsSMB\nY7zUuXNn/vvf/9KxY0fXUYzxmhUBc1o7duzg559/dh0jIOXl5ZGSksIVV1zh920fOHCAyZMn8+uv\nv/p926ZsOWMREJHpIpItIulF2mqJSLKIfCsiS0QkvMhzo0Rkk4hsEJEeRdoTRCRdRDaKyPMl/6WY\nkvbJJ58QGxtLkyZN2Lt3r+s4xqNevXrce++9PP300yxdutR1HBPkinMk8CrQ84S2kUCKqp4LpAKj\nAESkNXA90AroBUyW3+6oNQUYpKotgBYicuI6TYDZs2cP3bt3p3r16rbHGUDKly/PmDFjOO+881xH\nMWXAGYuAqn4CnLgb2A94zbP8GnBsvtpVwCxVzVfVLcAmIFFEooEaqvqFp9/rRV5jjDHGEW/PCdRT\n1WwAVc0C6nnaGwLbi/TL9LQ1BH4s0v6jp82YoLN7926GDBnCvn377NbRJuiV1BRRLaH1HFf0roxF\npxIaN77//nsaNgzNup2fn09GRgb16tXjl19+YcqUKXzyySe89dZbTqdrxsTEcM011zBs2DCee+45\nZzmMO2lpaaSlpfm0Dm+LQLaIRKlqtmeoZ6enPRNoVKRfjKftVO2nVJq35jVnZ+jQoXz99df07t3b\ndRQnJk6cyKOPPkq9evVo1aoVu3bt4v777+fPf/6z01wvvfQS3bp1Y8aMGU5zGHdO3EEeN27cWa+j\nuMNB4vl3zAJgoGd5APBekfYbRKSiiDQFmgErPUNGOSKS6DlRfHOR15gA98ADD/Duu+9ywQUXuI7i\nRG5uLtdddx25ubkcPXqUkSNHcsMNN7iORbly5ahYsaLrGCbInfFIQETeApKASBHZBjwEPAG8KyK3\nAlspnBGEqq4XkXeA9cARYIiqHhsquhOYAVQGPlDVxSX7pRhjjDlbZywCqtr/FE91O0X/8cD4k7Sv\nAtqeVTpjzGnVr1+fxYsXEx0dzcaNG6lZs6brSCbI2BXDJmBMnz7dPjDlLHXo0IGDBw8iIuTl5bmO\nY4KQFQHzOwcPHuTaa6+le/fubN261W/bffjhh6lduza33HILo0ePDsjPFsjLy+P77793HeMPwsLC\nbKqq8ZoVAfM7e/bsYdmyZeTl5bFmzRq/bTchIYEpU6bw4IMP8sILL7Br1y6/bbs4jt2qoVevXnTq\n1Ml1HGNKjN1K2vxBpUqViIyM9Pt2K1asyKBBgxg9erTft30mYWFhjBkzxnUMY0qcHQkYY0wIsyJg\njDEhzIqAMWXE4sWLadeuHbfddpvrKCaIWBEwpoz4/PPPiY2NZdGiRa6jmCBiRcCcVKVKlVi0aBGV\nKlVyHcWchTp16riOYIKMzQ4yxx09epQffvgBgMmTJ7N27VratnVzkXdeXh4LFy4kNjaWNm3aOMlg\nTCiwImCOe+GFF3jkkUfo1asXUVFRREVFOcvy/PPPs3DhQkTErxetGRNqbDjIHLdv3z5GjBjBzJkz\nXUfhyJEj9OrViyNHjriOYkyZZkXAmNPYuHEjmzdvdh3jjCpXrszChQupXLmy6ygmyNhwkDGn0b17\nd5o0acLf//5311FOKyUlhW3btnHOOefwzjvvuI5jgogVAWNOIT09nQMHDjBr1izq16/vOs5pxcbG\nEhsby44dO1xHMUHGhoOMOYnU1FQ6d+5M+/btiYiIcB3HmFJjRwLGnMSBAwe49NJLWbhwoesoxpQq\nOxIwxpgQZkXABKRvv/3WdQRjQoIVARNw/vWvf9G0aVP69evHgQMHmDZtGocPH3Ydy5gyyc4JmIBz\n1113AZCfn88//vEPHnvsMRo0aEDfvn39sv0vv/ySuXPn+mVbxrhmRcAErLCwMMaNG8eqVav8ut2x\nY8dSsWJF7r77br9u1xgXbDjIGI8DBw6QmJhIamoqd999N127dnUdySu//PILr7/+ut1ywxSLT0VA\nRO4SkbUiki4iM0WkoojUEpFkEflWRJaISHiR/qNEZJOIbBCRHr7HN6bk5OTksGXLFjIyMoL2w+Tr\n1KnD4MGDeeCBB1iyZInrOCYIeF0ERKQBMAxIUNXzKRxa+h9gJJCiqucCqcAoT//WwPVAK6AXMFlE\nxLf4pqRcfvnlPP744zRo0MB1FKfCwsJo3Lix6xheq1ChAk8++STx8fGoqus4Jgj4ek6gPFBNRAqA\nKkAmhW/6XTzPvwakUVgYrgJmqWo+sEVENgGJwAofM5gSsH79ejZs2EBcXJzrKMYYP/L6SEBVfwKe\nBbZR+Oafo6opQJSqZnv6ZAH1PC9pCGwvsopMT5sJENWqVcMOzowJLV4fCYhIBNAPOAfIAd4VkRuB\nE49BvTomHTt27PHlpKQkkpKSvMppjDFlVVpaGmlpaT6tw5fhoG7A96q6B0BE/g+4GMgWkShVzRaR\naGCnp38m0KjI62M8bSdVtAgYY4z5oxN3kMeNG3fW6/BldtA2oKOIVPac4L0cWA8sAAZ6+gwA3vMs\nLwBu8Mwgago0A1b6sH0TIsLCwnjkkUd48803+eyzz5gxYwZ79uzxaZ35+fns37+f/Px8ABYuXEin\nTp0IDw8/wyuDx1tvvUXv3r0ZNGgQhw4dch3HBChfzgmsBOYAq4FvAAGmAU8C3UXkWwoLwxOe/uuB\ndygsFB8AQ9SmL5himDRpEt26dWP27Nn85S9/4fHHH+ell17yaZ133HEHkZGRDBkyBIBVq1bRp08f\nli9fXhKRnRs0aBAVKlSgffv2zJs3j59//tl1JBOgfJodpKrjgBOPP/ZQOFR0sv7jgfG+bNOEnkaN\nGtGxY0fS09MpKCggKSmJgoICn9a5Y8cObr31Vn766Sfef/99li9fziWXXELt2rVLKLVbV199NVdf\nfTUA06dPd5zGBDK7YjjEHT16lDlz5nDw4EHXUZz529/+RvPmzbnxxhtdRzHG7+zeQSFu8eLFDB8+\nnJtuuonIyEjXcU7r8OHDPh8BnMqYMWMC/iMkjSkNVgRC3NGjR+nQoQMvvvii6yin1aRJE77++mvC\nw8OpUaOG6zjGlBk2HGSCQtu2bcnOzmbjxo1Ur17ddRxjygwrAibo1KxZkxdffJEBAwb4vK4NGzaQ\nm5tbAqmMCU42HGSCzogRI2jfvj39+/f3aT2XXXYZcXFxhIeHEx0dXULpjAkuVgRM0AkLC6NFixY+\nr6datWr861//KoFExgQvGw4ypoyLiYmhadOmvPnmm66jmABkRwImZOzZs4eFCxdy7rnnuo7iV59+\n+iljx47lu+++cx3FBCArAiZkvPjii7zxxhscPHiQ+Ph413H8JiwsjLAw+1M3J2fDQSbopaenU7du\nXZo1a3bamT4FBQV07dq11C44MyYYWREwQW/r1q20b9+e3Nxc8vLyXMcJWGvXrmXVqlWuY5gAY0Ug\nhGVkZLB06VLXMXwyY8YMxo0bR9WqVSlXzn6dT6Vv374cOXKEfv36uY5iAoz91YSwe++9l4yMDP72\nt7+5juK1OXPm0KdPn5Pe9iIrK4ulS5eSlZXlIFlg6dChA1OnTuXo0aOuo5gAY2eLQpiqMmLECPr2\n7es6ik8uuOACGjRoAEBubi4rVqwgLi6Ohx56iBUrVtCxY0fmzp0LQLly5di7dy+rV6+2IwdjsCJg\ngtj+/fvJyMj4XdvEiRNZtGgRqsr5559P3759WblyJTfddBMrV67kr3/9K3PnziUnJ4cePXo4Su7W\nvn372Lx5M3FxcURERLiOYxyzImCCUv369XnmmWfIycmha9eux9vz8/Pp2bMnixYtAiAxMZHDhw9z\n6NAhHn74YXr37k3NmjVdxQ4IgwcPJiUlhW7dujFr1izXcYxjVgRMUBIRBg8efMZ+UVFRvPrqq35I\nFBzy8vJIT0/nmmuuYceOHa7jmABgg6LGhIh69eoxatQorrzySjp37uw6jgkQdiRgTIgoX748999/\nPwALFy50nMYECjsSMMaYEGZFIAQdPHiQvn378sknn1ClShXXcUpUTk6O6wjGBBUrAiFoz549fP75\n56SkpHDZZZe5jlNikpKS+PLLL0lISHAdxZig4VMREJFwEXlXRDaIyDoRuVBEaolIsoh8KyJLRCS8\nSP9RIrLJ0z80J2kHiEqVKtGhQwdExHWUEvP222/z3Xffcc0117B3715WrlxJhQoVXMcyJqD5emJ4\nIvCBql4nImFANeB+IEVVnxKR+4BRwEgRaQ1cD7QCYoAUEWmuqupjBmN+p169eqSmppKbm0u3bt1c\nxzEmoHldBESkJtBZVQcCqGo+kCMi/YAunm6vAWnASOAqYJan3xYR2QQkAiu8Tm/MKVx00UWuIxgT\nFHwZDmoK7BaRV0XkKxGZJiJVgShVzQZQ1Sygnqd/Q2B7kddnetqMMcY44stwUBiQANypql+KyAQK\n9/hPHN7xarhn7Nixx5eTkpJISkryLqUx5g9q1arFhx9+SPv27Vm2bBnVq1d3Hcl4IS0tjbS0NJ/W\nId4OyYtIFPCZqsZ6HneisAjEAUmqmi0i0cBHqtpKREYCqqpPevovBh5S1T8MB4mInSooRZmZmSQm\nJpKZmek6inEoMzOTdu3a8dxzz/HnP/+ZatWquY5kfCQiqOpZzfbwejjIM+SzXURaeJouB9YBC4CB\nnrYBwHue5QXADSJSUUSaAs2Ald5u3xjjm4YNG3LbbbcxduxYnn/+eddxjCO+zg4aDswUkQrA98At\nQHngHRG5FdhK4YwgVHW9iLwDrAeOAENsd98Ytx577DEqVKjAkSNHXEcxjvhUBFT1G+CCkzx10nl5\nqjoeGO/LNo1vDh48yPLly13HMMYECLuBXIiZMGEC06ZNo3///q6jGGMCgN02IsQcOnSIW265haef\nftp1FGNMALAiYIwxIcyKgDHGhDArAsYYE8KsCBhjTAizImCMMSHMioAxxoQwKwLGGHbu3MmuXbtc\nxzAOWBEwJsRdcMEFpKSk0KlTJ9dRjANWBELAunXraNKkCfHx8Rw4cMB1HBNg+vTpw7Jly9i/f7/r\nKMYBKwIhYPPmzTRv3pyCggJeeOEFoqOjXUcyxgQIu3dQiKhSpQpffPEFBw8eJDw83HUcY0yAsCIQ\nQipVqkSlSpVcxzDGBBAbDjLGmBBmRcAYY0KYFQFjDABHjhxhzZo1rmMYP7MiYIwhIiKChIQEOnXq\nRGpqqus4xo+sCBhjqFKlCsnJyVx66aV2LUmIsSJgjDEhzIpAGbd37142btzoOoYJEuHh4QwePJjx\n48e7jmL8xK4TKONuv/12Vq9ezfDhw11HMUHg5Zdf5o033mDRokWuoxg/sSOBMu6XX35hwoQJDBs2\nzHUUEwSqVatGu3btWLp0Kc2bNyc3N9d1JFPKfC4CIlJORL4SkQWex7VEJFlEvhWRJSISXqTvKBHZ\nJCIbRKSHr9s2xpS8Cy+8kJ07d5Kbm2tFIASUxJHACGB9kccjgRRVPRdIBUYBiEhr4HqgFdALmCwi\nUgLbN8aUsPDwcMqVs4GCUODTT1lEYoDewP8Wae4HvOZZfg242rN8FTBLVfNVdQuwCUj0ZfvGGGN8\n42upnwDcC2iRtihVzQZQ1Sygnqe9IbC9SL9MT5sxxhhHvJ4dJCJ9gGxV/VpEkk7TVU/z3CmNHTv2\n+HJSUhJJSafbhDlRQUEBEydOJCMjw3UUY0wpSUtLIy0tzad1iKpX79GIyOPATUA+UAWoAfwf0AFI\nUtVsEYkGPlLVViIyElBVfdLz+sXAQ6q64iTrVm9zmUKZmZm0bNmS0aNHM3ToUKpVq+Y6kgkyDRo0\nIDk5mTZt2mCn74KDiKCqZ/XD8no4SFXvV9XGqhoL3ACkqupfgYXAQE+3AcB7nuUFwA0iUlFEmgLN\ngJXebt+cWc2aNbnvvvusABivtGvXjsTERKZMmeI6iilFpXGx2BPAOyJyK7CVwhlBqOp6EXmHwplE\nR4AhtrtvTOB6//33GTt2LDt37nQdxZSiEikCqroMWOZZ3gN0O0W/8YBdj26MMQHCJgIbY0wIsyJg\njDEhzIqAMeaMtmzZwoYNG1zHMKXAioAx5pQqVKjAtGnTaNq0KYmJiT7PSTeBx4pAGfTuu+8SHx9P\n3bp1XUcxQW7EiBFMmzaNdevWkZSUZDeUK4OsCJRBa9euZcCAASxfvtx1FBPkqlevTt++fWndurXr\nKKaUWBEoo2rUqGEXiRljzsiKgDHGhDArAsYYE8KsCBhjTAizImCMMSHMikAZs3nzZn744QfXMUwZ\ntXfvXg4fPuw6hilBVgTKmMsuu4zvvvvOPoTHlLj4+HiGDx/O7bff7jqKKUFef6hMabIPlfFedHQ0\nX3/9NdHR0a6jmDJo0aJFTJ06lUWLFrmOYk7Crx8qY4wxJvhZETDGFFtkZCQpKSmcd9555OXluY5j\nSoAVAWNMsV100UVs376dXbt2sX//ftdxTAmwImCMOSt169YlNjaWhg0b8tprr7mOY3xkRaAM2bt3\nLwUFBa5jmBDw6aef8uCDD9p05DLAikAZ8dlnnxEdHU1kZCTVq1d3HceUcSKCyFlNQjEByopAGbFr\n1y569uzJhg0brAgYv/n444+ZP3++6xjGB1YEjDFeufnmm2nbti133HGH6yjGB1YEjDFeiY2NZeTI\nka5jGB95XQREJEZEUkVknYisEZHhnvZaIpIsIt+KyBIRCS/ymlEisklENohIj5L4AowxxnjPlyOB\nfOBuVW0DXATcKSItgZFAiqqeC6QCowBEpDVwPdAK6AVMFjuzZEzQO3jwIPPmzePo0aOuoxgveF0E\nVDVLVb/2LOcBG4AYoB9wbPLwa8DVnuWrgFmqmq+qW4BNQKK32ze/eeyxx7j99tuJiIhwHcWEmMjI\nSG688UaGDh3KkiVLXMcxXiiRcwIi0gSIBz4HolQ1GwoLBVDP060hsL3IyzI9bcZHn332GePGjWPq\n1Kmuo5gQU7FiRV566SWuuOIK+vfvz+jRo11HMmcpzNcViEh1YA4wQlXzROTE2396dTvQsWPHHl9O\nSkqyWyOfQmZmJnl5eTRo0ICqVau6jmNC1NSpU+natSuzZ892HSWkpKWlkZaW5tM6fCoCIhJGYQF4\nQ1Xf8zRni0iUqmaLSDSw09OeCTQq8vIYT9tJFS0C5uR2795NixYtaNq0KS1atHAdx4SwihUrUqtW\nLXJycti8eTNxcXGuI4WEE3eQx40bd9br8HU46BVgvapOLNK2ABjoWR4AvFek/QYRqSgiTYFmwEof\ntx+yMjIyGDVqFBEREaxdu9aKgHGuVatWHD58mLZt27J7927XcUwxeX0kICKXADcCa0RkNYXDPvcD\nTwLviMitwFYKZwShqutF5B1gPXAEGGKfHOO9t99+m23btvHKK6+4jmIMAHFxcaxYsYIGDRrYR1AG\nEa+LgKouB8qf4ulup3jNeGC8t9s0v3fxxRfTs2dP1zGMMUHMrhg2xpgQZkXAGGNCmBWBINSzZ0+e\nfPJJ+zB5E/CysrLYunWr6xjmNCQQz82KiJ0zPo3o6GjS0tI499xz7Z7uJuC0aNGCsLAw8vPz2bRp\nE1WrVmXx4sV07tzZdbQyT0RQ1bN6U/D5YjHjP6rKl19+yeHDh4mIiLACYALSxx9/THZ2NgBRUVHc\ndttt7Nu3z3EqcypWBIJISkoK119/PZ07d7b7BJmAFR0dbUOVQcSKQJDYvn07aWlpdO7cmQULFriO\nY4wpI6yAqMyqAAAJ9UlEQVQIBIlhw4axY8cOhg0b5jqKMWft0KFDFBQUUK6czUUJNPYTCRL5+fmM\nGTOGm266yXUUY85KixYt6N+/v30MZYCyImCMKVXPPvss8+bNIzMzk+TkZCZMmMBPP/3kOpbxsOEg\nY4zfDBw4kMjISPbt2+fVHS9NybMjgSCwceNGm2Jngt6GDRvYv38/F198sesopgg7Eghwu3btIj4+\nnvj4eFq3bu06jjFeueyyy7jjjjsIDw8nMzMTuxg0cNgVwwEsJSWF++67j5ycHL777jvXcYwpEU89\n9RTPPfccO3fupE6dOnz44Ye0bdvWdawywZsrhm04KID997//JT4+ntTUVNdRjCkxI0aMYN68eaxf\nv5527dqxZcsW15FCmhWBANe4cWMaN27sOoYxJaZSpUpcfPHFtGzZkooVK5KRkcGePXtcxwpZVgSM\nMc707NmT6dOnc+211/LVV1+5jhOSrAgEqJ9//pmff/7ZdQxjStXQoUNJTk6mcuXKXHLJJezcudN1\npJBjRSAArVmzhsTERD744AMSEhJcxzGmVDVu3Jj//Oc/1KpVi6eeeoqYmBgGDBjgOlbIsCIQIHJz\nc5kwYQKNGjXi/PPPJy4ujnXr1nHllVe6jmaMX/Tt25f58+czYMAAkpOTXccJGTZF1LH09HSeeOIJ\n3n77bQCeeeYZ+vfvT/369R0nM8aNHTt2kJCQwI4dO1xHCTr2oTJBaO7cueTm5rJs2TIuvfRS13GM\nca5q1ar8+uuviAgxMTHExsbStGlT7rrrLtq0aUNYmL1tlSS/fzdF5ArgeQqHoqar6pP+zhAoZs6c\nyYcffkj37t2tABjjER4eTlZWFvv27WPLli1s2rSJZ599lk6dOjFw4EC2b99OQUEBTz/9NM2aNaN8\n+fKuIwc1v54TEJFywItAT6AN8D8i0tKfGbyVlpbm9WsLCgrIzs4mOzubkSNH0qhRIx5++GHuuece\nLr30Um699Va/ZyotgZgJAjOXZTq1SpUqERUVxYUXXkhMTAyrV69m3rx5LFu2jH379pGTk0PLli1p\n3rw5U6ZMYcqUKX6dWRQo36eS4O8Tw4nAJlXdqqpHgFlAPz9n8IovP/TJkycTGxvL+eefz5tvvsm9\n997LY489xq+//sp9991Ho0aN/J6ptARiJgjMXJapeI5l6t69O+np6aSlpbFs2TJ2795Nnz59SE9P\n56WXXiIqKoorrriCyZMnEx8fzzXXXMOMGTNK5T5Fgfh98pa/h4MaAtuLPP6RwsIQ0NavX/+Hu3gW\nFBRw+PBh3n//fcqVK0dcXBytW7fmxx9/5IcffmDr1q1s3bqVjRs38tVXXzFq1ChGjx59/PUDBgwg\nLCyMatWq+fvLMaZMiIyMZNKkSQCoKhkZGVxwwQX88MMPXHvttRw6dIg777yTQ4cO0bZtW5YtW8bh\nw4fp3r07LVu2pEaNGgCEhYUhIuzfv58tW7bQpEkTDh48SHZ2Nq1atTp2spUKFSoAcPToUY4ePYqq\nInJW52ADUsidYRkzZgxTp04lLi6OunXrnrH/xo0b2bhxIwAzZsygWbNm/Pjjj+zatYuCggKio6Op\nXbs2GRkZFBQUABAXF0dsbCwxMTE0a9aMPn360KdPn9+tNzw8vOS/OGNClIhw7rnnMn78ePbt28c/\n/vEPatSoQZ06dZg9ezaPPvootWvXJioqikmTJvHzzz8TFhZGfn4+VatWpW7dumzduvUP661atSq/\n/PILAO3bt6d+/fosWrQIgOeee466desSERFBVlYW1atXp127dsycOZOKFSv69ev3hV+niIpIR2Cs\nql7heTwS0BNPDotIaMwPNcaYEna2U0T9XQTKA98ClwM7gJXA/6jqBr+FMMYYc5xfh4NU9aiIDAWS\n+W2KqBUAY4xxJCCvGDbGGOMfAXnvIBH5k4h8JiKrRWSliHRwnekYERkmIhtEZI2IPOE6zzEico+I\nFIhI7QDI8pTne/S1iMwVkZoOs1whIhkislFE7nOVo0ieGBFJFZF1nt+h4a4zHSMi5UTkKxFZ4DrL\nMSISLiLven6f1onIhQGQ6S4RWSsi6SIyU0T8fhZYRKaLSLaIpBdpqyUiySLyrYgsEZFizT4JyCIA\nPAU8pKrtgIeApx3nAUBEkoArgbaq2hZ4xm2iQiISA3QH/ji9wY1koI2qxgObgFEuQgToxYn5wN2q\n2ga4CLgzADIdMwJY7zrECSYCH6hqK+BPgNPhYxFpAAwDElT1fAqH1G9wEOVVCn+vixoJpKjquUAq\nxfy7C9QiUAAcq2IRQKbDLEXdATyhqvkAqrrbcZ5jJgD3ug5xjKqmqGqB5+HnQIyjKAF3caKqZqnq\n157lPArf1Bq6zATHdyR6A//rOssxniPIzqr6KoCq5qvqfsexAMoD1UQkDKgK/OTvAKr6CbD3hOZ+\nwGue5deAq4uzrkAtAncBz4jINgqPCpzsSZ5EC+BSEflcRD4KhGEqEbkK2K6qa1xnOYVbgf842vbJ\nLk50/oZ7jIg0AeKBFW6TAL/tSATSScKmwG4RedUzTDVNRKq4DKSqPwHPAtso3Dndp6opLjMVUU9V\ns6FwZwOoV5wXObtYTESWAlFFmyj8BXwA6AaMUNX5IvIX4BUKhztc5hpN4ferlqp2FJELgHeAWMeZ\n7uf33xu/XMJ4up+fqi709HkAOKKqb/kjUzARkerAHAp/z/McZ+kDZKvq154hz0C5DDYMSADuVNUv\nReR5Coc8HnIVSEQiKNzjPgfIAeaISP8A/R0vVkF3VgRU9ZRv6iLyhqqO8PSbIyLTAyTXYGCep98X\nnhOxkapaqp8DeapMInIe0AT4RgqvX48BVolIoqqW6t20Tvd98mQbSOHwwmWlmeMMMoHGRR7HEABD\ni55hhDnAG6r6nus8wCXAVSLSG6gC1BCR11X1Zse5fqTwKPdLz+M5gOuT+92A71V1D4CIzAMuBgKh\nCGSLSJSqZotINFCs94BAHQ7KFJEuACJyObDRcZ5j5uN5UxORFkCF0i4Ap6Oqa1U1WlVjVbUphX80\n7Uq7AJyJ53bh9wJXqeohh1G+AJqJyDmeGRw3AIEw8+UVYL2qTnQdBEBV71fVxqoaS+H3KDUACgCe\noY3tnr81KLzI1PWJ621ARxGp7Nnxuhx3J6uF3x+1LQAGepYHAMXawQjUewfdBrzgucL4V+DvjvMc\n8yrwioisAQ4Bzv9QTqAExqH8JKAisNRzg63PVXWIv0ME4sWJInIJcCOwRkRWU/gzu19VF7vMFcCG\nAzNFpALwPXCLyzCqulJE5gCrgSOe/6f5O4eIvAUkAZGec6cPAU8A74rIrRTOFLy+WOuyi8WMMSZ0\nBepwkDHGGD+wImCMMSHMioAxxoQwKwLGGBPCrAgYY0wIsyJgjDEhzIqAMcaEMCsCxhgTwv4/CXMA\nzWf+KoQAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10f504d30>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "N = 100000\n",
    "N1 = np.random.binomial(N, 0.4)\n",
    "X1 = np.random.normal(-2, 1, N1)\n",
    "X2 = np.random.normal(2, 1.5, N - N1)\n",
    "X = np.append(X1, X2)\n",
    "\n",
    "plt.hist(X, histtype='step', color='k', bins=200);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define the likelihood model.\n",
    "\n",
    "Note that I mapped parameters bounded below by 0 to ones on (-inf, inf) and likewise for ones on (0, 1) to (-inf, +inf).\n",
    "\n",
    "Instead of naively calculating $\\log(f * \\exp(p_1) + (1 - f) * \\exp(p_2))$, I use the log-sum-exp trick.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def logpdf(params, data):\n",
    "    f = params['f']\n",
    "    mu1 = params['mu1']\n",
    "    mu2 = params['mu2']\n",
    "    sigma1 = params['sigma1']\n",
    "    sigma2 = params['sigma2']\n",
    "\n",
    "    p1 = np.log(expit(f)) + norm.logpdf(data, mu1, np.exp(sigma1))\n",
    "    p2 = np.log(1 - expit(f)) + norm.logpdf(data, mu2, np.exp(sigma2))\n",
    "\n",
    "    return np.logaddexp(p1, p2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From this, we can construct the log-likelihood for the model with some helpful logging output"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def ll(params, data):\n",
    "    out = np.sum(logpdf(params, data))\n",
    "    if isinstance(out, np.float64):\n",
    "        print(out)\n",
    "    return out"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The objective function (that needs to be minimized) is given by"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "data = X\n",
    "\n",
    "def objective(x):\n",
    "    params = dict()\n",
    "    params['f'], params['mu1'], params['mu2'], params['sigma1'], params['sigma2'] = x \n",
    "\n",
    "    out = -ll(params, data)\n",
    "    if np.isnan(out):\n",
    "        raise ValueError(\"Illegal probability encountered at {}\".format(params))\n",
    "    return out"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, we can use the `autograd` library to get the gradient and hessian of the objective function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "gradient = grad(objective)\n",
    "hess = hessian(objective)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For the optimization, we use `scipy.optimize` with some terribly wrong initial values:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-442770.185088\n",
      "-325866.959365\n",
      "-293026.634957\n",
      "-235160.946548\n",
      "-240483.273636\n",
      "-225903.288555\n",
      "-222531.149637\n",
      "-221215.701022\n",
      "-220459.306516\n",
      "-220348.844488\n",
      "-220213.02008\n",
      "-219992.734087\n",
      "-219757.599445\n",
      "-219622.435724\n",
      "-219534.412677\n",
      "-219525.28318\n",
      "-219525.039083\n",
      "-219525.010259\n",
      "-219525.008819\n",
      "-219525.008748\n"
     ]
    }
   ],
   "source": [
    "out = minimize(objective, jac=gradient, x0=[-3, -3, 3, 1, -1], method='L-BFGS-B')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The minimization has finished, and we can overlay the model on top of our data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x1033e9ef0>]"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAEACAYAAAC9Gb03AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xt8z3X/x/HHa+eZGXMYNmbMWTmPCCuSqFzqEtJBpRIq\nUldF/eK66rrq6qh0lSIdVAhFXSipkRZziBw2RjnPhm3M2PH7/v2x5VqL7Tu+332+3+9e99utm30/\n3/fn83lO89r7+/58Pu+3GGNQSinlubysDqCUUsq5tNArpZSH00KvlFIeTgu9Ukp5OC30Sinl4bTQ\nK6WUh7Or0IvIABFJEpHdIvL4ed5vKSLxIpIjIo+Uem+iiGwXkV9E5GMR8XNUeKWUUuUrt9CLiBcw\nA7gWaAuMEJFWpZqdAB4EXiy1b8Pi7Z2MMZcDPsBwB+RWSillJ3t69DFAsjFmvzEmH5gHDC7ZwBhz\n3BizCSg4z/7eQJCI+ADVgCOXmFkppVQF2FPow4GDJV4fKt5WLmPMEeBl4ABwGMg0xnxb0ZBKKaUu\nnlMvxopITYp6/5FAQ6C6iNzqzHMqpZT6Ix872hwGGpd4HVG8zR79gF+NMekAIrIY6AF8UrqhiOik\nO0opVUHGGCmvjT09+g1AtIhEFt8xMxxYWkb7kic9AHQXkQAREaAvkFhGYLf875lnnrE8g+a3Pofm\nd8//3Dm/vcrt0RtjCkVkPPANRb8YZhtjEkXk/qK3zTsiEgZsBIIBm4g8DLQxxiSIyELgZyC/+M93\n7E6nlFLqktkzdIMxZgXQstS2mSW+TgUaXWDfacC0S8iolFLqEuiTsQ4QGxtrdYRLovmtpfmt5e75\n7SEVGedxJhExrpJFKaXcgYhgHHQxVimllBvTQq+UUh5OC71SSnk4LfRKKeXhtNArVYl+/vlnBgwY\nwJgxY6yOoqoQLfRKVaLVq1fj7+/PzJkzy2+slINooVeqkjVp0sTqCKqK0UKvlFIeTgu9Ukp5OC30\nSinl4bTQK6WUh9NCr5RSHs6uaYqV5zIG1q+H336DBg2gd2/49dc9pKSk0K1bN/z8/KyOqJS6RFro\nq7Ddu+HeeyElBUQ289tv/vj51cGY+wkJSeSxxx5j4sSJVsesUrZt24aPjw+tW7e2OoryIHYN3YjI\nABFJEpHdIvL4ed5vKSLxIpIjIo+Uei9ERD4TkUQR2SEi3RwVXl28pCS46iq46SbYuRN8fG7npZdW\n0bPnQkSW0q3bC+Tk5Fgds0pZt24dPXr0oHPnzuzdu9fqOMqDlFvoRcQLmAFcC7QFRohIq1LNTgAP\nAi+e5xDTgWXGmNZAe8pYM1ZVjsxMuPZamDDhGPXrz2f//qKi0rdvX77+ehw//BDE118PITm5OpGR\nkQQHB7N+/XqLU3uepUuX0qdPH5544gkAsrOziYmJoWnTppw9e/ZP7VetWsVnn31Gfn5+ZUdVbs6e\noZsYINkYsx9AROYBg4Gk3xsYY44Dx0Xk+pI7ikgNoJcxZlRxuwLglGOiK3sdPnyYt956i6ioKLp1\n68aIEQW0bl2DtWsnsnXrVvbv30+1atWoUaMGAB07wrXXLmHRolvx83uF/v37s2/fPrp10w9jjrRq\n1Srq1KnDggULeP7558tse+zYMQYNGkS9evUwxnDLLbdUUkrlCewZugkHDpZ4fah4mz2iKPoFMEdE\nNovIOyISWNGQ6tIsWLCAVatWce+99zJ58hqSk+uwa9eNFBYWMmPGDDIyMkhJSaFRo/8t+9uq1XYa\nNPiVM2cexcdHL+U4S0REhF3tbDYbISEh9OzZk8LCQienUp7G2f+CfYBOwDhjzEYReQ14AnjmfI2n\nTp167uvY2NgqsZZjZYmJiWHdugR++ulG+vX7msTE/w0N1KxZ87z79Or1BbNnT+Tkyd8qK6bHeeqp\np5g3bx733HMPTz75ZIl3fEhLiyA7uwY5OUeZPXsJc+e+pr9UVZni4uKIi4ur8H72/FQdBhqXeB1R\nvM0eh4CDxpiNxa8XAn+6mPu7koVeOcMteHsXEB29ndWr08jLy0PkwstNVqt2mqCg19m+/RZAi/3F\niI+Pp2fPnqxbtw6AzMwgvv9+KPAMX36ZQZ06uaSnd2bMmFY0bNiaZ5+txgsvDLA2tHJZpTvA06ZN\ns2s/ewr9BiBaRCKBFGA4MKKM9ucqhzEmVUQOikgLY8xuoC+w065kyqGK1l2fQkDAv2jWrB0ffPAB\n2dnZ9O/fv8z9qlX7kBMnJvDbbxmVktMT1axZk4yMTN55B1566TZatfqRJ574hMzMHdx0003cf//9\ntGrVkVatpvDMMy1IT3+B7Gx9llE5TrmF3hhTKCLjgW8oGtOfbYxJFJH7i94274hIGLARCAZsIvIw\n0MYYcxp4CPhYRHyBX4G7nPXNqAs7fLgljRs35tFHL+e220ZecLimpDVr1iCSQ6tWS1i8+BqKbw5R\nFWSzCVu3jmf/fnjwwQWcPfsz//rXa0DRRdYTJ07w7bdfMm3aE/zzn9Cw4RlGjGjK6tVQ4rKJUhfN\nrgFBY8wKoGWpbTNLfJ0KnPdH0hizFeh6CRmVA2zZ0ospU2pw333j7Go/duxY6tSpQ5s2bZg161NW\nrryVX3+Fpk2dHNTDGOPF11/fwbFjOdx//3T8/Qs4WOLWhrp163Ls2DFsNhsBAQEANGw4jeuuu5LY\n2AZ8/z00bnyBgytlJ/18WAWkpwdz5EgUI0fav0+jRo2YNGkS1113HT4+ucTG/sqMGc7L6Ini439i\n27Z7CAxsxJNP/siLL07lwIEDf2rn5+d3rsj/Ljf3X9x001Guuw5O6Q3J6hJpoa8CNm1qTYsWWwgK\nuvhj9O+/hw8+gNOnHZfLk+3atYurrlpAXl4PFiyw8fTTj9k1XAbw3HPPsWfPHjIynqZPHxgxAmw2\nJwdWHk0LvYeLi1vNDz9E0br1hks6Tt26Z7jySpg/30HBPNzWrV4UFj7Npk2RtGwZVqF9Bw8ezJAh\nQxCB6dPh5EmYObOak5KqqkALvYcbPvxVRIRJk3pd0nEOHDjAsGGnmTPHQcE8TGJiImPHjmX69Oms\nWhXP6NFBNGjwPM2b/7Hdd999Z/cx09PTSU09xNy5MGNGEPn5bR2cWlUVWug9XHb2X7jvvkC6du1y\n0cfo1asXb731Fu+9N5TkZBsrVx7QCc9KWbx4MYmJiUyZMoWJE09Rq9Z+FiwY8oc2H3zwAQ888ACT\nJk0q93jt27dn9+7dtGrVivr1c3j66Syysl7BZtN/sqri9KfGg9lscPbsAAYOzL6k44wfP56PP/6Y\n7OxMAgMX0b//xzRt2pTXX38dmw4en3P55ZdTWNiBPXuuYMKEvVxxxRV/eL93796MHTv2D1NNXEhM\nTAzbtm3DGENhYSHDhuUgcprk5GudFV95MC30Hmz9evDyOkWzZgUOO2b16l8QFjaBcePG88QTT3Di\nxAmHHdtdJSUlkZycjDFCXt7LdOmymOBgx37iEYHg4MfYsWMIJ0/6O/TYyvNpofcgx44do0OHDrRp\n04Z9+/axeDEEBi53yLG9vb3ZvXs3hw8vx99f6N9/MkGXchuPB+nXrx/79u3D2/tOwJvmzeMdevy1\na9eSm5uLj8+vREWtYeHCyxx6fOX5dAYlD5KSksLZs2epVasWP/20jpkzq5GfPw+RIeXvXI6uXbvy\n3nvvISKsW+fPZ585ILCHyMvL48MPF9CjRx38/Poh4rinykaOHMno0aO59dZbAWjb9nO+/XYGO3dC\nmzYOO43ycNqj9zD+/v4EBgby3//+Sk5OIZ9++gSNHfBopYgwePBgbrzxRm65Rfjss9/nz1EAH30U\nyGWX2fD2TnDocd955x1uv/32cxe//fyyGTQokWfOO/+rUuenhd4DhYeH88knGbRosYe//GWww4/f\nvj34+kJBQQeHH9sdGRPIG29U46mn8rDZbJw8edKp5+vfP5m1a2HrVqeeRnkQLfQe6IMPPqB37+d5\n7rmeTjm+CAwdCnl5Nzrl+O7m7Nm7iYnJJybGly5durBlyxbatWvntPP5+xfy+OMwbVrRvfZPP/00\nr732mt4BpS5IC70Hys72ZtMmb/r1c97/3qFDITf3xio/fJOVBWfOjOXRR7Px9fVlzZo1JCcn0717\nd6eed9iwU6xZY+P9939kyZIlTJkyhePHjzv1nMp9aaH3QKtWQY8eXNLcNuVp3x6gkC1bqu71/AMH\nDnDllfMR+ZY2bZz7G+/w4f+t9bNr1y6aNg0jM/MFPv20Pm3bttU7oFSZtNB7oJUroZz1RC6ZCPj7\nL+Wrr6ruPd3ffRdPYmI/3n03ijp16jjtPLGxsaSlpTF8+HAAMjIyaN++PZddtoZt29qSkxPstHMr\nz6CF3gN9/z1cdZXzzxMZuZU33viVWbNmOf9kLig+vgl16x7g7rudO0zTv39/1qxZw/Tp0/+wvUWL\nGuTmfsLx40Oden7l/uwq9CIyQESSRGS3iPxpzVcRaSki8SKSIyKPnOd9LxHZLCJLHRFaXVh+fh1S\nU38fWnGu7dvnEBQUxebNGbz44otMmjSJlJQU55/YBRgDy5a1omXLZZZlmD9/Pjt2jGb37mswJqD8\nHVSVVW6hFxEvYAZwLdAWGCEirUo1OwE8CLx4gcM8jK4VWymys7vSuzd4ezv/XN7e0KLFHnbujOKp\np57ihx9+YNky6wpfZVqxAry9DWFhOyzN0aYNdO0KOTnaq1cXZk+PPgZINsbsN8bkA/OAP9ycbYw5\nbozZBPxpUhURiQAGAlXz830ly87uQolF4p2uZctkfvutLX5+flx2WdV5NP+VV2DgwERErE4CDz8M\nZ8/excGDh3j22WeZr4sGqFLsKfThQIlVLjlUvM1erwKPAVX8RrzKcfp010oZn/9ds2a/kpIShTFV\n566PrVshMRGuuOLPywJWhs2bN//hdd++YEwwL774PUuXLuXuu++2JJdyXU69N05EBgGpxpgtIhIL\nlNn/mTp16rmvY2Njia3MrqkHSE31obAwhMrsWAcE5FG//n6OHavE3y4We+MNGDsWfHwq/wGlcePG\n8fnnn9O7d+9z27y8IDDwI7Zvv5p+/Q6yY4e1w0nKeeLi4oiLi6vwfvYU+sNAyclSIoq32aMncKOI\nDAQCgWAR+dAYc8f5Gpcs9KriNmyoTlDQJry8+lXqeaOidpCSMgBYU6nntUJmJixaBElJRc8rVLYr\nr7ySK6+88k/bAwLmsWfPI/TqtaTyQ6lKU7oDPG3aNLv2s2foZgMQLSKRIuIHDAfKunvmXK/dGDPZ\nGNPYGNO0eL/vLlTk1aUrKvSOnVTLHk2abKew8FqMcYEBayf76CO49loIq9gysE7n5ZVGREQyu3d3\ntjqKckHlFnpjTCEwHvgG2AHMM8Ykisj9InIfgIiEichBYCIwRUQOiEh1ZwZXf5Sfn8+6dQFUr35p\ni4BfjJo1TwAnOX48stLPXZmMgbffhgcegNGjR3PvvfdSv359q2MBUL9+ffbvn0JSUu/yG6sqx64x\nemPMCqBlqW0zS3ydCpS5PpoxZjWw+iIyKjuMG/dvUlPvZeDA2pac38fnaw4evBzw3PlWvvwyg8xM\nL/71r+GsW/cT//3vf887jGKFhIQEsrPP0qlTTby8utKwYUMWLFjgMvmUtfTJWA+RlBRG587ZzJ5t\nzV2s3t7fcPiw82ZsdAUPP5xEQcGbRESEs3DhQvr06YOPj2vM9RMQEEDt2rUYPVoYPvxrYmNjSUpK\nsjqWchGu8VOqLllaWjM6dkwFoiw5v7f3T2RmNiAryzPnvklLg8OHL+c//znC6NGTrY5zQbfdBjEx\n/lx/vc5/o/5He/QeIi0tmmbN0iw7v0geDRrsYufOijxi4T4++gjCwzdSvXq+1VHKFBUFrVrh8Z+u\nVMVoofcAp05BVlZdGjdOt+T8q1cXXXoJD9/O9u0RlmRwJGMMe/bsITU1tfg1vP8+REXFWZrLXrfd\nBnv2dCc1NdXpq10p96CF3gMkJEDt2gcseYBn9OjR3HnnnXz66adERGxj+/YI3H2ho4ULF9KpUyea\nN29Obm4uP/8Mp09D3bruMeZ9yy1w9OjlvPXWJ1x//fVWx1EuQAu9B4iPh7p191hy7vr16zNx4kSu\nv/56goNPEBh4hoEDn2Lo0KEcO3bMkkyXKjMzk2HDhlFQUEB+fj5vvplNtWqf8cMPq13m4mtZatWC\nAQP8efjh1WzatIkmTZqwZ481Px/KNWih9wA//QT16u21OgYADRv+wrp1tUhKSmL79u1Wx7lkr776\nJu+9d4bU1Bd59913GTJkiNWR7HL77bB8eR0OHjxIREQE+/fvtzqSspAWejc1cuRI2rZtyxdfLGXd\nOtcp9JGRieTnX+3UFZcqU0JCHZo0Ocvx4wkMGjQI78qY/9kBBg4smnztzJna+Pt75p1Qyn5a6N3U\n119/TadOnVi+/DdCQyEw8JTVkRARkpNnk5PTnPx8z7i9b+fOGDp2/MXqGBXm7w833wzz5lmdRLkC\nLfRurHbt2uzb15DU1MV88MEHhIdbe2vjlClTePDBe7niinwyMrpYmsURjKnDoUPRtGvnHhdhSxs2\nDD77zOoUyhW4/pUlVaZDhxqRnz+fM2fOWP4RPSoqijFjxlBYCP/4R4ylWRyhsPAWmjXbTkBAntVR\nLkqfPrBvH0RHu8Z8PMo62qN3cykpTfD2Xm95kS9pwABIT49x+9ssCwpuQ+Qjq2NcNB8fuOkmSEvr\nY3UUZTEt9G4sJ6cap0/XxMvLtZbjbdYMfHyy2bvXfScwPXEinFq1mjNkSE1GjhxpdZyLdsstRYX+\n7NmzFBYWWh1HWUQLvRs7erQx9eodRMT1/gGHhiaQkGDNTJqOkJzcgzFjgnjuub8TGem+0y/37g0F\nBQ246abHGDdunNVxlEW00LuxtLRGhIVZs25peUJD15OQEGp1jItSUCDs3duNO++0Osml8/GBUaNq\nMHz4Io4cOWJ1HGURLfRuLC2tEfXqHSy/oQVq1tzK3r3VycyE/v374+/vz4wZM6yOZZedOxsRHHyM\nFi2sTuIYQ4fCjz82tDqGspBdhV5EBohIkojsFpHHz/N+SxGJF5EcEXmkxPYIEflORHaIyDYReciR\n4auigoICvvjiC3Jzc0lNdd1C7+2dR7t2J/n2W9i9ezejR49m717XeKirPPHxzWnRIt7qGA7Tuzek\npwdw+nQDq6Moi5Rb6EXEC5gBXAu0BUaISKtSzU4ADwIvltpeADxijGkLXAGMO8++qgIWLVrEgw8+\nyE03jaGgIAA/P9cs9AAxMSdYvrzoa3eYIwYgPR0SE8OJitpodRSH8faGbt2OcvToFezevZvly5eT\nmZlpdSxViezp0ccAycaY/caYfGAeMLhkA2PMcWPMJooKe8ntR40xW4q/Pg0kAp45YXklKSgooHfv\n3gwf/iLR0Vl8/vliOnXqZHWs8+raNZ0VK4qm+XUXb72VQePGO/H3P2t1FIfq3j2Fo0e7cf311zNm\nzBieffZZqyOpSmRPoQ8HSnYbD3ERxVpEmgAdgPUV3Vf92caNMHBgPY4fP35uPnhXExFxloAAyM93\nnw9x//znEWy2OfTv39/qKA512WXHOXkynGPHfBgyZAj5+a69gIpyrEr5PC0i1YGFwMPFPfvzmjp1\n6rmvY2NjiY2NdXo2d7VxI7j67d0icN118PHHsUCK1XHKtXMn5OXV4bPP7qN9e89aoemqq3oSHZ1M\ns2bPExGxm4MHXXfIT11YXFwccXFxFd7PnkJ/GGhc4nVE8Ta7iIgPRUX+I2PMkrLaliz0qmybNsGr\nr1qdomwHDx7kqqtymDUrFvjU6jjl+uADCAn5Cm/vblZHcbjQ0FCefTaU2bM74+s73eo46iKV7gBP\nmzbNrv3sGbrZAESLSKSI+AHDgaVltJdSr98Ddhpj9KfLQc6eDeHMmaL1QV3VlVdeyRNPPMHChePJ\ny2tHbq7rTNFwPoWFMHcu1KxZZl/ErQ0YULRIzdmzflZHUZWs3EJvjCkExgPfADuAecaYRBG5X0Tu\nAxCRMBE5CEwEpojIARGpLiI9gZHA1SLys4hsFpEBzvt2qob09Cg6dy4aGnFV//jHP3jrrbc4c+YY\n/v6bOXSopdWRLig/P59Fi07SsKGNgIBfrY7jNMHBRbda7tzpwj0E5RR2jdEbY1YALUttm1ni61Sg\n0Xl2/RFwj5Ua3Eh6elMGDbI6hf0CA1ezb19fmjd3zRWnxowZw/vvX0unThe8fOQxhgyBGTOa0ru3\n1UlUZdInY91QenoTOne2OoX9AgPj2L+/NZs3/8zSpWWN+lnj8OFsfHyup1q1L62O4nQ33ABJSZEU\nFLjHcw3KMfT/thuw2Wx8+umniAg2m4309KZ0caN1PXx89hIUFEB4+DUMHjwY42I31h850ovGjXfh\n65tldRSnq1cPGjY8zsGDrjuUphxPe/RuYP369UycOJGHH36Y1auTKSz0wZ0mVBSBfv3y6dhxstVR\nzuvAgWto02ad1TEqzeWX72Xv3susjqEqkRZ6N2CMITo6mqioKFJTIwgN/c2lL8SeT58+Z1mxwuoU\nf7ZtG+TmhtK48S6ro1Sayy7by2+/tUOnp686tNC7mbS0RoSG/mZ1jArr0SOHhAQA11o0fM4ciIhY\nhZeXaw0nOVPt2qcICjpJvOfM26bKoYXejXh7e7N9uz916+63OkqFBQUZevUCuM7qKOfk5RXdO9+o\n0Sq8vb355ZdfOHDgAN7enn+jWLNm2/j8c6tTqMqihd6NzJnzPn5+V/DCC7dYHeWiDBkCcJPVMc75\n6itDaGgqaWnx9OnTh1mzZrFgwQJatXKfuXkuVlTUNr780r0mnFMXT++6cSPBwS3x8YGYmECro1yU\nG2+E++67lpwcCAiwNktOTg7PP3+S/fuf4YEH+tOrVy9CQ91zRayLUbfuYXJyYNcuqAK/16o87dG7\nkU2bcPknYktq0KABy5YtIzMzk+DgYOrVA9jKt99anQwee+xVNm7057bbAnjllVeqVJEHWLNmNY0b\n/8KXnv/ogEILvVvZuBG3un8+JiaGM2fOcPz4cWrX/n2h8MUsXmxpLAA2b25L9+6HePfd16yOUulu\nv/12Ro0aRXz8E1roqwgt9G7k9x69O/H19S21utTnLF0KBQUX3MXpjIGdO7vRrdtO60JYKDQ0lAkT\nJgDfsXUrnDhhdSLlbFro3YQx7tejP7+DREXBDz9Yl+D32wqbNDlqXQiXkMtVV8F//2vj3Xff5ZVX\nXuHMmTNWh1JOoIXeTeTl1cMYiIiwOsmlu+kmzg3fvPDCC7Ru3Zq///3vlXb+t96Cdu3i3eZahzPd\ncAMsWHCGRx99lDfffJNvXeECinI4LfRuIju7FV26uM+F2LLk5y9g0SIbNhusXbuWLl26sGbNmko5\nd1oa/Pe/0Lp11ZnyoCzff/8ocXH+1KoVRtu2ba2Oo5xEC72bOH26lQcM2xTNUz9//jS8vE6cG76p\nrDte9u3bx1/+8hVNmmzi7NkjlXJOV7Z8+XKOHNlM9epHyMmJsTqOciK7Cr2IDBCRJBHZLSKPn+f9\nliISLyI5IvJIRfZV9snObul2F2JLExGeeuop+vbtS+fOu/jkk8o79+TJk+nXbwA//xzDgQOPk5CQ\nQIcOHSovgAsaMGAAl112GU2bbufMmb5Wx1FOVG6hFxEvYAZwLdAWGCEipR+xOAE8CLx4Efuqchjj\nOT3633XsuJtFi8Bmc+wze9OnTycsLIxRo0b9Yfsnn3zCDTfM5PLLa3PixLfs2rWL3rr6Br6+vuzf\n/yZnz/bTp2Q9mD09+hgg2Riz3xiTD8wDBpdsYIw5bozZBJS+aa7cfVX5UlP9ELHRsKHVSRzD29ub\nL754nRo1DpKW5tiPKb/88gs33ngjmzdvpqCggKysLAqLp2ncsCGGhx/2/HlsKmLKlCm8/fZY6tat\nz+nT51skTnkCewp9OHCwxOtDxdvscSn7qmJJSdUJCkryiAuxUDSMcs8995CV9SYHD/Zz+PGrVasG\nwJ133kloaCjjxo0jP78Zu3b5cfPNDj+dW6tVqxY33HA9Q4b4cvSojtN7Kr0Y6wZ27apO9epJVsdw\nmLp16zJw4ECqV19Geno7srOdM3VxSkoKY8eOJSUlhZMnx3D77afw93fKqdzejTdCaqoWek9lzwDp\nYaBxidcRxdvsUaF9p06deu7r2NhYYmNj7TyNZ0tKCiIoyHMKPYC/vz+pqXsQ+Zy9e3vg7e282x3P\nnKnD2bP9uPPOLKCW087jzmJjISurESdPut8U2FVJXFwccXFxFd7PnkK/AYgWkUggBRgOjCijfckB\nhgrtW7LQqyLGFPXoo6M9awWkiIgIEhISWL/em6lTI8nKuo0hQ4bw2GOP8X//9380bNiQOXPmOGRu\n+A0b+uDj8yEhIXp56EL8/aFOnV/YtCmM22+3Oo26kNId4GnTptm1X7lDN8aYQmA88A2wA5hnjEkU\nkftF5D4AEQkTkYPARGCKiBwQkeoX2rdC31kVt38/+PgY/PyOWx3F4dq1a8fdd7emTp1Ann56DZs3\nb2bBggX4+fkxf/58cnNzL/kcXbsOIi9vKM8+G0rjxo3L36EKq18/gYSEMKtjKCew6942Y8wKoGWp\nbTNLfJ0KnPeS/fn2VfZ78811GGPz2FWPRODBB4VFi9pRvXp1AMLDwx32/S5f3oI77wzksce0m1qe\nevU28eOP48jNRa9leBi9GOvivvrqKC1anOK9996zOorTjBgBGzZAbq5je9z5+UF89VVDHnvMoYf1\nWIGBWXh772Lq1O+tjqIcTAu9i8vIiKJXr2o0b97c6ihOExgIDzwAJ07c49DjHjgwgh49jtO0qUMP\n67FeeuklOnQ4xNy5mVZHUQ6mhd6FGQMZGc1o2jTD6ihO99BDkJXVj1OnajrkePn5YRw5MohRo35z\nyPGqgujoaCZMaMahQx2pVSuUI0d0PiBPoYXehe3bB97eedSqlWN1FKerXRtq1lzEpk3/m3MlMTGR\nxx9/nI8++qjCxzt69EEaNlxG3bqXfkG3KvnLX5rTpEkktWv34YSuSOIxtNC7qNzcXO66601stgR8\nfX2tjlPSGXTRAAAe/ElEQVQpateew+7dHcnMLLrz48MPP2Tt2rU8VsFB9sOHm5KdHUNk5FxnxPRo\nInDDDUJWVqzVUZQDaaF3URkZGaxbl8fQoc0YPLhq3P/t45NBly6rSEj467ltXS4wk1tWVhYvvvgi\nc+bMwZSYjauw0Ie4uKE0aPAiPj5nnZ7ZE91wA1roPYwWehdmTGdGjGjusbdWnk/79mvIzGxIYeGg\nMtt9+eWXvPfee4wdO5bjx4ueMVi4cCHffHMlISHHqVFjZWXE9Uh9+kBubhQnTjh2ZlFlHS30LsoY\nyM+/3O3noK8oH59CevV6n7y81zl7ttof3rPZbHz00UfMnTsXm81Gx44dCQ4OZuLEiaxdu5Zhw2Zy\n6tQNPPfcMY+ZAM4Kfn5Qvfo61qxxzhxEqvJpoXdR+/d7I5JNWBV8ULFBg2R8fBawcuVt2Gz/q9gJ\nCQlMmjSJiRMnsmnTJgA+++wzsrKyePfdOGy2j1m6tCZt2tS1KrrHCA7+ntWrtdB7Ci30LqiwsJDv\nvjuFr+9Wq6NYJijoWVJS0omP/+u5BTGMMTRr1oyoqKhz4/J9+vShefN+LFlyP35+/+Dqqy0M7UGC\ng38gPj6A4OC6iAiffPIJp0+ftjqWukha6F3Qyy+/zLRpX9Gs2Umro1hm/fq1LFhgyM3tSlbWS+Tn\nn7/dpk0wa9YdtG//A76+cyo3pAfz8TmJj08iMTGPM3HiRKZOncqbb75pdSx1kbTQu6D09HQiIv7C\nyy+XNUmo5wkNDWX27NmEhobSsmVLbryxN0uWZGKz1aZjR/j++1AKCqphjHDiRBibNt3OwIEwYMC3\ndOy4xur4HuWFF17giivSCA29i1deeYW//vWv51bqUu5HL6u7IGOEI0fCqtyF2GXLlpGenk6DBg3O\nbQsKMoSE3MW//nWUZ56pz9atX2Gz+ZKUlEVExFq2boVZsxI57nmTe1pq4MCBREcXzVNvs1mdRl0q\nLfQuKCMjlICAXOrWDbA6SqUKDg4mOPjPFwCLHuKBOnV2MmHCY9hshfTs2Y20tDTq17/BgqRVQ4sW\nEBICGzdanURdKh26cUFHj0YQHn7U6hguoVq1auTn5+Pv7098fDxeXoVERoYzffr0P8wvn56e/of9\nMjIyOHpU/w4v1eDBsGSJ1SnUpdIevQtKSYkgIiIFncYfatSoQVpaGiNGjODAgQMALFiwAJvtf3P0\nd+7cmblz59K/f38AmjVrxg033EBeXh4dOnSwLLsnGDwY7ruv6BOVcl929ehFZICIJInIbhF5/AJt\nXheRZBHZIiIdSmyfKCLbReQXEflYRPwcFd4TTZgwgW3bAoiMPGZ1FJfh7e2Nv78/S5YswdfXFy8v\nL3x8fJDip6Kuu+46kpKSWLx4MQBBQUH85z//YdasWYSGhloZ3e3FxEBaGmRk6Fq77kxKzhNy3gYi\nXsBuoC9whKJ1YIcbY5JKtLkOGG+MGSQi3YDpxpjuItIQWAu0Msbkich84L/GmA/Pcx5TXpaqICCg\nOiIZ/PprAQ0aBFodx2WkpqayadMm2rRpQ5MmTayOU6WMHg379v2Xq6/eyuTJk62Oo0oQEYwx5T4H\nbk+PPgZINsbsN8bkA/OA0rNsDQY+BDDGrAdCROT3Zzq9gSAR8QGqUfTLQl2AMa1o1Mhbi3wpYWFh\nDBw4UIu8BQYPhuTk1lbHUJfAnkIfDhws8fpQ8bay2hwGwo0xR4CXgQPF2zKNMd9efFzPZ7N1pnNn\nvZ9NuY5+/SA1tSFnzmjnw1059a4bEalJUW8/EmgIVBeRW515TndnTBct9MqlBAZCZORedu1qZnUU\ndZHsuevmMFBy1eaI4m2l2zQ6T5t+wK/GmHQAEVkM9AA+Od+Jpk6deu7r2NhYYmNj7YjnWYp69Hqt\nQrmW5s0TSUxsZ3WMKi8uLo64uLgK72fPxVhvYBdFF2NTgARghDEmsUSbgcC44oux3YHXii/GxgCz\nga5ALjAH2GCM+dOkGXoxFnJyIDDwDMePC7Vr68dk5TomTHiOt9+eRGZmAAFV6zk+l+awi7HGmEJg\nPPANsAOYZ4xJFJH7ReS+4jbLgN9EZA8wExhbvD0BWAj8DGwFBHjn4r4lz7d1K4jsplq18tsqVZmq\nVcumXr1Uvv76ArPLKZdWbo++smiPHmbMgAkTZpOVdSuBgdqjV67jgw8+YNSo7UAzpk/P45577iEo\nKMjqWFWeI2+vVJVkwwYQ2WR1DKX+5M4772TnzucJDh7JG2+8ydKlS62OpCpAp0BwEcYYfvwxj6Ln\n0ZRyPa1bexMeHkyjRkOtjqIqSHv0LmLOnEXs3VtAr16h+Pr6Wh1HqfMaPBgOH65i82d7AC30LmLH\nDn/q1TvKd999jY+PftBSrqmo0HexOoaqIC30LmLfvnrUq/eb1TGUKlO3bpCbW52UFF043J1ooXcR\n+/bVpW7dfVbHUKpMXl7QqNFGpk8/xJw5ukavu9BC7yL2769LnTr7rI6hVLmmTbuc/PwbmTlzptVR\nlJ200LuAY8fgzBk/QkJ0Dnrl+oYODSMnpy45OQ3Kb6xcghZ6F7BhA0RGHkOkaj8wptyDjw/07p1O\nRsZVVkdRdtJC7wISEqBJE+3NK/cRG5tOeroWenehhd4FaKFX7qZLl1Pk5DTm0CGrkyh7aKG3mM0G\n69ZBs2apVkdRym4+PoaaNddSvEyvcnFa6C22axfUqgU1apy1OopSFRIa+h0LF1qdQtlDH8G0SHJy\nMvHx8Rw82J8rrtC7F5T7CQlJ4Jdf4OhRqF/f6jSqLFroLTJp0iSSk5PJyqrHgAEpLFy4UBe+Vm7F\nyyufQYPg88/hgQesTqPKokM3FrHZbAwaNIiUlKZ8/PF4OnfuzKRJk6yOpVSF/PWv6PCNG7Cr0IvI\nABFJEpHdIvL4Bdq8LiLJIrJFRDqU2B4iIp+JSKKI7BCRbo4K7+46d+6Hv38zVqx4kX/961+0bNnS\n6khK2SUkJIStW7fyt7+1Z8MGG8f0pjGXVm6hFxEvYAZwLdAWGCEirUq1uQ5oZoxpDtwPvF3i7enA\nMmNMa6A9kIgCYNeuWnTv7kOfPj2tjqJUhbRu3Zp9+/YREuJLtWprmDp1q9WRVBns6dHHAMnGmP3G\nmHxgHjC4VJvBwIcAxpj1QIiIhIlIDaCXMWZO8XsFxphTjovv3pKSanHFFVanUOrihIWF8eGHH9Kp\n024WLtTLfa7MnkIfDhws8fpQ8bay2hwu3hYFHBeROSKyWUTeERFdDLVYUlItevSwOoVSF69NmzYM\nGxZCZmYkhw9bnUZdiLN/DfsAnYBxxpiNIvIa8ATwzPkaT5069dzXsbGxxMbGOjmedYzxYvfumnTv\nbnUSpS6Nn5+N8PCNPPpoNSIiPuPuu++mdevWVsfySHFxccTFxVV4P3sK/WGgcYnXEcXbSrdpdIE2\nB40xG4u/Xgic92Iu/LHQe7qsrEhCQ3OpXVuXDVTur7DwE+bPv48ePX7C19eXf/7zn1ZH8kilO8DT\npk2zaz97hm42ANEiEikifsBwoPQS8EuBOwBEpDuQaYxJNcakAgdFpEVxu77ATruSebiMjFa0apVh\ndQylLtnAgQO5554mBAW1JSbmVqvjqPMot0dvjCkUkfHANxT9YphtjEkUkfuL3jbvGGOWichAEdkD\nZAN3lTjEQ8DHIuIL/FrqvSorPb0tsbHp/PGDkFLuJyQkhP/7v8mkpsK2bW3o2lVnOnM1do3RG2NW\nAC1LbZtZ6vX4C+y7Feh6sQE9kTGQnt6Oyy7bZnUUpRxmxAi4+ea2dO36jdVRVCn6ZKwFfvsNjBEa\nNMi2OopSDtOjB+Tl+ZGWphPfuBot9BZYvRpq196OiNVJlHIcLy9o334727Z1tDqKKkULfSVLT0/n\noYcWkZr6GXXr1rU6jlIO1bHjNnbsaE9BgdVJVEla6CtZeno6OTkxrF37T7p102l/lGepWzedkJBM\nVq60OokqSQt9JTp9+jQrVuzEZqtGt241rI6jlFO0a/czH3xgdQpVkk5QUYn+/ve/8/77+URHRyFS\n2+o4SjmciHDw4IusXn0tmZkB1KxpdSIF2qOvVGfPniU6ejTjx19mdRSlnGLMmDE8+OBICgpW8Nln\nVqdRv9NCX8n27GlInz5Wp1DKOWrVqsUdd9yBt/fHOnzjQrTQV6Ls7BpkZwfSrp3VSZRyLm/vb0hO\nht27rU6iQAt9pTp8OJqmTY/gpX/rysOJFHDHHTBrltVJFGihr1SHDrWkZcuD5TdUys0VFhbSseNG\nPvjAkJdndRqlhb6SGAMHDrSiVSst9Mqz+fv7M3jwYCZOHAQk0b//DB544AEyMnS2Vqtooa8ku3YV\n/Vmvnv6wK8/m4+PD/PnzWbhwIW3bxrNv3zXExcWxefNmq6NVWXoffSUwxvDUU2vw80vX+W1UldGr\nVy+WLetFRAQ0bNjJ6jhVmvboK4ExhkWLshg40Idbb9WFGVTVERAAt90GR45cZ3WUKk0LfSXIzQXo\nzWuv3UDt2vpErKpa7r0XUlIGUFCgH2etYlehF5EBIpIkIrtF5LxrvorI6yKSLCJbRKRDqfe8RGSz\niJRegrBK+OkngCS0xquqqG1bqFbtED/8oLO1WqXcQi8iXsAM4FqgLTBCRFqVanMd0MwY0xy4H3i7\n1GEepgqvFbtypQA6nZ+quiIiPueLL3TZTKvY06OPAZKNMfuNMfnAPGBwqTaDgQ8BjDHrgRARCQMQ\nkQhgIFAlH51YvHgxzz+/kfr1f7E6ilKWqVMnnmPH/Nm40eokVZM9hT4cKHnz96HibWW1OVyizavA\nY4C5yIxubfv2o/j5tWPfvnlWR1HKMl5eNgYPPsQbb1idpGpy6u2VIjIISDXGbBGRWKDMqzFTp049\n93VsbCyxsbHOjFcpdu1qQkTEXvz9dYIbVbUNGHCEe+5pTmoqhIVZncY9xcXFERcXV+H9xJiyO9oi\n0h2YaowZUPz6CcAYY14o0eZt4HtjzPzi10lAH4rG5m8DCoBAIBhYbIy54zznMeVlcUcdOuyiXr0t\nfPPNMKujKGWZvn37EhkZSXr683TuXI+nn7Y6kWcQEYwx5d7OZM/QzQYgWkQiRcQPGA6UvntmKXBH\n8Ym7A5nGmFRjzGRjTGNjTNPi/b47X5H3VNnZBezaFUnTpjusjqKUpSZPnkxWVhaZmf/grbcgJ8fq\nRFVLuYXeGFMIjAe+AXYA84wxiSJyv4jcV9xmGfCbiOwBZgJjnZjZbXTp8jdyczdy+eX1rY6ilKX6\n9u3L7bffTmHhFqKjT/PRR1YnqlrKHbqpLJ44dBMc/DEPPTSQ556rZXUUpSy3d+9ebrvtNjZvrk7D\nhsvZs8cHb2+rU7k3Rw7dqItgs8GZM9dwzTVnrI6ilEto1qwZP/30Ew0b7gXSeP/9k1ZHqjK00DtJ\nQgJ4e2cSFVVgdRSlXMqDD45H5N9MnnwKD/sQ77K00DvJF19AYOA3VsdQyuU88sgjPPhgFAUFPnz1\nVQ7p6el42rCtq9FC7yRLlkBg4NdWx1DKJXl5Qe3a73HjjT/RoEED5s3TBwqdSQu9EyQmwqlT4Oe3\nzeooSrmkIUOGMGhQJnXqdOT6618hNTXV6kgeTQu9E3zyCQwfDiL6cVSp82ncuDGvvvoi06fX5Kef\nBupYvZNpoXcwmw3mzjWEha0kOzvb6jhKubRhwyAvL4CdO5tYHcWjaaF3sPh4KCzM4rXXRjFs2DDq\n19eHpZS6EG9v6N59OcuX98BmszqN59JC72Bz50KfPgeJienKm2++ib+/v9WRlHJpzZr9gojhvvtW\n8Pzzz+snYSfQQu9AubmwcCH07n3I6ihKuQ0R6Nx5MXPmRDNr1oesWrXK6kgeRwu9g9hsNv72t++p\nUyeFWrWyrI6jlFs5eXIB/v4H8Pd/1OooHsmp89FXJdu2beOtt05TrdrLrF2rT8MqZa+RI0eSl5dH\nv37p/P3vN3Py5E9WR/I4Wugd5ORJbwoLY2nXbiYQbXUcpdxG165d6dq1KwCzZy9h7tzLCQpajJ+f\nH4MGDUKk3Dm7VDm00DvAokWLmDTpEDVrNsfH57TVcZRyWy1azOP77/uwYcMbFBZuJj4+nrZt21od\ny+3pGL0DfPXVMk6eHMHMmR2tjqKUW7v99utp3foTQkPn0ahRE2x6z6VDaKF3gJSUVlSv7sPNNzfA\n29ubZcuW4a0TbStVYSNHjmTz5rE0bBhGerouv+kodhV6ERkgIkkisltEHr9Am9dFJFlEtohIh+Jt\nESLynYjsEJFtIvKQI8O7AmMMO3f25aqrEhGBt99+m3/84x+8/PLLVkdTyi15ecG770Ja2hhSUnyt\njuMR7Fkc3AvYDfQFjlC0huxwY0xSiTbXAeONMYNEpBsw3RjTXUTqA/WNMVtEpDqwCRhcct8Sx3C7\nFaZycnKIiLiKEyeW8OWX27n++qutjqSUxwgLe52WLUexenUNRODYsWN8+eWXtGjRgiuvvNLqeC7B\nkStMxQDJxpj9xph8YB4wuFSbwcCHAMaY9UCIiIQZY44aY7YUbz8NJALhFfg+XFpeXh6nTt3FlCn1\ntMgr5WB16szh6FFfXn45ldTUVKZPn85LL73EzTffbHU0t2NPoQ8HDpZ4fYg/F+vSbQ6XbiMiTYAO\nwPqKhnRVx48L+flDGTfO6iRKeR4vrwKuumoOf/ubN82a9efs2bNcc801ukjJRaiU2yuLh20WAg8X\n9+zPa+rUqee+jo2NJTY21unZLsWMGX74+s6nQYPRVkdRyiMFBOyiW7dvSEiYTW7uJ1V+MfG4uDji\n4uIqvJ89hf4w0LjE64jibaXbNDpfGxHxoajIf2SMWVLWiUoWele2d+9efvwxkbffvho/v1cBLfRK\nOUubNt+xYUNN1q0bSM+eX1odx1KlO8DTpk2zaz97hm42ANEiEikifsBwYGmpNkuBOwBEpDuQaYz5\nfcmY94CdxpjpdiVyccYY+vTpw9/+lkr9+ut49dUJVkdSymNt3boVEfD3H0tSUhd++62N1ZHcUrmF\n3hhTCIwHvgF2APOMMYkicr+I3FfcZhnwm4jsAWYCDwCISE9gJHC1iPwsIptFZICTvhenW716Nd7e\n3uTkhJKffzerV1/Nvffea3UspTzSE088Qdu2bbnrrrsQOU6rVv/HypUjKShoYXU0t1Pu7ZWVxR1u\nr5w7dy4rVqwgKGguISHw739bnUipquGFF15g165d1KnzKC+/HEhaWhS1a1udynr23l6pc91UUHp6\nY1auLFoAXClVOR5/vOg5zbS0NGbM+IyhQ8exYgX4+VkczE3oFAgVYAxs2DCKZ5+F0FCr0yhVNQUF\n/YPgYLjjDkhPP8mJEyc4ceIEBQU6PfiFaKGvgNWrm2Gz+XD33VYnUapq8vf3Jycnm23b2rBjx1Hq\n1fuc5s1bEhkZyaRJk6yO57K00Nvp0CFYsKAj3bu/W+Xv5VXKKiEhIezZs4d27aK58sqXCQzsyO23\nH+fdd2eRmppa/gGqKC30drDZ4JZbThEd/TW1ah0sfwellNOEhYVRvXp1Tp06QosWj/DTTzBrVlds\nNl2g5EK00JfDGMNdd21n8+Y9NGz4PnfeeafVkZSq8nr16sW2bdsYMKA7q1ZBWloQ8fETycmxOplr\n0tsry3D48GGefXY1M2f2Y+jQl3j33aeoUaOG1bGUUqV89NECxo2rTt26lxMf35CwsKrRh3Xk7JVV\n1gMPvMDs2f0ZPz6B+fP/rUVeKRc1ePAAnnxyGwcPfkrz5id5//2dnDlzxupYLkN79Bdw/Dg0bXqI\nm28+xpw5ukSgUu5g/vz5/PvfyWzefD+hodPZv/9JqlcPsjqW02iP/hJkZMCAARAensB11yVbHUcp\nZadhw4axadNTbN9ek+zs4QQH/8icOausjmU5LfSlpKXBVVdBePgeqlX7p9VxlFIXoW1bX7Ky2tG5\nczYPPdSDt9+GwkKrU1lHCz1w9uxZPv30U+bMSaB7d+jX7zRff30ZnTt3olevXlbHU0pdBF9fuOGG\nXzh9OoaHHoqndu0DREc/wH/+8x+ro1W6Kj9Gn5iYyBtvvMGiRQGkpT1J3br/5Nix1+jQoQM///xz\npedRSjmOMYaMjAwmTXqU5OTW7NhxFz4+h5g/vwNXXQUiMHr0aPbt28cLL7xA586drY5cIfaO0VfZ\nQm+M4bvvvuOuu/6P3Nx/4uvbhddfT6VGjV+5/PLLCQkJwd/fv9LyKKWcb/Hirxg79keM+RtNmtTk\nkUeE4cP9+ctfBtKnTx8mTHCv9SV09soyFBYW8tVXqxgx4icKC5czaZI/Tz/tT2BgU6Cp1fGUUk5y\nzTV9GDfuZ156qRmdOi1g1KhARA6xb992mjY9gjFFvXxPY1ePvnixkNcoGtOfbYx54TxtXgeuA7KB\nUcaYLfbuW9zOqT36M2fOsG3bNqpVi2TUqAQ2b46hadNUVq5sT1Ot7UpVKdHR0fj6+nLFFVcwduzL\njB27gZ07LyMgoBqRkZuoWXM9d9wRzebNa1m5ciWPPPIIo0e73pKhDhu6EREvYDfQFzhC0dKCw40x\nSSXaXAeMN8YMEpFuwHRjTHd79i1xDKcV+tOn4dZbP+Gbb2qQn9+b8PB4Jk8OZMyYPg45flxcnMsv\nZF4WzW8tzV/5NmzYwM6dO7n22mtJSkqievXqzJnzPidOhLN3b3sOHWpJamp9vL230qZNKgUFa+ja\n1Q8/vwwmTJhAmzausaShI4duYoBkY8z+4gPPAwYDJYv1YOBDAGPMehEJEZEwIMqOfZ0qO9tQu/Zp\nIIzrrz/LypWXc/ToEbp2/clh53DHH/SSNL+1NH/l69q1K127dgXg7bffZurUqXTp0uXc+wUFBcya\nNY99+8I5efJqlixpzaefRiCSx0cf/UJo6M8EBh5h3LgBLFz4b7Kzkxg1aiRDhgwhMjLSqm/rguwp\n9OFAySkbD1FU/MtrE27nvk4VEGAjL68Ba9d+TbdufYA9GGPw9fWtzBhKKTfi4+PDmDG3nXv91lsh\nGAPJyTl8/nk9fv01iuXLk3jyyVPYbP+msLAujz6azsSJRwgNPYa393G6dGlEz55NmTw50PJxf2dd\njHWpyxldurSmZ8+eVsdQSrkxEWjRIoDHHy8atjl+PIgtW7bQsmUBYWE+HD1ah+XL93P8uA8//pjF\n119/zObN0UyZYv1KRfaM0XcHphpjBhS/fgIwJS+qisjbwPfGmPnFr5OAPhQN3ZS5b4ljuMZ9nkop\n5UYcNUa/AYgWkUggBRgOjCjVZikwDphf/Ish0xiTKiLH7djX7rBKKaUqrtxCb4wpFJHxwDf87xbJ\nRBG5v+ht844xZpmIDBSRPRTdXnlXWfs67btRSin1Jy7zZKxSSinncKlJzUTkQRFJFJFtIvK81Xku\nhohMEhGbiIRanaUiROTfxX/3W0RkkYi4/CorIjJARJJEZLeIPG51nooQkQgR+U5EdhT/vD9kdaaL\nISJeIrJZRJZanaWiim8D/6z4535H8TNAbkNEJorIdhH5RUQ+FhG/C7V1mUIvIrHADcBlxpjLgJes\nTVRxIhIBXAPstzrLRfgGaGuM6QAkA09anKdMxQ/jzQCuBdoCI0SklbWpKqQAeMQY0xa4AhjnZvl/\n9zCw0+oQF2k6sMwY0xpoD7jNsLKINAQeBDoZYy6naBh++IXau0yhBx4AnjfGFAAYY45bnOdivAo8\nZnWIi2GM+dYYYyt+uQ6IsDKPHc49yGeMyQd+fxjPLRhjjv4+TYgx5jRFRSbc2lQVU9yxGQjMsjpL\nRRV/Yu1ljJkDYIwpMMacsjhWRXkDQSLiA1SjaPaB83KlQt8C6C0i60TkexHpUu4eLkREbgQOGmO2\nWZ3FAe4GllsdohwXekjP7YhIE6ADsN7aJBX2e8fGHS/0RQHHRWRO8dDTOyISaHUoexljjgAvAweA\nwxTd6fjthdpX6uyVIrISCCu5iaIfkqeKs9QqniOnK7AAF5tKspz8kykatin5nkspI/8UY8yXxW2m\nAPnGmE8siFjliEh1YCHwcHHP3i2IyCAg1RizpXjY1eV+3svhA3QCxhljNorIa8ATwDPWxrKPiNSk\n6BNsJHASWCgit17o322lFnpjzDUXek9ExgCLi9ttKL6gWdsYc6LSApbjQvlFpB3QBNgqIkLRsMcm\nEYkxxqRVYsQylfX3DyAioyj6KH51pQS6NIeBxiVeRxRvcxvFH7kXAh8ZY5ZYnaeCegI3ishAIBAI\nFpEPjTF3WJzLXoco+gS+sfj1QsCdLuj3A341xqQDiMhioAdw3kLvSkM3X1BcYESkBeDrSkW+LMaY\n7caY+saYpsaYKIp+iDq6UpEvT/F00o8BNxpjcq3OY4dzD/IV320wnKIH99zJe8BOY8x0q4NUlDFm\nsjGmsTGmKUV/99+5UZHHGJMKHCyuNVA0w647XVQ+AHQXkYDizmVfyriY7EoLj8wB3hORbUAu4DY/\nNOdhcL+Psm8AfsDKop8b1hljxlob6cLc/WE8EekJjAS2icjPFP3MTDbGrLA2WZXyEPCxiPgCv1L8\noKc7MMYkiMhC4Gcgv/jPdy7UXh+YUkopD+dKQzdKKaWcQAu9Ukp5OC30Sinl4bTQK6WUh9NCr5RS\nHk4LvVJKeTgt9Eop5eG00CullIf7fws18rQIm+CuAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10fce5b38>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "params = dict()\n",
    "params['f'], params['mu1'], params['mu2'], params['sigma1'], params['sigma2'] = out.x\n",
    "xs = np.linspace(-5, 7, 200)\n",
    "plt.hist(data, histtype='step', color='k', bins=200, range=(-5, 7), normed=True)\n",
    "f = lambda xs: np.exp(logpdf(params, xs))\n",
    "plt.plot(xs, f(xs))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we're interested in the errors on our parameters, we can get them from the fisher information:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-0.399848771353 0.0104203415992\n",
      "-1.99967054954 0.00859753196565\n",
      "1.99448598671 0.0105927661859\n",
      "0.00393739343292 0.00578702226849\n",
      "0.404863029427 0.00503073097642\n"
     ]
    }
   ],
   "source": [
    "hs = hess(out.x)\n",
    "cov = np.linalg.inv(hs)\n",
    "err = np.sqrt(np.diag(cov))\n",
    "\n",
    "for x, e in zip(out.x, err):\n",
    "    print(x, e)"
   ]
  }
 ],
 "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.5.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}