Mixture_model_example.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Mixture model example."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import scipy.optimize as sco"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Polynomial model\n",
    "def line(theta, x):\n",
    "    return np.polyval(theta, x)\n",
    "\n",
    "def lnlike(pars, x, y, yerr):\n",
    "    Y, V, P, slope, intercept = pars\n",
    "    theta = [slope, intercept]\n",
    "    \n",
    "    model1 = line(theta, x)\n",
    "    model2 = Y\n",
    "    ivar1 = 1./(yerr**2)\n",
    "    ivar2 = 1./(V + yerr**2)\n",
    "    lnlike1 =  - 0.5 * ((y-model1)**2*ivar1)\n",
    "    lnlike2 =  - 0.5 * ((y-model2)**2*ivar2 - np.log(np.sqrt(ivar2)))\n",
    "    ll = np.sum(np.logaddexp(np.log(1-P) + lnlike1, np.log(P) + lnlike2))\n",
    "\n",
    "    # Uniform priors over P and variances\n",
    "    if P > 0 and P < 1 and V > 0:\n",
    "        return ll\n",
    "    else:\n",
    "        return -np.inf"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Simulate some data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAEGCAYAAABsLkJ6AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3deXxU9b3/8ddnJgsRKKLsKmVX0QKRlDqKOhrXatUqat2LvSK2Uttba7H+er3e1uKttnpvtSpt9WqlVitqrddbVCQqJS5hCQIuIAiK7EIICEkm5/v7YzLDTJiwZJkzy/v5ePAgOTlnzneG8P2c7+e7mXMOERHJPwG/CyAiIv5QABARyVMKACIieUoBQEQkTykAiIjkqQK/C7A/evTo4QYMGOB3MUREssrcuXM3Oud6Nj+eVQFgwIABVFVV+V0MEZGsYmYrUx1XCkhEJE8pAIiI5CkFABGRPKUAICKSpxQARETylAKAiEieUgAQEclTCgAiInlKASADhcNhwuGw38UQkRynAJCBampqWLVqFZWVlX4XRURymAJAhqmsrGThwoWsWLGC8vJyBQER6TAKABmmoqICz/MAqK+vp6Kiwt8CiWSg0tJSBg0apAekNlIAyDDhcJhAIPrPUlRUpL4AkWbUSm4/CgAZJhQKMWLECAYOHMjMmTMJhUJ+F0kko6iV3H6yajnofDF//ny/iyCSsWKtZM/z1EpuI7UARCSrqJXcftQCEJGso1Zy+1ALQEQkTykAiIjkKQUAEZEMks6lYBQAREQySDqXgvE1AJjZx2b2rpktMLMqP8siIuK3dE9yy4QWwMnOuVHOuTK/CyIikk7N0z0tTXJbuWk7zrl2v38mBAARkbzUPN3TfCmYPocfw4THqgjfXcHrSze2+/39ngfggJfMzAEPOeemNj/BzCYAEwD69++f5uJJLok9aWnpAMkEpaWlLFiwAIDy8vL4pLYRI0ayoT7I4LP+hdvfjtCt5HNuOHkIw/t+qd3L4HcAON4595mZ9QJeNrP3nXOvJ57QFBSmApSVlbV/G0jyRk1NDTU1NVRWVmr2qPhu3bp18a/r6+t5ZearLA/048DLf8Pmjdtp7F7Cv48dyMVfPYwDijqmqvY1ADjnPmv6e72ZPQuMAV7f81Ui+y/WueZ5XtLTlkg6NG99VlZWsnbt2vjPPQePr+xC3bOLGHFoN+67rJQzj+pDQbBjs/S+9QGYWWcz6xr7GjgdWORXeSS3aQVJ8VNNTQ1vvvkmpaWlQPT3MbFTt7D3YIb06sKfr/0af/ve8Zwzol+HV/7gbydwb2C2mVUDbwP/65z7h4/lkRymfRbEL7HWZ11dHdXV1fzbnfdSVdcXKyiKn1O/5kNm3TMJW78UM0vbZDDfAoBzbrlzbmTTn6Occ3f4VRbJfVpBUjpaql3KwuEwl19+ebz16Zzj57fexIJPNtN39GlJ19fV1cVbpumaDOZ3J7BI2mgFSekozfuYDj30UCKRCAWFRdTaAWABcNEgYM5x5YDtFA85np++9b/x1wgGg4TD4bT2VykAiIi0UWIfU11dHcuWLUvI8RsEAtFB70CnTsWccWo5QHxjGzPjvvvuIxQKMWXKlN36qxQARER8kGr+SGlpKTU1NUybNo1QKJS0SxkWwHmNCa/gCBr06tuXTp06xa8BGDFiRNLrxO6Xrh3PFABERPag+fyRVCmarv2H0//r17Nl61YKDujGppcfwjU2xF+jqKiI6dOn7/YknyotGQqFmD17NhUVFYTD4Q7tr1IAEBFpQarKPindU1/P9Xf9iS3DzqZz6Tl8Z0x/xo8dyDknz2bx4sV4nkevXr1SVv57EgqF0jJQQQFARKQFiZX9jh07qKioYOwJJ2LBApzn4SxIfc8j+MmZR3DZ1/rTraQQgN/97neMHTsWz/PYsmWLn29hjxQARERSKC0tZd26dZhZvEP3d8/M4nnG0PvSKdR/9BY//JfL+fFV36C4IJh0baqJh5k49FgBQESkmcTUT2FhIQ0N0Xz+p1Wv0OfwUv78/35E+PAfYWYpr09nR25baDloEclriRO4YjNwE5/gY5V/lGPeE3fTafNHLVb+kD0TD9UCEJG8lWoC145gZ55/vxYrKMI1RgAHCev2eJ63TymdbJh4qAAgInkrcZmGnXV1LP1oOTjH6mXvUdJ3KPVrl3HKySfx0ksvxa+JzdjNBQoAIpKT9rYBUDgcZs26jWAWfcjHwPMAh3kN7Ph0Cc45Zs2alXTdtddem7Epnf2lPgARyUl7WlBtQ20dy3YcwM6dX0TTOwHj4m9PJBCI5vWDwWB85E9jY2M8319SUsJVV12VvjfRwRQARCTnxHL7K1asoLy8PB4Elq3fxuTpCzlm0v2srnoJYks2eB4vP/tnBg8ezMCBA7nvvvviy4cXFxczZMiQjO/QbQ2lgEQkK8VSPLGlGhLX02k+Dv/xZ1/kjx8W8sp76yguCDAkspJPk9brgc2bN7Nz5854Jf+Vr3wlLcsx+EkBQESyUk1NDevWrWPt2rU45xg7diyzZ89utjibwwUKeGZNN/p0+ZzvnzKEq44bwNJ3D2TsE/fhed6uRdxInrSVruUY/KQAICJZJ3H4ZozneTz22GNcf8P3Wd9QxEHhq4lEGvny0V/lR1ecw7jRh1FSFJ2x26NpnH5NTQ2TJ0/m+uuvz/hJWx1BAUBEMlbzZZdjElM8iX7/hz/Q6ADPI1Awl189/Fd+cNk5BAO7T9pKHKf/wAMPpLzP/pYr2/geAMwsCFQBq51z5/hdHhHJDC3tjBUOh6mpqUl5TWOkkdjOK+ZFqP9kEcHAN/Z6r/2ZtJXOHbs6WiaMAroReM/vQohIZki1HEMsNw+7cv9xgSAEC8GMoqKi+JDNjkrntFSubORrADCzQ4GzgT/4WQ4RyRyx8fsHH3xwfChmrDL/5z/nUL3wXdasWRM/34A+ffowcMAAKipmMXLkyA4dshnrYE4sV7Yyl7DGRdpvbvY0MAXoCtyUKgVkZhOACQD9+/cfvXLlytbdbNWbUFfb+sKKSIcbdtp4lq5cDUBJp2IO7d2DSGMjD//qJ6ztegS3/vYJlr/8eHyD9dh5Mx/9T0Klw9NWztLzrqemdjvTfj05ffct7gr9j23VpWY21zlX1vy4b30AZnYOsN45N9fMwi2d55ybCkwFKCsra320qquFzj1afbmIdKzKqmqWrVod/76uvoFvXXQ+vcdezE3V29n4RSP9hpWy4pVp8bXZAoEA995+E6GxJ6a1rPNf+Wta7wfA9o3t/pJ+poCOB841s4+BvwCnmNnjPpZHRNKo9IxLGRQ6h8qqagAqKueSlJAw47GNA7mrspajehYy7ZsH88/JJzJy+BD69u6BmeF5Hj+47e74a8j+8S0AOOducc4d6pwbAHwLeNU5d4Vf5RGR9Kmsqmbhkg9Zseozyi+ZSGVVNeHQ6OjCbBgEghx02kQuOKmUf1zWk0fPO5jjDyvGzJg/4wkmjf9WfK2e+oYIFZVz/X1DWcr3YaAikn8qKufiedEKfMfOOs78zmTGXHULXcdcwBfv/5PjTgrz5K++Td+uwZTXh0OjCQQMz3MUFRZEg4fst0wYBopzrkJzAERyS/MUT6JoBb6r+tm6cT0z772J2rem01izljee/wt//9tzLb52qGwkI4YPY2D/fsx88kFCZSM75D3kuowIACKSWxJTPGO/eU1SENha57HQBtNn3L/RaUAp0YGcDudFkl5j+osz93iP+TOeYHnlC6r820ApIBFpd4kpHs9zPPb0Cww4/GgeXrCNJxZ9wbYGx0knHMfx5X254brrACgoCBKJ7Fqh88Kvl/tS9nyiACAibRIedy0AFU//ftex0GjMLN5R+/snnufvHEvxIUfy9aElXHdMZ47uVURl1a5hn4FAgKJCo9Fz/O6Xk5lwxYXpfSN5SAFARNqkpnYbNVtrqayqjqdjjh09gqJuPanbsh6I7qp1RP0HPHL1SRz2pV3VTuLoncZGj5//+HpumXRNet9AHlMfgIi0WukZl7Jg0Qfx4ZxvvF3N3z74gsOnzKPz1y6KrtMDlBQX8vPLjk+q/GHXaB5Ao3l8oBaAiOyTVKmedRs2xb+uq2/git+9hh3Tk+1LZvP5rMeIrcx57+03peysjY3mqdlay7Tf3qEO3TRTABCRfdI81VNZVc3a9bsCgGcBBh81ivMHrObG/3w06dr5i95v8XXnz3iiw8ose6YUkIjsVfOZu0+9Oo/zb5+GY9dGKxee/w1e/fGJbF/5ro8llf2hACAiuwmPuzae8oHkYZ076xu47vevEzl4SPznJZ2K+dFV50avbRoBFFNQEOSqcZrnmYkUAERkNzW121i1eg2VVdU0eo7Cw76CFRSCBSBQwOnDe7Jl5lRwHmbJOf5Q2UhGHjWMwsIC+vbuwevT/6DcfoZSH4BIHmjegZty7H7TsSmTb2Dhkg/xPEf4ookcOX4KWw4cRr9vTsZbMoM//8dEKt+Zx1NNu2IFAkE2bU7eolF5/eygACCSB6qqlxBpbIx34KYaux879tBf/k5sv/X6hgYaPn2XBy/7GqcNOp9g4JsAFBeYFmPLAQoAIjmusqqa7V/sAOC488Zz83evjj/hl18ykZlPPsh3b72TBYs/BAcrVj8PwQJobMBw3HhiP84cUpL0mqGykcx+9mEqKucSDo1WiidLKQCI5LjHnn4h6fu7HngsaS39aS+9zceRA5uG7DvMeZwaPoFXZr2Gc44f3HY3XzliyG6VfKhspCr+LKdOYJEctzZhshbArn3ADQIFPLEswJal84hN2iosCDL4kB7x3bm04UruUgAQyXF9eh6827EDBh3DgSdeyfBrpnBynwZobIj/7JpLzuWqcedoiYY8oBSQSI5INbIHoPToI5JPtAClZ1/JTy4Yw1lDOvHOvO48+6ff43mOkk7FXDXuHC3RkCcUAERyRKqRPZ/URHhy7jpim64AnH/+eTzzkxPjk7Vaquw1lDP3KQCIZKnEJ/7KqmoWLPoAgPJLJvLQA/cx5fVNfNFzOJHiw7FgENcYoaRTMTd/+9ykmbqgyj5f+dYHYGadzOxtM6s2s8VmdrtfZRHJJgceeSIHHnli0mzdxE7aHXUNfP+RN9i8s5GdM37Ng2cfxMgjB2v/XNmNny2AOuAU59w2MysEZpvZ/znn3vSxTCIZr66unvqGSNIT/2nnjtt1gvM4tmeEl578GZ7ncfG331TFLyn51gJwUduavi1s+uP2cIlI3pv6+HR21tXjxabqEn3if+ntJfHvA4EA3sbl8XM0jFNa4uswUDMLmtkCYD3wsnPurRTnTDCzKjOr2rBhQ/oLKZJBpr84c7djZgEODHwR/764qJALv16uYZwdpPSMSxkUOofKqmq/i9JmvgYA51yjc24UcCgwxsyOTnHOVOdcmXOurGfPnukvpEiGWLyhgciXQ7sdd16EdSs+BIivzDnhigsZMXyY8v7trPm+CNkeBDJiFJBzbouZVQBnAot8Lo5Im7U0Jn9/Oed4Y1UdU+dtY/Yn9XTufyr9xrzP5wtnsXPnzqTzIHllTo3saX+J+yLEUmvZHFz9HAXU08wObPq6BDgVaHnfOJEsEhuhU3rGpUkbq+yrhkbHM+99wVlPbOCqv33OB5si3HxcV+aM783qZ+/g1ScfoNlITkDpno6Wa5vY+9kC6As8amZBooHoKefcC3u5RiTtSs+4dL9mxMbSBJ7nMIM+vXokTc5K1LylcOIlN1DbZzSdRp7Fmm0eQw8q4FenHsh5w0ooLthV44fKRjJkYH+WLl8FEL/P9Kl3ZfUTaabLtRnSvgUA59xCoNSv+4vsi8TKPLZ08t7+0yemCZyDNes2tnhtbPbuiLO/A0eezobDvs7O1R9yzOcf8vCFYwgPKCaQ6lEf+PCN5/Y7OEnb5VJqLSP6AEQyVWtyvrE0Qey6xGtvufO++LGa2m0sXLI0et7q9ZRsK2bH8nfAa6Ty7ScpCT9IYOCe75VLlZGkn1YDFdmD1uR8Y5ulTLzywviSC7Fra2q3sXL1Gj6u8Vi62TUFCQeNEXYsexO8RkBj9yU91AIQ2YPW5HwTO30PKOlErx7dmfbbO4g4qF6yDOd5wGd7fI1c6GCUzKcAILIX+5tmieX1CwoK6NWjO3+85xe8XzCY2//r4abhmrtSQ4FAgLIRR/JO9RKccxQXFTH+km/El2QW6UgKACLtKLHTOKb8kol0L7+Wxk+X7XZ+cVEh995+E9+99U515kraKQCItFHiUM7LJ92aVPkDuEg9n8+4HwAzIxAMAtCrR/f4sE115oof1Aks0grhcdfGK/5Yx+4f/1HFjsNCEAi2eJ0BntdIY2MjW2pq01RakdTUAhBphVie/8E/TY927DqPf5nwPcCB14iZEQwEaPS8hE3YATNcDi0lINlNLQCR/RTN8y9lxarPuP6nd0ZH9TgHXgQaI0A01dPoNSZV/mZw08Qrc2opAcluagGI7IfPdzRyx1NziKf5m8btAxQXFlLfUI9zEAwEaIjsWrO/e7eu/O9j/02obCTnnxGmonIu4dBoPf2Lr9QCENkHK7dE+NmsLRz3yHrm2jDMkv/rnH9mmFl/fYiRRx3OwP79uO+On8Sf9Es6Fccrf4jOLbhl0jWq/MV3agHkoPZaijjXpfqcmh+bv7aeqXO38Y+PdlIYgPOPOIBrLzmRk14+kDXrNsave+GVN7j5+quTRvN85YghetKXjKYAkINiHZQtrUApUak+p5rabWzZuo3/fv4dZtcP5O3P6vlSsTFxdBfGj+pMr87RET7Tp97F8eePJ5bi9zy3W4duqGykPn/JaAoAOaY1q1fmo1SfU10jVL8XXarhB5MmcfQ1d/Kzc8u45KgD6FKUnPJpvhxzcVGhOnQl66gPIMekWr1SohLH7jf/nKb8dQ6X318RHaLpHOZFuPjAj/hOaZfdKv+YD994jlFHH65tFyVrqQWQYxKXItYww2SJKZ/o5xTAcw7PgrzDMEaNKmLtrMfxXPSJvvz4sr2+pmbwSjZTAMgxubZjUXtJTPmccvFEzv7Xu+h16RR2rlrE6Sd8ldsu/hrDexZSOmOIPjvJGwoAPunInZz0VLr75zursiqe8tlZ38Abb87l+9d9m/GjTqFf111LN+izk3yiPgAfxJ5GV6z6jPJLJlJZVe13kXJK88/3l9Pf5vnawVhBEViAwsJCpn3vJG49oVtS5S+Sb3wLAGZ2mJnNMrP3zGyxmd3oV1nSTR217SOxUzdR4ue7o66Bu56eQ7cvD2fYcacz4LB+vPbXBzk1NCrdxRXJOH6mgCLAj5xz88ysKzDXzF52zi3xsUxpoY7a9pFqHP+a2kZWdj4CKyjCNUYIFhTy88uP53tn98Qu/w+fSyySWSxppcJUJ5jdAExzzm3u0IKY/Q24zzn3ckvnlA3o5qpuG9u6G+z4HIKFrSxd+3v0lUVUrW7gshMGExr0Jb+Lk3Uql2/l5N9U09AIxQUBXph0NP17dGbTDg8HLF6xgfnLN3H52AH6fCU3NDZAyUGtutSueXGuc263YW370gLoA7xjZvOAh4EZbm9RY38LZzYAKAXeSvGzCcAEgBH9Strztr66+tSjudrvQmSJVMGy4sMaGhrBc1Af8fjr/M8Zf0IJvTsH6dMlwLGHHAJjD/G55CKZba8tAAAzM+B0YDxQBjwF/NE591GbC2DWBXgNuMM598yezi0rK3NVVVWtu9HSl6Fzj9ZdK+1ub+sVxX4+ZfINjP3mNXieo6RTMTOffJCy0hH85m/vcMsPJ+EaG7FgAT/8xT387KKvcWAnjWuQHLV9Iww9rVWXmlmrWwA455yZrQXWEs3ddweebsrZ39yqEkULVQhMJ5pi2mPlL7llb+sVxX7+2NMvJHWY3z29krXv9mF17aH0O/9m3JKXmPaL7xL+msbsi+yvvQYAM/s+cDWwEfgD8GPnXINF18NdCrQqADS1Kv4IvOec+01rXkOy097WK0r8+UOPNz0XmOFZkDneUE7oGuT2k7pxysALCNiFPr2L/KHVZXPXvrQAegAXOOdWJh50znlmdk4b7n08cCXwrpktaDr2U+fci214Tclw4XHXsmr1mt2GwSYGgMRhnM45ivoMpajPYE494+v87KIxHNO3yJey5yutLpu79powdc79W/PKP+Fn77X2xs652c45c86NcM6Navqjyj/H1dRuY2ddPdEGYHRbxKdeeDn+lOmco+eQERAsILqFOtSvXUbk/VncfFzXjK38W5qTkO00aTG3qcdM0qayqpoFiz5gzbqNFBYWMPHKC5n55IMArFy9ll8/9w7nP7WRXyztR/eTribQqUvTlY6GDJ8wV1O7jVWr1+RcBalJi7lNAUDSJrHyiEQa6X9IX+oijur3lvHxJ2v48Y2TWPXhYq7q/Qk1FY/g7ayNn5/JE+Zy+Sk5NmkRMvvfQFpHAUDS5uDu3eJfe57HvC3FXHb/a01r8Hu4SANbXvw1JesX43nRDdUDgQCnn3RsRq+3n8tPybHVZbXnQW7SaqDSoRJHkMxf9H7Sz/7x5hLCp5/FixXRNfjNHHV1dRzcvVt8qYziokL+/V+vy+iKJ9eX9tAKqblLLQDpULHc+CMz5vLayrqkn31jWAl/v3EsI4YPoW/vHoCxZt1GfnDb3QwecFjWPHXqKVmylVoA0ibNR77ExoqHx11LTe12qpcswzmP70z4HgeFvw0WAOdRVFTIpMu/AUSfMKf89mF+eud9QDSNMv7ic7ll0jVpfS9toadkyUYKANImsTHiBQUFRCIRKquqKR01gjWFh/LJ1o9wngc4XKSexqonGTrwUCKRyG4b4aQzjaKJTSJRCgCyXxIrz8QZuzEnXTSRL51wJZ/P+T9cpD7p2i1baqjbuTNlmiSdW1lqYpNIlAKA7JfEyjNx9EtMQ0MD26tnJFX+ZkZs0cFUM39j9jeN0pon+b0tQyGST9QJ3EGyfWZoqvI3H+++s6ALFiyEplm9AAVBo27jJ0nXJa44WxAMtlt6pzWTr3J5yKbI/lIA6CDZPjM0VfmTtlrcWc+vHp9B91OuYcyp34if45xjTwuMj7/kG+3yxN3ayVea2CSyiwJAB8j2maGpyl8XcTT2PRorKCS6Ro9j58fz2fzKVPoFt8av3dP2EsVFRVw1ri3rB+7S2id5DdkU2UV9AB0gVeWU6RVNYj69efnvfHoOa6p7s377ofQ77ya2vv4otRs+A6JP/H9/+Q3MopV/cVEh995+E/MXvc+8d9/n7QWL4/dor6d/aNuoIQ3ZFIlSAOgA2TgzNLFzN1r+AJ5zeBbkbTeM8oMKuPvULpzQfxxvzh3G8eePjz/tO+fo06sHnYqLkkbwVFZVJ+3m1V5P/5DeUUMiuWqftoTMFNm0JWTpGZdmTeVUesalLFj0AQCdiovpU/5tdhw4iLpPF1G/5kPqP1nM/bffyIQrdm2+MuyE81m6fBVAfKvGVO8zmz4HyS0597vn15aQsv+yIc0Q+w+ytXZ7/NjO+gY2rl5B59UfcGnZYTzw+lsAXPeTOwDiQeDDN56LDwUNh0a3+B8sGz4HyT0a7rtvFABywL6Mh29+TqpJXAA4j23VL7ENeOjd5DEC01+cmdQKCJWNJFQ2MveetCTrZWM/nB80CigH7MuQ08RzttV73DW9Es9Zi+dD8vh9gAu/Xr7bOdk+4klyk4b77hsFgCy3LxVw4jknXTSRkf9RQaU3NLow2x50Ki6iqLCAYDDIQ/95a9LTf4wmVkkm0nDffeNrCsjMHgbOAdY75472syzZal+auk/PfJum/VVoaGigduHLbF+5CLxI/JyDu3fj8y01OBdduqFPr4OZPvWuvf7HycYRT5If1P+0d373AfwPcB/wmM/lyGipcvyxY1Mm3xA/FjCLV8AnjbuWnd0GMPy8ifzj80FYQSGusYFOxYUUr1nAho3rdl0XCPCj664kHBq9107d5jQcUyR7+RoAnHOvm9kAP8uQDVKtXhk79tyMivh5DZEIz/yjgo1fGsbyvuVs27iWLVULueXCMfzx3SOo3byRyd8bz8TJdyS9fiBg8Uq/NRW4nrREspPfLYC9MrMJwASA/v37+1ya9Es1nO27t94ZH7d/94N/Sjr/7gce5eHlXdg88/e4SAN1bz/JV899kO+/+D8ATPntw7st13D/HZNz9slda/+LtCzjO4Gdc1Odc2XOubKePXv6XZy0S5XjX7dhU/znqSby9Vj1atNyzG63jtnE0RFmtNi5myuyfVE+kY6U8QEgV+3rctGJnapFhQUc3L0ba9dv3HVCoICCngOTrhn+5d7xr4PBQNJrJI6O+Odzj+R05a8hqiJ7pgDgk319Mg2VjWTooP5079aVe2+/ieqVm3DExu8bI086kz/fc2v8/OKiIs46+XisaY1+Y/ex/vNnPMHyyhdyNu0ToyGq2a30jEsZFDpHgbsD+RoAzOwJoBI43Mw+NbPv+FmedNnTk2nzlkFlVTUfffwJm2tq+e7/u5u/fBSMb8LSqVMRv7vxAi46+RhGHX04A/v3Y9ZfH2LT5pp4aijS2Ji3FZ8mA2Uvtd7Sw+9RQJf6eX+/7GnsfuKIn9JRI7jn2UpiqzU0Rho4pU8D333yQd58Z17ScM3mI3FyaWx+aztyNUQ1e2kph/TI+FFAuailyVOJI35OumgiA6/8JVsjQ7FgAS7SgOFx2vCenHzsKE4+dlSLr59rFV9bNnHXENXspAmG6aE+gDRJTO2EykZS0qkTxcVFSdPUn531Tvypp6G+nq4bF/PspLE8eMfN8Q1XfnDb3fvUHM6VPL9SAflJSzmkh1oAbbSv6Yma2m28t3QFXYYeT9nI4Qwd1J+arbUALFhbz9R523hySWKHrWPCcX059tBiZm2uiY/dz7fmcDamAjT3oH2o9dbx1AJoo30ZzRN7iq2rq2f7FzuYv+h9qhd/wIpVn3HCuImcdc8bvLGqDj76Z/yaQCDAps01QH53Zmbje9fcA8kWCgBt0FJ6ovlInsSnWICttdvjT/SNkQaOCyzl3lHr2Lx8V4VRUBCMV3b53BzOtveulJU/9nVejSRTCqgNUqUnbrnzPqqql9CrR/d4p2U4NBozSzlrt6S4kJsuDFFROTfp59dccpGd5qsAAAt3SURBVG5SZZfPzeFseu/ZmLLKBW0ZKJDP1AJog1Tpic/WbWD7FzviT4B/e20eL30xgO6nXrfb+vt9e/eIP9UmvlZ7b6Au6ZONKatsp1ZX6ykAtEGobCRFhYUEg0Huvf0mAJatWBX/+Y66BsY/9Dp/Wridqy+7iEGHD4//rLioKGm9/WxLdUhq+ndMP834bj2lgFohcS3++oYGPM9x/S2/5NrLL0headOMi8rHcNvFvenXNUjPD8L89M5FwK4Zukrz5J58/Hf0c19ozRloPbUAEuzr2iOxUR6PPf1C/MnD8xzPzPsMKygCDCzAkIFf5oP/e4R+XYOA0gOSm/xOwajV1XoKAE329Ze4sqqaBYuiQzgfefL5pLz+pg+qOOzorzKg/yHMee6PdCkpTBoOGCobyexnH+aXk2/QL6rkjExIweTKxMd0y9sA0PxpP9UvcaqhZYm/3HUNEYr6DIl/b85j4pmlrKh8HiBlQAmVjeSWSdfoF1Vyhlq22SsvA0Cqp/1Uv8TNJ/S8v7GBmWsKd72Qc1x+xtdS/vJnwlORSDooBZO98rITOLFy3rGzjseefoEH7ryVEcOH8d7SFRw+ZEDStosnXzyRUyb9iiXBwXyx8vP46wQCAYb26ZZy4TV1TEk+yceO71yQlwGg+cSsh598nqvGncP8GU9QMuhY3n1vGZ0PKImfX1ffwJuVlfQ/5GO6bZ7PpqaKvbiokHBoNLdMuma3e+TaipwiknvyMgCEykbSp9fBrFkX3VqxsdGjonIu776/jJ119QBsrd226wLnUbzsVb5Y04lIJMLgAYcRiUT2WrHrqUhEMlle9gEATJ96V1LufkRpKXf/+aUWz1+7bj1Ll69ixarP+PSzdXqqF5Gsl7cBIJai6T9oKBf8+G5+uKA3m/qM2adr1akrIm2RKYvX+b0n8Jlm9oGZLTOzyR19v9jQzzlV1by9uo7RN9yHXXQP7zQO5OKjDqDqN1fHN2C//IKzsN33Uwc01E1E2iZTlgz3rQ/AzILA/cBpwKfAO2b2vHNuSUfcr3L+kvionhPGTaTXJb+gz+CjuHFMF575r58xe9YX/OLp3/PhG8/Fr1n84XJqttYy+Xvj2bS5hkeeen6fcv8iIi1J3Pq1/JKJvg6d9bMTeAywzDm3HMDM/gKcB7R7ANhR38g9f18AGODwIvWcULCUR8efTElhgMduW59yKdnmnbipRvuIiOyPTFoy3M8U0CHAJwnff9p0rN399Nl3mbmhM9C0UptznHpkT0oKA76vY5JPMiXvKeKnTJo57WcLIFWGfbcdU8xsAjABoH///q260YQTB1FX0cADTd8HAsamdZ/B9o1UvPZ6QjRuoOK11wkd2SFxKO/VbNlCTe12Kme/Tqg0ujR2+IroMtoVj9/tZ9FE0iZ05CGMOHwQNbXbmfbrydH6ZvvGvV9Y3LX9C+Oc8+UPEAJmJHx/C3DLnq4ZPXq0a605c+a4QCDgAFdSUuLmzJmzx+PSvlr6nEeNGuUGDhyoz12kAwFVLkWd6mcL4B1gqJkNBFYD3wIu66ibhUIhRowYQU1NDdOmTSMUCsWPz549m4qKCsLhcPy4tK+Kigo8zwOgvr6eiooKABYuXIjneZSXlzNz5kx9/iJp5FsAcM5FzOwGYAYQBB52zi3uyHvOnz8/5fFQKKSKp4OFw2ECgQCe51FUVEQ4HE4ZFPTvIJI+vi4F4Zx7EXjRzzJIerTUAmseFEQkffJyLSDxR/MWWEtBQUTSQwFAfNVSWk5EOl7ergUkIpLvFABERPKUAoCISJ5SABARyVMKACIieUoBQEQkTykAiIjkKQUAEZE8pQAgIpKnFAAkrrS0lEGDBlFZWel3UVIKh8NaL0ikHSkACACVlZUsXLiQFStWUF5enpFBoKamhlWrVmVk2USykQKAAC2v158psiFAiWQbBQABdq3XD2Tk0syZHqBEspECgAC7lmYeOHBgRu7MlekBSiQbaTloicvkpZm1d4BI+1MAkKyRyQFKJBspBSQikqd8CQBmdpGZLTYzz8zK/CiDiEh7ydY5Kn61ABYBFwCv+3R/EZF2k61zVHwJAM6595xzH/hxbxGR9pTNc1Qyvg/AzCaYWZWZVW3YsMHv4oiIJMnmOSodFgDM7BUzW5Tiz3n78zrOuanOuTLnXFnPnj07qrgiIq2SzXNUOmwYqHPu1I56bRGRTJHNc1Q0D0BEpI2ydY6KX8NAv2lmnwIh4H/NbIYf5RARyWe+tACcc88Cz/pxbxERicr4UUAiItki2yaEKQCIiLSTbJsQpgAgItIOsnFCmAKAiEg7yMYJYQoA4qtsy5mKtCQbJ4QpAIivsi1nKtKSTN9VLxVNBBPfxHKmnudRXl6eNf9pRFqSbRPC1AIQ32RjzlQklygAiG+yMWcqkksUAMQ32ZgzFckl6gMQX2VbzlQkl6gFICKSpxQARETylAKAiEieUgAQEclTCgAiInlKAUBEJE8pAIiI5CkFABGRPKUAICKSp8w553cZ9pmZbQBWtvLyHsDGdixONtB7zg96z/mhLe/5y865ns0PZlUAaAszq3LOlfldjnTSe84Pes/5oSPes1JAIiJ5SgFARCRP5VMAmOp3AXyg95wf9J7zQ7u/57zpAxARkWT51AIQEZEECgAiInkqLwKAmZ1pZh+Y2TIzm+x3eTqamR1mZrPM7D0zW2xmN/pdpnQws6CZzTezF/wuSzqY2YFm9rSZvd/0b53ze2qa2Q+bfqcXmdkTZtbJ7zK1NzN72MzWm9mihGMHmdnLZra06e/u7XGvnA8AZhYE7gfOAoYDl5rZcH9L1eEiwI+cc0cCxwLfy4P3DHAj8J7fhUij/wL+4Zw7AhhJjr93MzsE+D5Q5pw7GggC3/K3VB3if4Azmx2bDMx0zg0FZjZ932Y5HwCAMcAy59xy51w98BfgPJ/L1KGcc2ucc/Oavq4lWjEc4m+pOpaZHQqcDfzB77Kkg5l9CTgR+COAc67eObfF31KlRQFQYmYFwAHAZz6Xp905514HPm92+Dzg0aavHwXOb4975UMAOAT4JOH7T8nxyjCRmQ0ASoG3/C1Jh7sXuBnw/C5ImgwCNgCPNKW9/mBmnf0uVEdyzq0G7gZWAWuAGufcS/6WKm16O+fWQPQBD+jVHi+aDwHAUhzLi7GvZtYFmA78wDm31e/ydBQzOwdY75yb63dZ0qgAOAZ4wDlXCmynndICmaop730eMBDoB3Q2syv8LVV2y4cA8ClwWML3h5KDzcbmzKyQaOU/zTn3jN/l6WDHA+ea2cdEU3ynmNnj/hapw30KfOqci7XsniYaEHLZqcAK59wG51wD8AxwnM9lSpd1ZtYXoOnv9e3xovkQAN4BhprZQDMrItpp9LzPZepQZmZEc8PvOed+43d5Oppz7hbn3KHOuQFE/31fdc7l9JOhc24t8ImZHd50qBxY4mOR0mEVcKyZHdD0O15Ojnd8J3geuLrp66uBv7XHixa0x4tkMudcxMxuAGYQHTXwsHNusc/F6mjHA1cC75rZgqZjP3XOvehjmaT9TQKmNT3YLAfG+1yeDuWce8vMngbmER3pNp8cXBLCzJ4AwkAPM/sUuA24E3jKzL5DNBBe1C730lIQIiL5KR9SQCIikoICgIhInlIAEBHJUwoAIiJ5SgFARCRPKQCIiOQpBQARkTylACDSBmb2VTNbaGadzKxz01r1R/tdLpF9oYlgIm1kZr8AOgElRNfnmeJzkUT2iQKASBs1LcXwDrATOM451+hzkUT2iVJAIm13ENAF6Eq0JSCSFdQCEGkjM3ue6DLUA4G+zrkbfC6SyD7J+dVARTqSmV0FRJxzf27af3qOmZ3inHvV77KJ7I1aACIieUp9ACIieUoBQEQkTykAiIjkKQUAEZE8pQAgIpKnFABERPKUAoCISJ76/+lDgpWunaV5AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "np.random.seed(42)\n",
    "\n",
    "# Parameters\n",
    "V = 2**2  # The variance of the outlier population\n",
    "Y = 2  # The mean of the outlier Gaussian distribution\n",
    "P = .15  # The probability that a star is an outlier (15% of stars will be outliers).\n",
    "slope = .5  # Line parameter\n",
    "intercept = 0.  # Line parameter\n",
    "\n",
    "nstars = 200\n",
    "x = np.random.uniform(0, 10, nstars)\n",
    "y_uncert = .1\n",
    "y = intercept + slope * x + np.random.randn(nstars)*y_uncert\n",
    "yerr = np.ones_like(y)*y_uncert\n",
    "\n",
    "# Add outliers.\n",
    "nout = int(nstars * P)\n",
    "outliers = np.random.choice(np.arange(nstars), size=nout)\n",
    "y[outliers] = np.random.randn(nout)*np.sqrt(V) + Y\n",
    "\n",
    "xs = np.linspace(0, 10, 100)\n",
    "plt.plot(xs, intercept + slope*xs)\n",
    "plt.axhline(Y, color=\"C1\")\n",
    "plt.fill_between(xs, Y - np.sqrt(V), Y + np.sqrt(V), alpha=.2, color=\"C1\")\n",
    "plt.errorbar(x, y, yerr=yerr, fmt=\"k.\");\n",
    "plt.xlabel(\"x\")\n",
    "plt.ylabel(\"y\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Test the likelihood function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-150.38688879464635\n"
     ]
    }
   ],
   "source": [
    "p_true = [Y, V, P, slope, intercept]\n",
    "print(lnlike(p_true, x, y, yerr))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Optimize the model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-450.8723236497227\n",
      "[ 2.07600245e+00  5.03943111e+00  6.31049597e-01  5.01141131e-01\n",
      " -2.08063826e-03]\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/rangus/Applications/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:15: RuntimeWarning: invalid value encountered in log\n",
      "  from ipykernel import kernelapp as app\n",
      "/Users/rangus/Applications/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:15: RuntimeWarning: invalid value encountered in log\n",
      "  from ipykernel import kernelapp as app\n",
      "/Users/rangus/Applications/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:15: RuntimeWarning: invalid value encountered in log\n",
      "  from ipykernel import kernelapp as app\n",
      "/Users/rangus/Applications/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:15: RuntimeWarning: invalid value encountered in log\n",
      "  from ipykernel import kernelapp as app\n"
     ]
    }
   ],
   "source": [
    "p_init = p_true + np.random.randn(len(p_true))*.1\n",
    "\n",
    "# Calculate maximum-likelihood values\n",
    "nll = lambda *args: -lnlike(*args)\n",
    "print(lnlike(p_init, x, y, yerr))\n",
    "result = sco.minimize(nll, p_init, args=(x, y, yerr))\n",
    "print(result.x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAEGCAYAAABsLkJ6AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3deXyU5b338c9vsiIgLqyyCCi4E1IiOmo1SBV3bUVxV+wj4kL1ObVW2z61eup2tEfrLrVarZRqwVqrHpeDRkWjEoQgi8qmyA6KIUGyTO7r+WMywyRMIGSZe5bv+/XyZTJzL9ck4f5d1+/azDmHiIhknoDfBRAREX8oAIiIZCgFABGRDKUAICKSoRQAREQyVLbfBdgV3bt3dwMHDvS7GCIiKWX27NkbnXM9mr6eUgFg4MCBlJWV+V0MEZGUYmZfxXtdKSARkQylACAikqEUAEREMpQCgIhIhlIAEBHJUAoAIiIZSgFARCRDKQCIiGQoBYAkVFxcTHFxsd/FEJE0pwCQhCoqKlixYgWlpaV+F0VE0pgCQJIpLS1l3rx5LF++nNGjRysIiEiHUQBIMiUlJXieB0BtbS0lJSX+FkgkCRUWFjJ48GBVkNpIASDJFBcXEwiEfy25ubnqCxBpQq3k9qMAkGSCwSDDhg1j0KBBzJgxg2Aw6HeRRJKKWsntJ6WWg84Uc+bM8bsIIkkr0kr2PE+t5DZSC0BEUopaye1HLQARSTlqJbcPtQBERDKUAoCISIZSABARSSKJXApGAUBEJIkkcikYXwOAmX1pZp+a2VwzK/OzLCIifkv0JLdkaAGMcs4Nd84V+V0QEZFEapruiTfJrbqunudnfU11XX2731/DQEVEfFJRUUFFRQWlpaUEg8HtJrlt6LIfx9z9FhuratktL4vThu3Trvf3OwA44A0zc8DjzrnJTQ8wswnABIABAwYkuHiSTiI1LS0dIMmgsLCQuXPnAjB69OjopLZDioKsr8snv+BUXljVmeIDujHh2MEEB+/d7mXwOwAc7ZxbbWY9gTfN7DPn3LuxBzQEhckARUVFzo9CSnpoWtsS8dO6deuiX9fW1vLcS6/xty/z2TL6V3Q1OL1gHyYcO5gDe+/eYWXwNQA451Y3/H+9mf0TGAm8u+OzRHZdpHPN87xGtS2RRGja+iwtLWXt2rXR9z1nPLeyK3vnruenxwxi/NED6dOtU4eXy7dOYDPrbGZdI18DJwLz/SqPpDetICl+qqio4MMPP6SwsBCAGW+9jYvJZ+T02JdLgvvywc3H86tTDkrIwx/8HQXUC5hpZuXAx8ArzrnXfCyPpDHtsyB+ibQ+a2pqKC8v56Ibbmfq112x7JzoMbXrlvLgDZeyYE54NHyiJoP5FgCcc8uccwUN/x3inLvdr7JI+tMKktLR4u1SVlxczIUXXhhtfTrnmHLfLezTLZ8Bh4xodH5NTU20ZZqoyWB+dwKLJIxWkJSO0rSPqV+/foRCIazTHmzK6wn2FbhwEAjgOH6PjRw/7nR+Vb7tAZ+VlUVxcXFC+6sUAERE2ii2j6mmpoYlS5Y05PgdYBCw8JdAXl5eNL0TGfNvZjz00EMEg0HuvPPO7fqrFABERHwQb/5IYWEhFRUVTJkypckELofDcM4RfeLjyLIAPfv0IT8/P3oOwLBhwxpdJ3K/RO14pgAgIrIDTeePNE3RvPLa66zI7k+fsf+PLWuX07NHdxa/+ACuPhS9Rm5uLtOnT9+uJh8vLRkMBpk5cyYlJSUUFxd3aH+VAoCISDPi5eNj0z3VNbVccNtT5BWdzcjjTmDCsftxwsG9KBrxFgsWLMDzPHr27Bn34b8jwWAwIQMVFABERJoR+7DfunUrJSUlDB1+BJaVg/PqIZBFUfAYfvvTIIcP3Ct63iOPPMIxxxyD53l89913fhV/pxQARETiKCwsZN26dZhFcvrwwKtzyN88jF7j/hPvs7d46LYbOefU0dudG2/iYTIOPVYAEBFpIjb1k52dQyhUB8DamdM4+weH8Mhjv6Rn1182e34iO3LbIhn2AxAR8U3sBK7IDNz/nfEWnheu9Uce/mGOfz50G0vn73hOSapMPFQLQEQyVtNO3n0GDGRLXg+Wz9qIZefg6kMYLpoCAvA8r0UpnVSYeKgWgIhkrNhlGrZW17B08Res/XQmq/7nUXJ3705WwDjhhBManROZsZsO1AIQkbS0sw2AiouLWVdZCxYA58AMPA9wuPo6ajetwTnH22+/3ei8K664ImlTOrtKLQARSUvNLajmnOPdLzbwmdeH779dG16jx4yrJ11HIGBAuJYfSfvU19djFn69U6dOXHLJJYn9IB0oc1oAKz6Emkq/SyEiCVA6ZyHz5pXjeY7Rx49ixtN3U1RwEC8v85j8qcfcTxey7v1/gNew0brzmPrUZPbr34dQfT03XTmOq255AM9z5OXm0K9Xd0L19Uz5w00Eu1fB4jcT/6HyusKAI9v1kpkTAGoqoXN3v0shIu2keOwVAFRUVlGxuZIpD95OsKgAgJI5S6KjeGrrQtzzykLWfnEAq6vqGbpXNqNylzA18vBvsGlzFdW1dcx47jGCRQUcNmw4JaWzKQ6OiF7XV1s2tvsllQISkZRUUVnFF8u+onzB5yxfsZpjfnw5pWXlABQHR4Q3ADLDsyxKvSH075bFk6fvxesX9mDSmUdG0z2RjYIgHCxKSmcDECwq4OZJlyfHw7+DZE4LQETSRmlZOfMWfhGt5QN4nuOZaS8z4T8fY73XlT1GX4lXU8XxRx3Ob88dSUGv3OixwaIChh08lIrNldx0zXiuuvkOPM+Rm5NNcXBEvFumJQUAEUlahWPO3y69A1BSOrvRwz9i8t/+hQfgeWTnvMfzf3mYHx/3g7jXnvP61OjXj/51Wtz77Gq5Uo3vAcDMsoAyYJVz7jS/yyMiySG2lj963MRobr547BVUVFbFHGlE1t736uujX7v6EJ/NmwvNBIBYscGgteVKRcnQB3AdsMjvQohIcigeewXFY69oVMuPzc1XVFaxdsMmwg9+IBDAsrIBIy83h4YRmx2WzmmuXKnI1wBgZv2AU4En/CyHiCSPisoqVqxaw957dot21EYe5q+9P4fyRUtZu249kZp+AEev7nswaEAf3v7H4xQccgCDBuzTYTXzcAdz43KlKr9TQPcDNwJdmzvAzCYAEwAGDBiQoGKJiB+G/vAsFi9bAcD1t9zLfgP7EwqFuPfO23ilcl8ef+IvuIalGyLycnN4YfI90Yf9rqRzWiO2A1l9AK1kZqcB651zs82suLnjnHOTgckARUVF2/f6iEhaKC0rZ8nyFdHva2rrGHPyqYQKzubn86rJCXzPmB8eztR3nyXS/xsIBLj/1hsS/hDu6CCTKH6mgI4GzjCzL4G/A8eb2bM+lkdEEqhwzPkMDp4WHbtfUjqbmEU3cRjTKgbz/soarirqwszxvXh24tEMO3gIfXp1x8zwPI/rb7k3eg3ZNb4FAOfczc65fs65gcB5wFvOuYv8Ko+IJE5kJM3yFasZPW4ipWXlHDVyBAQCgEEgi0GnX8PvzzuCD8b34sajdqdn5ywgXPueNP686Fo9qd4R6ye/+wBEJAPFjqTZWl3DiRNuof/p19H18B/z/cJ3Of2UE3n+nkvJybK450c6YjNx8lZ7SoZhoDjnSjQHQCS9NE3xxIou1dCgat3XfP7UL6n8aDr1lRt48bkpPDX1hWavHemI7cjRPpkgKQKAiKSX2BRP7Bo9AAs31PH8NwPoed6d5A8c3vCqw6sPNbrG9Fdn7PAec16fyrLSl/XwbwOlgESk3cWmeDzP8fS0l6nveSCPf1LFeytq6JxjXH36SApO35szz/8/AGRnZxEKbVuh8+xTRvtS9kyiACAibRJZlrlk2p+2vRYcgZlFO2on/+0l/uWOpN+QQ7jxqK5ceGhnuuUHKC3Lip4TCATIzTHqPccjd9zEhIvOTuwHyUBKAYlIm0Rm7sameYYNO4y8XoOi3zuvnuKcxcy8rBdXF3WlW3740RM7eqe+3uN3P59IaMUsPfwTRAFARFqtcMz5zJ3/eXQ45ysz5/BfH2zmsAeXs1vByeH9doFOeTn87KwgedmNR/Wk07IKqUgpIBFpkXipnnUbvol+XV1bx6WPvEPXI3tR/dFUvn1/2yie5mbrptOyCqlIAUBEWiSy9WJpWTnBogI+KCtn7YZvo+87C3Bq8eGcMXQtY+9uPIRzzvzPmr1uuiyrkIqUAhKRnWo6c/e+F2dx5h0vNlq64bJzz+DpCUfzxadz/Suo7BIFABHZTmRN/ohGM3dr6rjtufeh815ElmTulJ/HhPNOD5/bMAIoIjs7i0vGap5nMlIAEJHtxI7s2bTVY123g7DsXLAAgewcLjy8N9+W/AWcw6xxjj9YVEDBIUPJycmmT6/uvDv9CeX2k5T6AEQyQNMO3Lhj9xteu/Oma6NbHh57zkT2Of92rPeBDDh5IvWfzeC5e2/gnZgWQSCQxTebKhrdT3n91KAAIJIBysoXEqqvj3bgNu3QhW2dvPc981J0vf1QXR0Dtizizxf+kKF7XwZcBoQ3Y9RibKlPAUAkzZWWlbPl+60AHHXmeG68+tLtNjW/+td3M3fBYsCxfNW/IZAN9XUYjosP783QvXMaXTNYVMDMfz5JSelsioMjlOJJUeoDEElzz0x7udH39zz6TKNNzR/814csdb3AOXCOgPM47fijMQPnXLMbrgSLCrh50uV6+KcwBQCRNLc2ZrIWEF2fBwxnWby+Kpuqz94nMqInOzuLfr33jg7x1IYr6UsBQCTN9e6x93avddrvcPY49mJG/+y/OHtQPS5mKebLx53BJWNP0xINGUABQCRNNB27H1F46IGNX7AAp5x7Me/98Wre+PmxjD/5iOjDvlN+HpeMPU0brmQIdQKLpImmI3ucc5SurOWx99cQHrcTzulcdO5Z/PXaY6LnNbcej4Zypj8FAJEUFTuWv7SsnLnzPwdg9LiJ3HrvA/xpoVHbdR/yex6CZWXh6kN0ys/j6gtO3+5aethnJt9SQGaWb2Yfm1m5mS0ws1v9KotIKtnjoGPZ46BjG83Wje2k3VpTyx3Pv8/367+i7t+38tgpe1Fw0H5K58h2/GwB1ADHO+eqzCwHmGlm/+Oc+9DHMokkvZqaWmrrQtEa//HjJnLkmJ9sO8A5TuxXz7S/34HnOU65YKIe/BKXby0AF1bV8G1Ow39uB6eIZLzJz06nuqYWz/Oir1XX1PFh+aLo94FAgO9WLW001l/DOCUeX0cBmVmWmc0F1gNvOuc+inPMBDMrM7OyDRs2JL6QIklk+qsztnvNzNjL2xT9Pi83h7NPGa1hnB2kcMz5DA6eFndyXKrxNQA45+qdc8OBfsBIMzs0zjGTnXNFzrmiHj16JL6QIknAc47Xl25lfa8jtnvPefWs+XoFQHRlzgkXna1hnB2g6b4IqR4EkmIUkHPuOzMrAU4C5vtcHJE2i7faZmtUhxwvfPY9T3xSxbLv6ulXcBL9F89lw6IPqa6uiR4Xmd0buzKnRva0v9h9ESKptVQOrr4FADPrAdQ1PPw7AT8C7varPCLtKTImv3DM+XTr2mWXA8F31R5/nbeFp8u3sHGrx2E9c3jwpN05ef98si+7j9Kyco4+a3yjHblA6Z6OFtnEPl1WQfWzBdAHeNrMsginop53zr28k3NEEq5wzPm7tGl5JE3geeHNUnr37N5o2eVYTVsKR1/0Czb3DRLaN8jWkKN43zwm/KALwX65jXbZChYVsP+gASxeti3107tnd6ZPviela6TJLt02sfctADjn5gGFft1fpCViH+aRpZN39o8+Nk3gHKxZt7HZcyMthUPPmoR36Gls2mcM1V99yph9e/K7C47gwO458W4BwBfvvbjLwUnaLp1Sa0nRByCSrFqT841NE0REzr35roeir1VUVlG+cAnOebBqPZ021bF12Wzw6nlp1vNcN/Ix6L7je6XTw0gST4vBiexA5GEOLc+vRzZLmXjx2dG0TeTcisoqvlq1juXf5/N5RTbO88LNhPoQW5d8DF49oLH7khhqAYjsQGtyvrErcu7WKZ+e3fdkyoO3s6XOUb5oKa6+Hlas3OE10qGDUZKfAoDITuxqmiWS18/OzqZn9z154L/+k5Lqgfzxsb+Ea/wxE94DgQBFww5iVvlCnHPk5eYyftzp0SWZRTqSAoBIO4rtNA4zzrjoKvYYdQVu9WcYrtF6J3m5Odx/6w1c/eu71JkrCacAINJGsUM5L5j060adv+BwoVo2vfkwEF62ISsrC4Ce3feMDttUZ674QZ3AIq0Qu/tWuGN3Lff+cxZbBo2CQFaz5xngefXU19fzXUVlgkorEp9aACKtEMnzP/T0tHDHrufxi+uuBRx49eGafiBAvefFbMIOmOHSaCkBSW1qAYjsonCefzHLV6xh0m/uDo/qcR7mhaBhc3Uzo96rb/TwN4MbJl6sVTolaagFILILln8X4uZn3ydciXfRcfsAuTk51NbV4hxkBQLUhbat2b9nt6688swDBIsKOGtMMSWlsykOjlDtX3ylFoBIC3yyppYrX/mW459Zz9JOBzRalwfgrJOKefsfj1NwyAEMGrAPD93+y2hNv1N+XvThD+G5BTdPulwPf/GdWgBpqL2WIk538X5Osa95zjFjeQ2Pz66ibE0t3fKMaw7vwqUFxRS8vidr1m2Mnvfy/77HjVdd2mg0z2EH7q+aviQ1BYA0FOmgbG4FSgmL93OqqKziu81buO35j3inZhBLN4Xo2zWL3x67O+MO3o3OueFG8/TJ9zRajtnz3HYdusGiAv38JakpBZRm0m3Hoo4S7+f05gdzKV+0lC+/Xs3vfnEdW1cu4o9j9uCdS3ty+fAu0Yc/bFuOOSIvN0cdupJyFADSTLzVKyUsdux+05/T//vb+1z4yDsNi7N5mAtxWpclnHnAbmQHLO71vnjvRYYfeoC2XZSUpRRQmkm3HYvaU2zKJ/xzCuA5h2dZfJZ7AMcflcc/3vkrnoO8nBxGBYt2ek3N4JVUpgCQZtJtx6L2ErtGz6hzJ3LcNXfT84K7qf5yDueeeBS3nHsE+3TN4vN/DdHPTjKGAoBPOnInJ9VKt//5zni/DK9hWH5NbR0L5s3ht5Mu54LDRtEtb1smVD87ySTqA/CBOmo7VtOf781/+5DnNu2HZeeABcjNzWHK1cVcVdS10cNfJNP41gIws/7AM0BvwAMmO+f+uMOTNi6Gp05t3Q23fgtZze+vmkglr31NZPxgbW0NJU/dSnBNf59LlXrmLvgCgOGHDG30eqOfb00N3716B38/vh8fj92LT76q4oIf7kdw/ZPw74QXWaT16utg5v3tekk/U0Ah4OfOuU/MrCsw28zedM4t9LFMCVE8tBs5WRCqh9ysAMVDu/ldpJRUvqaWslV1XNBpM8HBuwOwtc6xf+8u5GQFqKv3yMkOMHb4XhzSI4dDiocw3ucyiyQTa7RSYbwDzK4FpjjnNnVoQcz+BTzknHuzuWOKiopcWVlZ626w+E3o3L2VpWt/HdkHkAlKy8o55seX43mOTvl5PPDQg5TWDWLGlzXkZxt1s6cR+uJd/n7fb/TzlfSwZSMMOaFVp5rZbOfcdsPaWtIC6A3MMrNPgCeB193OosauF24gUAh8FOe9CcAEgAEDBjR9O2Wps7Hl4gXL2HH8W2vq+MUzMxk0qh/XH9GVS4btxl5XTwIm+VhqkeS30x4w59xvgCHAn4HLgMVmdoeZ7dceBTCzLsB04Hrn3OY495/snCtyzhX16NGjPW4pSSB2UtaO3o/XYb61zmNrz0Ow7FywAIHsHCadFeT98T25/oiu7NWp+Q1ZRGSbFvUBOOecma0F1hLO3e8JTGvI2d/Y2pubWQ7hh/8U59wLrb2OpJ6drVcUef+ZaS83mrH7++c/YM0nvfi2uj99T78et+gN/n7XdRxzuNI8IrtqpwHAzH4GXApsBJ4AfuGcqzOzALAYaFUAsPB6un8GFjnn/rs115DUFDspa/S4idstoxD7/uPPNtQLzPAsi09sKKf0zmHCiC6M3GccZuf59Ckyh1aXTV8taQF0B37inPsq9kXnnGdmp7Xh3kcDFwOfmtnchtd+5Zx7tQ3XlCRXPPYKVqxas916RbEBIDa/75wjt/cQcnsN5swzTuN3545k/72SYzhvptDqsumrJX0Av2368I95b1Frb+ycm+mcM+fcMOfc8Ib/9PBPcxWVVVTX1EY3VMnNyeb5l9+M1jI95+gy8DAsK4fwFupQu3Yxoc9LuKaoS9I+/HfWp5GqNGkxvWkapCRMaVk5c+d/zpp1G8nJyWbixWcz47nHAPhq1Tp+/4+POfHZDfxheV/2OOZCAvmdo+fWJfnKphWVVaxYtSbtHpBaXTa9KQBIwsQ+PEKhegb07UNVnYuuwf/bG35G1dcLmdh3JRXvPY1XXRU9PplXNk3nWnJkdVlI7t+BtI4CgCTM3ntum/HseR4fbszlooe3rcHvQrWs+/e9uNXzo7XOQCDAiccdmdTr7adzLTmyuqz2PEhPWg1UOlTsCJI58z9r9N7bZYs4+ZSTmdawBr8Z1NTUsvee3aJ7GuTl5vC7/7gyqR886b4HgyYtpi+1AKRDVVRW8dWqNTz2Shkzltc0eu/MAzrx3DXHMOzgIfTp1R0w1qzbyPW33Mt+A/unTK1TtWRJVWoBSJs0HfkSGStePPYKKiq/p3zhEpzzuPrqa9n72EvAAuA8cnNzuPqC04FwDfPOB5/kV3c9BITTKOPPPYObJ12e2A/TBqolSypSAJA2iYwRz87OJhQKUVpWzrBhh7EmbyArVn0Rzu/jcKFaQrOfZ8igfoRCoe0WwUtkGkUTm0TCFABkl8Q+PGNn7EYcd85Edj/2Ur6d+W9cqLbRud9VbKampiZumiSRW1lqYpNImAKA7JLYh2fs6JeIuro6tsx9tdHD38yILCAbb+ZvxK6mUVpTk9/ZMhQimUSdwB0k1WeGxit/0/HuldY5PGO3YVYvQHaWUbNxZaPzYlcPz87Karf0TmsmX6XzkE2RXaUA0EFSfWZovPI3WoO/upb7pr7BXqMuI3jC6dFjnHPsaLOI8eNOb5cad2snX2lik8g2CgAdINVnhsYrf3XIURNZgx8DHNVfzuHbGU/QK7BtG4cdbRWUl5vLJWPbsn7gNq2tyWvIpsg26gPoAPEeTsn+oInNpzct/x3/+IDVc3rxzdb+9D3jeja/91c2b1gDhGv8/37zPczCD/+83Bzuv/UG5sz/jE8+/YyP5y6I3qO9av/QtlFDGrIpEqYA0AFScWZobOduuPwBPOfwLIsyhnJSrxwm/KALR/Ydx4ezD+Los8ZHa/vOOXr37E5+Xm6jETxN9+1tr9o/JHbUkEi62umm8MkklTaFT6VN3wvHnM/c+Z8DkJeXR+8fXU51t0HUfD2P2jWfU7tyIQ/fej0TLjo7es7QH57F4mUrAOiUn9dsOiWVfg6SXtLub8+nTeGlFVIhzRD5B7K5ckv0tZraOr75eimdVy7kwsMH8PB7swC48pe3A0SDwBfvvRgdClocHNHsP7BU+DlI+tFw35ZRAEgDLRkP3/SY8D+QxXie1/hA51E17w2qgEc/bTxGYPqrMxq1AoJFBQSLCtKvpiUpLxX74fygUUBpoCVDTmOP2Vzjcfe0D/B2kv1rmh48+5TR2x2T6iOeJD1puG/LKACkuJY8gGOPOe6ciQy/rYSP3NDwwmw7kJ+XS25ONllZWTx+968b1f4jNLFKkpGG+7aMrykgM3sSOA1Y75w71M+ypKqWNHWfn/ExkUxPXV0dWz59k6qvPgUvFD1m7z278e13FTgXXrqhd8+9mT75np3+w0nFEU+SGdT/tHN+9wH8BXgIeMbnciS1eDn+yGt33nRt9LWAWfQBfNzYK6jeYxAHnn4lb3w7GMvOwdXXkZ+XQ+6auXy/cd228wIBfn7lxRQHR+y0U7cpDccUSV2+BgDn3LtmNtDPMqSCeKtXRl578fWS6HF1oRAvvFbCht2HsqzPj6jauIaK2fP49diR/OnTA6nctJGbrhnPxJtub3T9QMCiD/3WPMBV0xJJTX63AHbKzCYAEwAGDBjgc2kSL95wtqt/fVd03P69j/210fH3Pvo0Ty7rwqYZf8KF6qj5+DlGnPEY1776FwDufPDJ7ZZrePj2m9K25q61/0Wal/SdwM65yc65IudcUY8ePfwuTsLFy/Gv2/BN9P14E/m6r3irYTlmt13HbOzoCDOa7dxNF6m+KJ9IR0r6AJCuWrpcdGynam5ONnvv2Y216zduOyCQTXaPQY3OOXjfXtGvs7ICja4ROzri/RefSuuHv4aoiuyYAoBPWlozDRYVMGTwAPbs1pX7b72BeSu+xRFZf98oOO4kpt736+jxebm5nDzqaKxhjX7DtrvmnNensqz05bRN+0RoiGpqKxxzPoODpylwdyBfA4CZTQVKgQPMbKWZ/dTP8iTKjmqmTVsGpWXlLP3yazZVVHL1b+5l6pJAdBOW/PxcHrnuJ4wd9QOGH3oAgwbsw9v/eJxvNlVEU0Oh+vqMffBpMlDqUustMfweBXS+n/f3y47G7seO+CkcPoz7XyyNztitD9Uxqncd1zz3GB/O+qTRcM2mI3HSaWx+aztyNUQ1dWkph8RI+lFA6ai5yVOxI36OO2cigy++g4q6IVhWNi5Uh+Fx4sE9GHXkcEYdObzZ66fbg68tm7hriGpq0gTDxFAfQILEpnaCRQV0ys8nLy+30TT1F0tmRWs9dbW1dN6wgOnXHsNjt98Y3XDl+lvubVFzOF3y/EoFZCYt5ZAYagG0UUvTExWVVSxavJwuQ46mqOBghgweQMXmSgDmr6/l8U+28LcFsWc4rjy6D0f1z+OdTRXRsfuZ1hxOxVSA5h60D7XeOp5aAG3UktE8kVpsTU0tW77fypz5n1G+4HOWr1jND8dO5IQ/vEfJl9W4ZR9EzwkEAnyzqQLI7M7MVPzsmnsgqUIBoA2aS080HckTW4sF2Fy5JVqjrw/VcaQt5r6CdXy3dNsDIzs7K/qwy+TmcKp9dqWs/NHSeTXSmFJAbRAvPXHzXQ9RVr6Qnt33jHZaFgdHYGZxZ+12yuGynYIAAAtvSURBVMvhxrFBSkpnN3r/8nFnNHrYZXJzOJU+eyqmrNJBWwYKZDK1ANogXnpi9boNbPl+a7QG+Mp7cyipGcSeP5qw3fr7fXp1j9ZqY6/V3huoS+KkYsoq1anV1XoKAG0QLCogNyeHrKws7r/1BgCWLF8Rfb+6po5LHn2HP31SxfnjzmHQgYdE38vLzW203n6qpTokPv0eE08zvltPKaBWiF2Lv7auDs9zXHXzHVxx4U8arbTpzDhr1EhuObcnA7plc+eS4/jVok+BbTN0leZJP5n4e/RzX2jNGWg9tQBitHTtkcgoj2emvRyteXie419z1mDZuYCBBdhv8ECWvv4UA7qF46zSA5KO/E7BqNXVegoADVr6R1xaVs7c+eEhnE8991KjvP66RbPod+gRDNy3Lx+8+Ge65mc3Gg4YLCpg5j+f5I6brtUfqqSNZEjBpMvEx0TL2ADQtLYf74843tCy2D/umroQub33j34fwOOqkwpY/sFLAHEDSrCogJsnXa4/VEkbatmmrowMAPFq+/H+iJtO6Fn+XYj31uduu5BzjPvRyLh//MlQKxJJBKVgUldGdgLHPpy3VtfwzLSXefSuXzPs4KEsWrycA/Yf2GjbxVHnTuSk6++hnMFULt22GUsgEOCgvnvEXXhNHVOSSTKx4zsdZGQAaDox68nnXuKSsacx5/WpdBp8JJ8uWkLn3TpFj6+prePd9z+kf/+VdK2Yx6aGB3tebg7FwRHcPOny7e6Rbityikj6ycgAECwqoHfPvVmzLlybr6/3KCmdzaefLaG6phaAzZVV205wjrzFM9i6Oo9QKMR+A/sTCoV2+mBXrUhEkllG9gEATJ98T6PcfdGIH3Df1DeaOdqxdt06Fi9bwfIVq1m5ep1q9SKS8jI2AERSNPsOGsyFN/+B68p7sb7XyBadq05dEWmLZFm8zu89gU8ys8/NbImZ3dTR94sd+rlgQx3H/scjZI17gLerBzJmcD7v3X1JdAP2C39yMrb9fuqAhrqJSNsky5LhvvUBmFkW8DBwArASmGVmLznnFnbE/UrnLIyO6vnh2In0GPd79hp4MOOHd+a1x25jTslmDp72J75478XoOQu+WEbF5kpuumY832yq4KnnX2pR7l9EpDmxW7+OHjfR16GzfnYCjwSWOOeWAZjZ34EzgXYPAHX1Ho+8OhcwwFEfqmWkfcEzl4+iW16A6b9fHXcp2aaduPFG+4iI7IpkWjLczxRQX+DrmO9XNrzW7n45fR6vru0MNKzU5hynHtqTbnkB39cxySTJkvcU8VMyzZz2MwDEy7Bvt2OKmU0wszIzK9uwYUOrbnTZUQP5cd8t0e9jt1vUjN3EiZf3VFCQTJNMM6f9DAArgf4x3/cDVjc9yDk32TlX5Jwr6tGjR6tuNKzfHvx0zPBo1I1M4ILkisbprLmWVrJ0hokkUrIsXudnH8AsYIiZDQJWAecBF3TUzYKFB8edmRtZobOkdDbFwRG+/0LSVXMtrWTpDBPJRL4FAOdcyMyuBV4HsoAnnXMLOvKezc3MDRYV6MHTweKtjZRMnWEimcjXpSCcc68CrybkZnldYcvGnR8nHSJ4UF+GHTCYisotTPnDTQQP6gvVFY2DQuH++h2JNCeva7tf0pzbrt81aRUVFbmysjK/iyHtqLCwkIqKCqZMmUIwGPS7OCJpycxmO+eKmr6ekYvBSfKYM2eO30UQyVgZuxaQiEimUwAQEclQCgAiIhlKAUBEJEMpAIiIZCgFABGRDKUAICKSoRQAREQylAKAiEiGUgCQqMLCQgYPHkxpaanfRYmruLiY4uJiv4shkjYUAASA0tJS5s2bx/Llyxk9enRSBoGKigpWrFiRlGUTSUUKAAJASUkJnucBUFtbS0lJib8FaiIVApRIqlEAECCcXgkEwn8Oubm5SZdqSfYAJZKKFAAEgGAwyLBhwxg0aBAzZsxIuqWZkz1AiaQiLQctUcm8NHMkQGnvAJH2owAgKSOZA5RIKlIKSEQkQ/kSAMzsHDNbYGaemW23TZmISCpJ1TkqfrUA5gM/Ad716f4iIu0mVeeo+BIAnHOLnHOf+3FvEZH2lMpzVJK+D8DMJphZmZmVbdiwwe/iiIg0kspzVDosAJjZ/5rZ/Dj/nbkr13HOTXbOFTnninr06NFRxRURaZVUnqPSYcNAnXM/6qhri4gki1Seo6J5ACIibZSqc1T8Ggb6YzNbCQSBV8zsdT/KISKSyXxpATjn/gn80497i4hIWNKPAhIRSRWpNiFMAUBEpJ2k2oQwBQARkXaQihPCFABERNpBKk4IUwAQX6VazlSkOak4IUwBQHyVajlTkeYk+6568WgimPgmkjP1PI/Ro0enzD8akeak2oQwtQDEN6mYMxVJJwoA4ptUzJmKpBMFAPFNKuZMRdKJ+gDEV6mWMxVJJ2oBiIhkKAUAEZEMpQAgIpKhFABERDKUAoCISIZSABARyVAKACIiGUoBQEQkQykAiIhkKHPO+V2GFjOzDcBXrTy9O7CxHYuTCvSZM4M+c2Zoy2fe1znXo+mLKRUA2sLMypxzRX6XI5H0mTODPnNm6IjPrBSQiEiGUgAQEclQmRQAJvtdAB/oM2cGfebM0O6fOWP6AEREpLFMagGIiEgMBQARkQyVEQHAzE4ys8/NbImZ3eR3eTqamfU3s7fNbJGZLTCz6/wuUyKYWZaZzTGzl/0uSyKY2R5mNs3MPmv4Xaf9nppm9n8b/qbnm9lUM8v3u0ztzcyeNLP1ZjY/5rW9zOxNM1vc8P892+NeaR8AzCwLeBg4GTgYON/MDva3VB0uBPzcOXcQcCRwTQZ8ZoDrgEV+FyKB/gi85pw7ECggzT+7mfUFfgYUOecOBbKA8/wtVYf4C3BSk9duAmY454YAMxq+b7O0DwDASGCJc26Zc64W+Dtwps9l6lDOuTXOuU8avq4k/GDo62+pOpaZ9QNOBZ7wuyyJYGa7A8cCfwZwztU6577zt1QJkQ10MrNsYDdgtc/laXfOuXeBb5u8fCbwdMPXTwNntce9MiEA9AW+jvl+JWn+MIxlZgOBQuAjf0vS4e4HbgQ8vwuSIIOBDcBTDWmvJ8yss9+F6kjOuVXAvcAKYA1Q4Zx7w99SJUwv59waCFfwgJ7tcdFMCAAW57WMGPtqZl2A6cD1zrnNfpeno5jZacB659xsv8uSQNnAD4BHnXOFwBbaKS2QrBry3mcCg4B9gM5mdpG/pUptmRAAVgL9Y77vRxo2G5sysxzCD/8pzrkX/C5PBzsaOMPMviSc4jvezJ71t0gdbiWw0jkXadlNIxwQ0tmPgOXOuQ3OuTrgBeAon8uUKOvMrA9Aw//Xt8dFMyEAzAKGmNkgM8sl3Gn0ks9l6lBmZoRzw4ucc//td3k6mnPuZudcP+fcQMK/37ecc2ldM3TOrQW+NrMDGl4aDSz0sUiJsAI40sx2a/gbH02ad3zHeAm4tOHrS4F/tcdFs9vjIsnMORcys2uB1wmPGnjSObfA52J1tKOBi4FPzWxuw2u/cs696mOZpP1NAqY0VGyWAeN9Lk+Hcs59ZGbTgE8Ij3SbQxouCWFmU4FioLuZrQRuAe4CnjeznxIOhOe0y720FISISGbKhBSQiIjEoQAgIpKhFABERDKUAoCISIZSABARyVAKACIiGUoBQEQkQykAiLSBmR1uZvPMLN/MOjesVX+o3+USaQlNBBNpIzP7PZAPdCK8Ps+dPhdJpEUUAETaqGEphllANXCUc67e5yKJtIhSQCJttxfQBehKuCUgkhLUAhBpIzN7ifAy1IOAPs65a30ukkiLpP1qoCIdycwuAULOub817D/9gZkd75x7y++yieyMWgAiIhlKfQAiIhlKAUBEJEMpAIiIZCgFABGRDKUAICKSoRQAREQylAKAiEiG+v/iWihdaDsr2wAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "xs = np.linspace(0, 10, 100)\n",
    "plt.plot(xs, result.x[4] + result.x[3]*xs)\n",
    "plt.axhline(result.x[0], color=\"C1\")\n",
    "plt.fill_between(xs, result.x[0] - np.sqrt(result.x[1]), result.x[0] + np.sqrt(result.x[1]), alpha=.2, color=\"C1\")\n",
    "plt.errorbar(x, y, yerr=yerr, fmt=\"k.\");\n",
    "plt.xlabel(\"x\")\n",
    "plt.ylabel(\"y\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}