MLE and Gradient Ascent Algorithm.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We study Table 13.2 on p 240 of Lange and Chapter 13 question 11 (on page 240) and Chapter 14, Example 14.10.1 on p246.\n",
    "\n",
    "The table contains counts of the number of days a given number of deaths occurred. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "freq = np.array([162,267,271,185,111,61,27,8,3,1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYEAAAEBCAYAAACe6Rn8AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAHmhJREFUeJzt3X9Q0/f9B/AnCQRNglJae7SnCYijRZyj4NZeFVzn2E3PXrfWQSCl3dma2Wp/KEvpTEs5JAqz7LyuVURbtDRF7LSt1Vt/MB0odv2RK2exwV5VQKejKrQ1ScMP8/n+4YjmK9DA/OQD+Twfd94l73c+n88rbzmevD8/wwRBEEBERLKkkLoAIiKSDkOAiEjGGAJERDLGECAikjGGABGRjDEEiIhkLFzqAobDbrdLXQIR0ZiUlpY2YPuYCgFg8C/yQxwOB5KSkq5xNWMXx8Mfx+MyjoW/UBiPof6A5u4gIiIZYwgQEcmYKLuDvF4vioqKcPToUahUKpSUlECv1/v6t27dir179wIA5s6di+XLl0MQBGRkZCAuLg4AkJKSgvz8fDHKIyKi/xIlBOrq6tDT04Pa2lo0NTWhtLQUGzduBACcPHkSu3fvxhtvvIGwsDDk5ubil7/8JcaPH4/k5GRUVFSIURIREQ1AlN1Bdrsd6enpAC79Rd/c3Ozri42NxZYtW6BUKqFQKNDX14fIyEgcOXIEHR0dyMvLw5IlS3D8+HExSiMioiuIEgJOpxNardb3XqlUoq+vDwAQERGBmJgYCIKAsrIyTJ8+HfHx8Zg0aRJMJhOqq6vxhz/8AWazWYzSZM9msyEuLg7JycmIi4uDzWaTuiQikpAou4O0Wi1cLpfvvdfrRXj45U11d3dj1apV0Gg0eO655wAAM2bMgFKpBADMmjULHR0dEAQBYWFhfut2OBwjqsnj8Yx42VCxZ88eFBYWwuPxAADa2trw8MMP4/Tp01i4cKHE1UmLPx+XcSz8hfp4iBICqamp2L9/PxYsWICmpiYkJib6+gRBwKOPPorbb78dJpPJ1/7iiy8iOjoaS5YsQUtLC26++earAgDAiM/XDYVzff9X8+fP9wVAP4/Hg5deekn2My/+fFzGsfAXCuMx1HUCooRAZmYmGhsbYTAYIAgC1qxZg6qqKuh0Oni9Xnz88cfo6enBgQMHAAArV66EyWSC2WxGfX09lEol1q5dK0ZpkrHZbLBYLGhvb4dOp4PVaoXRaAxqDe3t7cNqJ6LQJ0oIKBQKFBcX+7UlJCT4Xn/++ecDLldZWSlGOZKz2WwwmUxwu90ALu2G6Z8FBTMIdDod2traBmwnInnixWJBYLFYfAHQz+12w2KxBLUOq9UKtVrt16ZWq2G1WoNaBxGNHgyBIBgtu2GMRiMqKyuh1+sRFhYGvV6PysrKoO+WIqLRgyEQBIPtbpFiN4zRaERrayuOHDmC1tZWBgCRzDEEgoC7YYhotGIIBAF3wxDRaDXmnicwVhmNRv7SJ6JRhzMBIiIZYwgQEckYQ4CISMYYAkREMsYQICKSMYYAEZGMMQSIiGSMIUBEJGMMASIiGWMIEBHJGEOAiEjGGAJERDLGECAikjGGABGRjDEEiIhkjCFARCRjDAEiIhljCBARyRhDgIhIxhgCREQyxhAgIpIxhgARkYwxBIiIZIwhQEQkYwwBIiIZYwgQEckYQ4CISMbCxVip1+tFUVERjh49CpVKhZKSEuj1el//1q1bsXfvXgDA3LlzsXz5cng8HpjNZpw/fx4ajQZlZWWIiYkRozwiIvovUWYCdXV16OnpQW1tLfLz81FaWurrO3nyJHbv3o3t27ejtrYWBw8eREtLC2pqapCYmIjXX38dv/nNb7BhwwYxSiMioiuIEgJ2ux3p6ekAgJSUFDQ3N/v6YmNjsWXLFiiVSigUCvT19SEyMtJvmYyMDHz44YdilEZERFcQZXeQ0+mEVqv1vVcqlejr60N4eDgiIiIQExMDQRDw5z//GdOnT0d8fDycTieioqIAABqNBhcuXBhw3Q6HY0Q1eTyeES8bijge/jgel3Es/IX6eIgSAlqtFi6Xy/fe6/UiPPzyprq7u7Fq1SpoNBo899xzVy3jcrkwYcKEAdedlJQ0opocDseIlw1FHA9/HI/LOBb+QmE87Hb7oH2i7A5KTU1FQ0MDAKCpqQmJiYm+PkEQ8Oijj+KWW25BcXExlEqlb5n6+noAQENDA9LS0sQojYiIriDKTCAzMxONjY0wGAwQBAFr1qxBVVUVdDodvF4vPv74Y/T09ODAgQMAgJUrVyInJwcFBQXIyclBREQEysvLxSiNiIiuIEoIKBQKFBcX+7UlJCT4Xn/++ecDLvfCCy+IUQ4REQ2CF4sREckYQ4CISMYYAkREMsYQICKSMYYAEZGMMQSIiGSMIUBEJGMMASIiGWMIEBHJGEOAiEjGGAJERDLGECAikjGGABGRjDEEiIhkjCFARCRjDAEiIhljCBARyRhDgCRhs9kQFxcHhUKBuLg42Gw2qUsikiVRHi9JNBSbzQaTyQS32w0AaGtrg8lkAgAYjUYpSyOSHc4EKOgsFosvAPq53W5YLBaJKiKSL4YABV17e/uw2olIPAwBCjqdTjesdiISD0OAgs5qtUKtVvu1qdVqWK1WiSoiki+GAAWd0WhEZWUl9Ho9wsLCoNfrUVlZyYPCRBLg2UEkCaPRyF/6RKMAZwJERDLGECAikjGGABGRjDEEiIhkjCFARCRjDAEiIhljCBARyZgo1wl4vV4UFRXh6NGjUKlUKCkpgV6v9/tMZ2cnDAYD3nnnHURGRkIQBGRkZCAuLg4AkJKSgvz8fDHKIyKi/xIlBOrq6tDT04Pa2lo0NTWhtLQUGzdu9PUfOHAA5eXlOHfunK+tvb0dycnJqKioEKMkIiIagCi7g+x2O9LT0wFc+ou+ubnZf6MKBaqqqhAdHe1rO3LkCDo6OpCXl4clS5bg+PHj16SW/oeXJCcn8+ElRET/T0AzgcbGRlRVVaGnp8fX9uqrrw76eafTCa1W63uvVCrR19eH8PBLm5s9e/ZVy0yaNAkmkwnz58/Hp59+CrPZjJ07d171OYfDEUjJAIA9e/agsLAQHo8HwKWHlzz88MM4ffo0Fi5cGPB6QpHH4xnWWIY6jsdlHAt/oT4eAYXA2rVrsWrVKsTGxga0Uq1WC5fL5Xvv9Xp9ATCYGTNmQKlUAgBmzZqFjo4OCIKAsLAwv88lJSUFVAMAzJ8/3xcA/TweD1566SWYzeaA1xOKHA7HsMYy1HE8LuNY+AuF8bDb7YP2BRQCN910E+68886AN5iamor9+/djwYIFaGpqQmJi4g8u8+KLLyI6OhpLlixBS0sLbr755qsCYLj48BIioqEFFALXX389CgsLMX36dN8v5uzs7EE/n5mZicbGRhgMBgiCgDVr1qCqqgo6nQ7z5s0bcBmTyQSz2Yz6+noolUqsXbt2BF/Hn06nQ1tb24DtREQUYAhMnjwZAPzO5hmKQqFAcXGxX1tCQsJVn9u3b5/v9cSJE1FZWRnQ+gNltVr9HmgO8OElRERXCujsoOXLl2PGjBmIjIzErbfeiuXLl4td1zXBh5cQEQ0toBAoLy/Hrl27EBERgbfeegtlZWVi13XNGI1GtLa24siRI2htbWUAEBFdIaDdQZ988gm2b98OAHjwwQeRlZUlalFERBQcAc0E+vr64PV6AWDA0zaJiGhsCmgmsGDBAuTk5OAnP/kJDh8+jAULFohdFxERBUFAIbB48WLMmTMHx48fx6JFiwI675+IiEa/IUPgjTfewO9+9zuUl5f7dgF98cUXAICVK1eKXx0REYlqyBDov03E1KlT/dp5TICIKDQMeWC4/06gn3/+OX7729/6/h06dCgoxRERkbiGnAnYbDZs3LgR33zzDd5//31f+0BX/xIR0dgzZAgYjUYYjUZUVFRg6dKlwaqJiIiCJKDrBOrr68Wug4iIJBDQKaITJ07Etm3bEB8fD4XiUm7MmTNH1MKIiEh8AYXAddddh5aWFrS0tPjaGAJERGNfwE8W+/LLL/HVV18hPj5+zD9lh4iILgkoBKqrq7Fnzx7MnDkTr7zyCubPn4+HHnpI7NqIiEhkAYXAnj17YLPZEB4ejt7eXhgMBoYAEVEICOjsIEEQfA+Kj4iIQEREhKhFERFRcAQ0E0hLS8Pjjz+OtLQ02O123HbbbWLXRUREQRBQCBQUFOCf//wnjh8/jvvuuw9z584Vuy4iIgqCgELg/PnzOHjwIE6cOIGzZ88iJSUFEydOFLs2IiISWUDHBJ588kkkJCTAbDZj8uTJeOqpp8Sui4iIgiCgmQAA5OTkAABuvfVWvPvuu6IVREREwRPQTGDq1KnYvXs3Ojo6sG/fPkRHR+PEiRM4ceKE2PUREZGIApoJHD9+HMePH8ff/vY3CIIAACgsLERYWBheffVVUQskIiLxBHzFcFdXF06ePInJkycjJiZG7LqIiCgIAtod9Pe//x0GgwEVFRXIzs7G22+/LXZdREQUBAHNBLZu3Ypdu3ZBo9HA6XTiwQcfxD333CN2bUREJLKAZgJhYWHQaDQAAK1Wi8jISFGLIiKi4AhoJqDT6VBaWopZs2bh008/hU6nE7suIiIKgoBmAiUlJZgyZQoOHTqEKVOmYPXq1WLXRUREQRDQTGDp0qV45ZVXxK6FKOhsNhssFgva29uh0+lgtVphNBqlLosoaAIKgaioKNTV1fk9Yzg+Pn7Qz3u9XhQVFeHo0aNQqVQoKSmBXq/3+0xnZycMBgPeeecdREZGwuPxwGw24/z589BoNCgrK+OpqCQqm80Gk8kEt9sNAGhra4PJZAIABgHJxg+GgNPpxKlTp7Bt2zZf2w9dJFZXV4eenh7U1taiqakJpaWl2Lhxo6//wIEDKC8vx7lz53xtNTU1SExMxGOPPYa9e/diw4YNeOaZZ0b6vYh+kMVi8QVAP7fbDYvFwhAg2RgyBF577TW88sorUCqVeOKJJ5CRkRHQSu12O9LT0wEAKSkpaG5u9utXKBSoqqrCfffd57fMww8/DADIyMjAhg0bhvVFiIarvb19WO1EoWjIENizZw/effddOJ1OPPXUUwGHgNPphFar9b1XKpXo6+vzPZ1s9uzZAy4TFRUFANBoNLhw4cKA63Y4HAHV8P95PJ4RLxuKOB5AbGwszpw5M2C7nMeGPxv+Qn08hgwBlUoFlUqFmJgY9Pb2BrxSrVYLl8vle+/1en0BEMgyLpcLEyZMGPBzSUlJAddxJYfDMeJlQxHHA1i3bp3fMQEAUKvVWLdunazHhj8b/kJhPOx2+6B9AZ0iCsB347hApKamoqGhAQDQ1NSExMTEgJapr68HADQ0NCAtLS3g7RGNhNFoRGVlJfR6PcLCwqDX61FZWcnjASQrQ/55/tVXXyE/Px+CIPhe9ysvLx90uczMTDQ2NsJgMEAQBKxZswZVVVXQ6XSYN2/egMvk5OSgoKAAOTk5iIiIGHL9RNeK0WiE0WgMib/2iEZiyBBYv36977XBYAh4pQqFAsXFxX5tCQkJV31u3759vtfjx4/HCy+8EPA2iIjofzdkCPzsZz8LVh1ERCSBgI8JEBFR6GEIEBHJGEOAiEjGGAJERDLGECAikjGGABGRjDEEiIhkjCFARCRjDAEiIhljCBARyRhDgIhIxhgCREQyxhAgIpIxhgARkYwxBIiIZIwhQEQkYwwBIiIZYwgQEckYQ4CISMYYAkREMsYQICKSMYYAEZGMMQSIiGSMIUBEJGMMASIiGWMIEBHJGEOASGI2mw1xcXFQKBSIi4uDzWaTuiSSkXCpCyCSM5vNBpPJBLfbDQBoa2uDyWQCABiNRilLI5ngTIBIQhaLxRcA/dxuNywWi0QVkdwwBIgk1N7ePqx2omuNIUAkIZ1ON6x2omtNlBDwer0oLCxEdnY28vLy0NbW5te/Y8cO3HvvvcjKysL+/fsBAN988w1uv/125OXlIS8vD9u2bROjNKJRxWq1Qq1W+7Wp1WpYrVaJKiK5EeXAcF1dHXp6elBbW4umpiaUlpZi48aNAICzZ8+iuroaO3fuRHd3N3JzczF79mx88cUXWLhwIZ599lkxSiIalfoP/losFrS3t0On08FqtfKgMAWNKCFgt9uRnp4OAEhJSUFzc7Ov7/Dhw7jtttugUqmgUqmg0+nQ0tKC5uZmHDlyBPfffz9iYmLwzDPP4MYbbxSjPKJRxWg08pc+SUaUEHA6ndBqtb73SqUSfX19CA8Ph9PpRFRUlK9Po9HA6XRi6tSpmDFjBu68807s3r0bJSUleOGFF65at8PhGFFNHo9nxMuGIo6HP47HZRwLf6E+HqKEgFarhcvl8r33er0IDw8fsM/lciEqKgozZ87E+PHjAQCZmZkDBgAAJCUljagmh8Mx4mVDEcfDH8fjMo6Fv1AYD7vdPmifKAeGU1NT0dDQAABoampCYmKir2/mzJmw2+3o7u7GhQsXcOzYMSQmJuKZZ57Be++9BwD48MMPkZycLEZpRER0BVFmApmZmWhsbITBYIAgCFizZg2qqqqg0+kwb9485OXlITc3F4IgYMWKFYiMjER+fj5WrVqFmpoajB8/HiUlJWKURkREVxAlBBQKBYqLi/3aEhISfK+zsrKQlZXl1z9lyhRUV1eLUQ4REQ2CF4sREckYQ4CISMYYAkREMsYQICKSMYYAEZGMMQSIiGSMIUBEJGMMASIiGWMIEBHJGEOAiEjGGAJERDLGECAikjGGABEBAGw2G+Li4pCcnIy4uDjYbDapS6IgEOUuokQ0tthsNphMJrjdbgBAW1sbTCYTAPDRlyGOMwEigsVi8QVAP7fbDYvFIlFFFCwMASJCe3v7sNopdDAEiAg6nW5Y7RQ6GAJEBKvVCrVa7demVqthtVolqoiChSFARDAajaisrIRer0dYWBj0ej0qKyt5UFgGeHYQEQG4FARGoxEOhwNJSUlSl0NBwpkAEZGMMQSIiGSMIUBEJGMMASIiGWMIEBHJGEOAiEjGGAJENKr0381UoVDwbqZBwOsEiGjU4N1Mg48zASIaNXg30+BjCBDRqMG7mQYfQ4CIRo3RdDdTuTxpjSFARKPGaLmbaf+xiba2NgiC4Ds2IUUQiH2gXJQQ8Hq9KCwsRHZ2NvLy8tDW1ubXv2PHDtx7773IysrC/v37AQCdnZ1YvHgxcnNz8eSTT+L7778XozQiGsVGy91MR8uxiaCEkSCC9957TygoKBAEQRA+++wzYenSpb6+r7/+Wli4cKHQ3d0tfPfdd77Xq1evFnbu3CkIgiBs2rRJqKqqumq9n3766Yhr+uKLL0a8bCjiePjjeFzGsRCEsLAwAcBV/8LCwoJah16vH7AOvV4/rPUM9btTlJmA3W5Heno6ACAlJQXNzc2+vsOHD+O2226DSqVCVFQUdDodWlpa/JbJyMjAoUOHxCiNiOgHjZZjE8E4UC7KdQJOpxNardb3XqlUoq+vD+Hh4XA6nYiKivL1aTQaOJ1Ov3aNRoMLFy4MuG6HwzGimjwez4iXDUUcD38cj8s4FsCyZctQWFgIj8fjaxs3bhyWLVsW1LGJjY3FmTNnBmy/VnWIEgJarRYul8v33uv1Ijw8fMA+l8uFqKgoX/u4cePgcrkwYcKEAdc90odd8EEZ/jge/jgel3EsLv2eufnmm2GxWNDe3g6dTger1Rr0YxPr1q3zu3gOuHSgfN26dcP6P7Lb7YP2ibI7KDU1FQ0NDQCApqYmJCYm+vpmzpwJu92O7u5uXLhwAceOHUNiYiJSU1NRX18PAGhoaEBaWpoYpRERBcRoNKK1tRVHjhxBa2urJFcsB+NAuSgzgczMTDQ2NsJgMEAQBKxZswZVVVXQ6XSYN28e8vLykJubC0EQsGLFCkRGRuKRRx5BQUEBduzYgeuuuw7l5eVilEZENKb0P/ZTLKKEgEKhQHFxsV9bQkKC73VWVhaysrL8+m+44Qa8/PLLYpRDRESD4MViREQyxhAgIpIxhgARkYwxBIiIZCxMEARB6iICNdS5rkRENLjBTrsfUyFARETXFncHERHJGEOAiEjGQj4EfujZBnLS29sLs9mM3NxcLFq0CP/4xz+kLmlUOH/+PObOnYtjx45JXYrkNm3ahOzsbNx777144403pC5HMr29vcjPz4fBYEBubm5I/2yEfAjU1dWhp6cHtbW1yM/PR2lpqdQlSWb37t2Ijo7G66+/js2bN2P16tVSlyS53t5eFBYWYty4cVKXIrmPPvoIn332GWpqalBdXY3//Oc/Upckmfr6evT19WH79u1YtmwZ1q9fL3VJogn5EBjq2QZy8+tf/xpPPPGE771SqZSwmtGhrKwMBoMBN954o9SlSO7gwYNITEzEsmXLsHTpUvz85z+XuiTJxMfH4+LFi/B6vXA6nb67IIei0P1m/zXUsw3kRqPRALg0Jo8//jiefPJJiSuS1q5duxATE4P09HRUVlZKXY7kurq6cPr0aVRUVODUqVN45JFH8O677yIsLEzq0oJOrVbj3//+N+bPn4+uri5UVFRIXZJoQn4mMNSzDeTozJkzeOCBB3DPPffg7rvvlrocSe3cuROHDh1CXl4eHA4HCgoKcPbsWanLkkx0dDTmzJkDlUqFqVOnIjIyEp2dnVKXJYmtW7dizpw5eO+99/D222/j6aefRnd3t9RliSLkQ2CoZxvIzblz57B48WKYzWYsWrRI6nIkZ7PZ8Nprr6G6uhpJSUkoKyvDpEmTpC5LMmlpaThw4AAEQUBHRwe+//57REdHS12WJCZMmOB70uHEiRPR19eHixcvSlyVOEL+T+KBnm0gVxUVFfjuu++wYcMGbNiwAQCwefNmHhQlAMBdd92FTz75BIsWLYIgCCgsLJTtcaPf//73WLVqFXJzc9Hb24sVK1ZArVZLXZYoeMUwEZGMhfzuICIiGhxDgIhIxhgCREQyxhAgIpIxhgARkYwxBGhM+OijjzBr1iycOXPG1/b8889j165dI17nqVOnkJWVdS3Ku8rFixfx0EMPIScnB99++62v/emnn8bdd9+NvLw85OTk4NFHH8XJkyeHvf5vvvkG77zzjm+d/dfCEA0XQ4DGjIiICPzpT3/CWDir+ezZs+jq6kJNTQ0mTpzo12c2m1FdXY2amhosXrx4RLfvOHr0KPbt23etyiUZC/mLxSh03HHHHfB6vbDZbLj//vt97adOncLKlSuxY8cOAEBWVhb+8pe/4M0330RbWxu6urrw7bffIjc3F++//z5OnDiBsrIy3HDDDejs7MTSpUvR2dmJuXPnYtmyZThz5gyeffZZdHd3IzIyEqtXr8bFixfxyCOPIDo6GhkZGViyZIlv+7t378a2bdugUqkQFxeH4uJiPPvss2htbUVhYSGKi4sH/U6zZs1CREQE2traoFKprtruTTfdhPLycjQ3N8PlciEhIQFr165FRUUFWlpaUFtbCwCora3Fli1b4HQ6UVRUhFtuuQVPPPEEnE4nPB4PzGYzbr/9dpH+Z2gs40yAxpSioiJs3boVra2tAX1+3LhxePnll/GrX/0K9fX1qKiogMlkwt69ewEAbrcb69atQ01NDQ4cOICWlhaUlZUhLy8P1dXVeOihh/D8888DuPTX/csvv+wXAF1dXfjrX/+Kbdu2oaamBlFRUaitrcVzzz2HadOmDRkA/a6//np0dXUNuF2n04kJEyagqqoK27dvR1NTEzo6OrB06VLccccdyM7OBgAkJyfj1Vdfxf33349du3ahvb0d586dQ0VFBcrLy+HxeIY50iQXnAnQmHLddddh1apVePrpp5GamjrgZ67cXTR9+nQAQFRUFKZNmwbg0r1g+m8Gduutt/ruEfPjH/8YJ06cwJdffolNmzZhy5YtEAQBERERAIDJkydDpVL5bevkyZOYNm2a7061P/3pT3Hw4MFh3Yb59OnTiI2NHXC7/TdxW7lyJdRqNdxuN3p7e69aR3JyMgDghhtugMfjwY9+9CMYjUasXLkSfX19yMvLC7gekheGAI05v/jFL/DBBx/gzTffhNlsRmRkJM6fP4+LFy/C5XLh1KlTvs/+0G2Qjx07BpfLhcjISBw+fBjZ2dmYOnUqFi9ejNTUVBw7dgyffPIJAEChuHriPHnyZBw7dgxutxtqtRoff/wx4uPjA/4ujY2NGDduHGJjYwfcbkNDA86cOYP169ejs7MTH3zwAQRBgEKhgNfrHfR7Hj16FC6XC5WVlfj6669hMBhw1113BVwXyQdDgMYki8WCf/3rXwCASZMmYfbs2Vi0aBF0Oh30en3A65k4cSJWrFiBzs5OLFiwANOmTUNBQQGKiorQ3d0Nj8cDi8Uy6PIxMTF47LHH8MADD0ChUECn0+GPf/zjkLekXrduHTZv3gyFQgGNRuN7atVA2508eTI2bNiArKwsqFQqTJkyBV9//TV0Oh2+/PJLbN26dcBtxMXF4aWXXsJbb72FiIgIPP744wGPCckLbyBHRCRjPDBMRCRjDAEiIhljCBARyRhDgIhIxhgCREQyxhAgIpIxhgARkYwxBIiIZOz/ABgwPePVP3I5AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "plt.style.use('seaborn-whitegrid')\n",
    "\n",
    "## proportion of different number of deaths\n",
    "plt.plot(np.arange(freq.size), freq/np.sum(freq), 'o', color='black')\n",
    "plt.xlabel(\"Number of Deaths\")\n",
    "plt.ylabel(\"Proportion\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Problem 1\n",
    "\n",
    "Confirm Lange's claim in question 11 that \"A Poisson distribution gives a poor fit to these data . . .\" by computing the MLE of the Poisson parameter $\\lambda$ (this can be done by hand using calculus) and plotting the resulting probability mass function over the observed data. They should be a bit discrepant (thus justifying the claim)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Poisson probability mass function is\n",
    "$$ f(x|\\lambda) = \\frac{e^{-\\lambda}\\lambda^x}{x!}$$\n",
    "\n",
    "The MLE maximizes the log likelihood. The log likelihood for the sample is\n",
    "\\begin{align*}\n",
    "\\log f(\\vec{x}) &= \\log \\prod_{i=1}^n f(x_i|\\lambda)\\\\\n",
    "&= \\sum_{i=1}^n \\log f(x_i|\\lambda)\\\\\n",
    "&= \\sum_{i=1}^n (-\\lambda + x_i\\log \\lambda - \\log x!)\\\\\n",
    "&= -n\\lambda + \\sum_{i=1}^n(x_i) \\log \\lambda - \\sum_{i=1}^n \\log x!\n",
    "\\end{align*}\n",
    "Differentienting, setting equal to $0$, and solving for $\\lambda$, we obtain\n",
    "$$ \\widehat{\\lambda} = \\bar{X} = \\frac{1}{n}\\sum_{i=1}^n X_i $$\n",
    "In other words, the MLE is the sample mean."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "import math as mt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYEAAAEBCAYAAACe6Rn8AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3XmcTfX/wPHXnc0wxr7LEplPI0sikSW7bH2VxI+yZ+vL2GIYsi/JrpBKUTQUohTVlyayhCJlfCJKqZGxZZv9/v74XJrdGHPvmZn7fj4e9zH3nvM5577nGPd9z+ecz/tjs9vtCCGEcE8eVgcghBDCOpIEhBDCjUkSEEIINyZJQAgh3JgkASGEcGOSBIQQwo15WR3AnTh48KDczyqEEJlQu3ZtW2rLc1QSAKhdu3amtgsPDycwMDCLo8m55HgkJcfjX3IsksoNx+PgwYNprpPuICGEcGOSBIQQwo1JEhBCCDcmSUAIIdyYJAEhhHBjkgSEEMKNSRIQQgg3JknAjZw/Dy+/DJUrw/jxpa0ORwiRDUgScAOHDkHfvnDPPRAcDDYbbNhQiC++sDoyIYTVJAnkUrGxsG4dNGoEtWpBaCj06gU//mge5cvH8MILEB1tdaRCCCtJEshlzp6FqVOhYkXo0gX+/BPmzYMzZ2DpUnjgAfD1hfHjIzh+HObMsTpiIYSVJAnkEvv3Q48eUL48vPQSVK8On3wCx4/D8OFQqFDS9g0bXuPpp2HaNDh1ypqYhRDWkySQg0VHw3vvQb16ULcubNwI/fvDsWOwdSu0awce6fwLz58Pnp4wdKjrYhZCZC+SBHKgP/803/YrVIDnnoOLF2HRItPls3gxKJWx/dxzD0yebM4YNm92bsxCiOxJkkAOYbfD7t3Qtav58J82DerUMd/4w8NhyBAoUODO9zt0qLlOMHQoXLuW9XELIbI3SQLZXFQUvPOO+cBv0MB86A8ZYvr6P/kEWrdOv8vndry9zQXj336D6dOzLGwhRA4hSSCbOn0axo2DcuWgd2/T/79smenymTfPDPjKKo0aQc+e5k6hY8eybr9CiOxPkkA2YrfDV19Bp05w771mdG/DhvC//8GRIzBgAPj5Oee9Z882+37hBROHEMI9SBLIBq5fh+XLoWZNaNrUJIJRo+DkSXPHT7NmZpSvM5UoATNmwPbtZmCZEMI9OGWOYaWUB7AEqAlEA/201icSrR8OdHW8/FRrPVkpZQP+AI47lu/RWo91RnzZxalT8Npr8NZbcOmSSQJvvgndukHevK6Pp39/WLECRoyAtm2hYEHXxyCEcC1nTTTfEfDVWtdXStUD5gL/AVBKVQK6A48AdmCnUmojcB34TmvdwUkxZQt2O3z5pbmV85NPzEXdp54yF3sbNnT+N/70eHqai8R168LEibBggXWxCCFcw1ndQQ2BrQBa671AnUTrfgce11rHa60TAG8gCqgNlFVK7VBKfapURu92zznCwqBqVWjVCvbuNRd+f/313xo/ViaAm+rUgYEDTZI6dMjqaIQQzmazO+EqoFLqTWC91vozx+vTQCWtdVyiNjbgFcBfaz1AKdUYKKm1/kAp1RCYr7V+OPF+Dx48aM+XL1+mYoqKisLX1zeTv1HW6NTpXi5e9CQo6ByPP/4PefJYdwU2veNx+bIH7dpVpnz5GN5777e7ugU1p8gOfx/ZhRyLpHLD8bh+/Tq1a9dO9Wums7qD/gH8E732SJYAfIEVwBVgsGPxASAOQGu9SylVVill01on+aQMDAzMVEDh4eGZ3jYraG0GdS1YAEFBZYAyWbZvu93O1atXiYyMJDIyknPnzqX6PPHrNm3a8O6776a5z3nzoFcvL/buDaRv3ywLNduy+u8jO5FjkVRuOB4HDx5Mc52zksA3QAdgneOawJGbKxxnAJuA7VrrlxNtMxE4D8xWStUETidPADnZ2rWmu6dz59u3jYmJ4fz58xn+QI+MjCQ6jZrQXl5eFC9enGLFilGsWDFq1qzJ33//zXvvvccLL7xAvXr1Ut2uRw9zwXrMGOjYEYoWvZvfXgiRXTkrCWwEWiqldgM2oLdSagRwAvAEHgPyKKXaONqPBWYB7yml2mHOCHo5KTaXs9vh/fehQYN4DhzYcttv7JcvX05zX4UKFbr1oV6+fHlq16596wM+8Yf9zecFChTAluxiw5UrV6hcuTJDhgxh3759eKTS32OzwZIl8OCDMHasuYVVCJH7OCUJOC74Dky2OPFY1LQ62No5Ix6rHTliRuIGBr7Gf/4TdGu5r69vkg/uypUrp/uBXqRIEby9ve86Hn9/f0aOHElwcDDvvPMOffr0SbVdtWowbBjMnQt9+phqpUKI3MVZZwIikdBQ8PS0Ex4+lXHjxtG/f3+KFStGvnz5UnxLd5UOHTqwadMmgoODeeqppyiUfMIBh4kTTfyDBpk5C7zkL0aIXMUN7vuwlt1uPkSLFPmewoXjCQ4OpkKFCvj5+VmWAABsNhuLFy8mMjKSKVOmpNnO399czD50yIwhEELkLpIEnGz/fjMy+Ny5RQQFBeHv73/7jVykdu3a9OvXj8WLF3P06NE023XqZMY2jB8Pf/3lwgCFEE4nScDJQkPBwyMWP78vGTJkiNXhpDB9+nT8/PwICgoirTEjNhu8+qopaz1qlIsDFEI4lSQBJ0pIgNWr40hI2MILL3SnSJEiVoeUQvHixZkyZQpffvklmzZtSrNdlSoQHAxr1sCOHS4MUAjhVJIEnGjXLvj7by+8vTcwYsQIq8NJ06BBg3jggQcYPnw4N27cSLNdcLApcT14MMTEuDBAIYTTSBJwojfeuAJco0+f4pQsWdLqcNLk7e3NwoUL+fXXX5k7d26a7fLmNd1Cx46ZEcVCiJxPkoCTxMXBhx+CzfYJ48YF3X4DizVv3pxOnToxY8YMfv/99zTbtW0LTz4JU6eaKSmFEDmbJAEnWb/+IlFR/jRrdo7y5ctbHU6GzJkzB7vdzosvvphuu5slpocNc0FQQginkiTgJNOnnwAuM39+a6tDybCKFSsyZswY1q5dS1hYWJrtypeHl16Cjz6CLVtcGKAQIstJEnCCv/66wJEjVahY8RDVq1exOpw7Mnr0aMqXL8/QoUOJi4tLs93w4RAYaCbDSedashAim5Mk4AQjRmwDCjF6dDmrQ7lj+fLlY+7cufzwww8sT6dqnI+PKTB36hTMnOnCAIUQWUqSQBa7cuUK69d74+39D/36VbI6nEzp1KkTTZs2Zfz48Zw/fz7Ndk2aQPfu8PLLcPx4ms2EENmYJIEstnjxCmJj29C+/Q2yoOCnJWw2GwsXLuSff/5hwoQJ6badMwd8feG//zV1koQQOYskgSwUFRXF7Nk/AX4MHZp9xwVkRPXq1Rk8eDCvv/46h9KZbLhUKZg2DT7/3NwSK4TIWSQJZKEVK1Zw+fLjFC0aTaNGVkdz9yZPnkzhwoUZOnRomnWFwJSZrlXL3DJ65YoLAxRC3DVJAlkkNjaWWbOWYrO1o3t3Hzw9rY7o7hUuXJgZM2awc+dO1q5dm2Y7Ly9TZvqvv2DyZBcGKIS4a5IEssjq1av5/fda2O156NrVunkCslrfvn2pVasWo0aN4tq1a2m2e+QR6NfPDCQ7ciTNZkKIbEaSQBaIj49n5syZ+Ps/T4UK9lw1DaOnpyeLFy/mzJkzzLzNvaAzZ0KhQqbAnFwkFiJnkCSQBTZs2MDPP0dy/XoDunSxYeGEYU7RoEEDunfvziuvvMIvv/ySZruiRWH2bFM9ddUqFwYohMg0SQJ3yW63M336dEqWHEx8vAddu1odkXPMnj0bb29vRo4cmW67Xr2gfn148UW4cME1sQkhMk+SwF369NNPOXz4MEWKDCIgAB580OqInKNMmTJMmDCBTZs2sW3btjTbeXiYi8Tnz0NIiAsDFEJkiiSBu3DzLKBs2TocO1aarl3JdV1BiQ0bNoz77ruPoKAgYtKZVaZmTRg6FF5/3cyxLITIviQJ3IWvvvqKPXv28Oij87DbbXTpYnVEzpUnTx4WLFiA1ppXX3013baTJ5uBZIMGQXy8iwIUQtwxSQJ3wVwLKMnvvzegRg2oWtXqiJyvXbt2tG3blkmTJhEREZFmuwIFzOxjBw+aMwIhRPYkSSCT9u3bx//+9z/69JnM3r2594JwaubPn09UVBTjxo1Lt12XLtC8OYwbB2fPuig4IcQdkSSQSTNmzKBw4cLky9cTINd3BSUWEBDA8OHDefvtt/n222/TbGezwWuvwfXrMHq0CwMUQmSYJIFMOHLkCJs3byYoKIiNG32pWxcq5cyq0Zk2fvx4SpUqxZAhQ0hISEiznVLmdtFVq+Drr10YoBAiQ7ycsVOllAewBKgJRAP9tNYnEq0fDtzsQPlUaz1ZKZUXeA8oAVwBemqtzzkjvrs1Y8YM8ufPT5s2QUyaZPq+3Y2/vz+zZ8+mR48erFq1il69eqXZNiQEVq82I4m//54cW2JbiNzIWWcCHQFfrXV9IBiYe3OFUqoS0B14FKgPtFJK1QAGAUe01o2AVcB4J8V2V44fP866desYPHgw27YVwmaDZ56xOiprdO/enfr16xMcHMzly5fTbJcvHyxaBD/9BAsXujBAIcRtOSsJNAS2Amit9wJ1Eq37HXhcax2vtU4AvIGoxNsAnwEtnBTbXZk1axY+Pj4MHz6C99+HRo2gbFmro7KGh4cHixYt4u+//2bq1Knptn3iCejQASZNgj/+cE18Qojbc1YSKAAk/moYr5TyAtBax2qtI5VSNqXUHOB7rfXPyba5AhR0UmyZdvr0aVatWkW/fv04d64k4eG41V1BqalTpw59+/Zl4cKFHDt2LN22CxdCQoKZpF4IkT045ZoA8A/gn+i1h9Y67uYLpZQvsALzYT84lW38gUup7Tg8PDxTAUVFRWV625umTZsGQMeOHXn11Ug8PYtSvfpxwsNz3miorDgeN/Xs2ZO1a9fSr18/li9fji2dYdPPP1+URYtKsHz5aRo1Srs0tatl5fHI6eRYJJXrj4fdbs/yR0BAQKeAgIB3HM/rBQQEfJZonS0gIGBbQEDAmGTbjAwICJjkeN41ICBgafL9HjhwwJ5ZR48ezfS2drvdHhERYff19bX37t3bnpBgt1eqZLe3anVXu7TU3R6P5BYsWGAH7Js2bUq3XVSU3R4QYLdXrmy337iRpSHclaw+HjmZHIukcsPxcHx2pvp57azuoI1AlFJqNzAfGK6UGqGUegJz0fgxoI1S6ivHoz6wFHhAKbUL6A9kqzmq5s+fT0xMDMHBwRw4ACdPSldQYoMHD6Zq1aoMHz6cqKioNNvlyWPGDvzyiyk7LYSwllO6gxwXfAcmW5y4w9g3jU07OyOeu3Xx4kWWLFlC586dCQgIYORIc5vjk09aHVn24e3tzcKFC2nZsiXz5s1LdzRxixZmcN2MGdC9O1Su7MJAhRBJyGCxDFi8eDFXrlxh3LhxJCTA2rXQpo2ZRUv8q0WLFjz11FNMnz6dP25zC9C8eSaRjh3rouCEEKmSJHAbV69eZeHChXTo0IEaNWrwzTdw5ox7lYm4E3PnziUhIYHRt6kTUaaMuUvogw/g0CEXBSeESEGSwG0sW7aMCxcuEOKYISU0FPLmNfe9i5QqVqzI6NGjef/999m5c2e6bUeMMGdTL73kouCEEClIEkhHVFQUc+fOpXnz5jzyyCPExZlvru3bQ/78VkeXfY0ZM4Zy5coxZMgQ4tOZTKBQIVNX6OOPYd8+FwYohLhFkkA6VqxYQURExK2zgB074Nw5uSvodvLly8fcuXM5fPgwb7zxRrpthw6FYsVgwgQXBSeESEKSQBpiY2OZPXs29erVo0mTJoC5IOzvby4Ki/Q9/fTTNGnShJCQEC6kM+N8/vzm4vAXX0BYmAsDFEIAkgTStGbNGn777TdCQkKw2WzExMD69dCxo7kmINJns9lYtGgRly5dYsJtvuYPGgSlS5uzAbvdRQEKIQBJAqmKj49n5syZ1KxZk3bt2gHw+edw6ZJ0Bd2J6tWrM3jwYJYtW8bhw4fTbJc3L4wfDzt3mjMCIYTrSBJIxYYNG9BaM27cuFt1cEJDoUgRM9BJZNzkyZMpXLgwQ4cOxZ7O1/y+faFCBZMM5GxACNeRJJCM3W5n+vTpKKXo1KkTYKZH3LQJOnUCHx+LA8xhihQpwvTp0/n6669Zt25dmu3y5DG3iu7fb+4WEkK4hiSBZD799FMOHz5McHAwnp6ejmVw9ap0BWVWv379ePDBBxk1ahTXrqVdObRHD7jvPnNtIJ0ZK4UQWUiSQCI3zwIqVKhA9+7dby0PDYWSJeGxxywMLgfz9PRk8eLF/PHHH8yaNSvNdl5eMHky/PADfPihCwMUwo1JEkgkLCyMPXv2MHr0aLwdE+H+8w9s2WKmkHScGIhMaNiwId26deOVV17h5MmTabbr0gUeeMB0DcXFpdlMCJFFJAkkMn36dEqWLEnv3r1vLdu8GaKipCsoK8yePRsvLy9GjhyZZhtPT5gyBbSGNWtcGJwQbkqSgMO3337Ll19+yciRI8mbaCBAaCiULw/16lkYXC5RtmxZQkJC+Oijj/j888/TbPfkk1CrlpmPODbWdfEJ4Y4kCThMnz6dwoULM3Dgv9MgXLgA27aZLgoPOVJZYvjw4VSuXJmgoCBi0/iEt9lg2jQ4dQreftvFAQrhZuSjDThy5AibN28mKCgIf/9/p0besMH0S0tXUNbx9fVl/vz5HDt2jFdffTXNdm3aQP36MHWq6Y4TQjiHJAFgxowZ5M+fnyFDhiRZHhoKVaqYrgmRddq3b8/jjz/OpEmTOHv2bKptbp4N/PEHLF/u4gCFcCNunwSOHz/OunXrGDx4MEWKFLm1PCLCVA3t2tV8IImsY7PZWLBgATdu3GDYsGFptmvWDJo2NdNQpjO8QAhxF9w+Cbz88sv4+PgwYsSIJMs//NAMWJKuIOdQSjFhwgRCQ0P54IMP0mw3dSqcPWsmpxdCZD23TgK///47q1atol+/fpQsWTLJutBQqF4dqla1KDg3EBwcTJ06dRg0aBARERGptmnQwFwfePllM2ZDCJG13DoJvPLKK9jtdl588cUky0+fhm++kXmEnc3b25tVq1Zx9epVBgwYkGaBualTzZ1aCxa4OEAh3IDbJoGzZ8/yxhtv8Nxzz1G+fPkk627WOZMk4HyBgYHMmDGDzZs3s3LlylTb1K5txg7MnWuSgRAi62QoCSilWiqlPlNKbb/5cHZgzjZ//nxiYmIIDg5OsS40FOrUMcXMhPMNGzaMxo0bExQUxOnTp1NtM3kyXLkCc+a4ODghcrmMngnMB+YCgxI9cqyLFy+yZMkSOnfuTEBAQJJ1x4/DwYNyQdiVPDw8ePvtt4mPj6dv374kpFJCtHp182+ycCH8/bcFQQqRS2U0CZzWWn+pE3FqVE62ePFirly5wrhx41KsW7vW/HzmGRcH5eYqVarE3Llz+fLLL1m6dGmqbSZNMgPH0ilEKoS4QxlNAn8rpZYppQYopforpfo7NSonunr1KgsXLqRDhw7UqFEjxfrQUGjYEMqVsyA4N9e/f39at27N6NGjOX78eIr1AQHQsycsWWIGkQkh7l5Gk8Ap4C+gFFDa8ciRXn/9dS5cuEBISEiKdT/+CD/9JF1BVrHZbLz11lv4+PjQq1cv4uPjU7R56SUzfmP6dAsCFCIXylAS0FpPBg4AN4BDjtc5TlRUFHPmzKF58+Y88sgjKdaHhppCcU8/bUFwAjCVRhcvXszu3buZO3duivUVK8Lzz8Obb5oCc0KIu+OVkUZKqZlAFWAX0FMp1UhrPSqd9h7AEqAmEA3001qfSNamOLAbqK61jlJK2YA/gJv9AHu01mPv9BdKz9tvv01ERASrV69Osc5uN9cDmjUzs4gJ63Tv3p2NGzcyYcIE2rZtS7Vq1ZKsDwmBFSvMvANSZVSIu5PR7qDGWuuntdYLgE5Ao9u07wj4aq3rA8GYO4tuUUq1Bj4HEn/cVga+01o3cTyyNAHExsby8ssvU69ePZo2bZpi/XffwYkT0hWUHdhsNpYuXUrBggXp0aNHipLTZcrA4MGwapWZfEYIkXkZTQLejm/3ADYg9aGd/2oIbAXQWu8F6iRbnwC0ABIP/akNlFVK7VBKfaqUUhmMLUO2bNnCb7/9RkhICLZUKsKFhoK3txmUJKxXokQJXn/9db7//numTZuWYv2YMZA3r7ljSAiReRnqDgLWAt8opfYCjzhep6cAcDnR63illJfWOg5Aa/0FmCJiifwFzNRaf6CUagi8BzycfMfh4eEZDDnRm8fHs3z5cpRSVKpUKcU+EhLgvffuo0GDKM6e/YM0qhvnKlFRUZk6lq50//3388QTTzB9+nSqVauWoluoe/fiLF9ejC5dTqJU9F29V044Hq4ixyKpXH887HZ7hh4BAQHVAgICng4ICKiWgbbzAgICnkn0+o802v0aEBDg63ieLyAgwCfRuj8DAgJsidsfOHDAnhlbtmyxA/a1a9emun7XLrsd7Pb33svU7nOko0ePWh1Chly4cMFetmxZe2BgoP3GjRvJ1tntBQva7R073v375JTj4QpyLJLKDcfD8dmZ6ud1ut1BSql+jp8zgW7AQ0A3pdSM2+SWb4C2jm3rAUcykI8mAsMc29TEDFC7XbdThlSuXJkBAwbQqVOnVNeHhoKvLzzxRFa8m8hKhQsX5q233iI8PJzx48cnWwcjR8JHH8GBAxYFKEQOd7trAr87fh4DdKLHsdtstxGIUkrtxpScGK6UGqGUSu9jdhbwmFIqDJgH9LrNe2SYUoqgoCA8PT1TrIuLMwXj2reHRDNLimykdevWDBw4kHnz5rFz584k64KCoGhRmDDBouCEyOHSvSagtd7mePqw1vq/N5crpVYBq9LZLgEYmGxxisShta6Y6PlFoN3tQ85aYWGmFo3cFZS9vfLKK2zbto1evXpx+PBh8ufPD0CBAuYi8ejRsGuXGe0thMi423UHvaCU+gt4Xin1p+PxF1DWNeE5X2go5M8PbdtaHYlIT/78+Vm5ciWnTp1KMf/DCy+YsR3jx5vxHkKIjEs3CWitX9NalwYma63LOB6ltdbNXRSfU8XEwPr10LGjud1QZG+NGjVixIgRLFu2jG3btt1ani+fGUAWFgbbc3yRcyFcK6PjBHLl9+QvvoCLF6UrKCeZNm0agYGB9O3bl0uXLt1a3r+/KfonZwNC3JmMJoELSqkgpdTjSqlWSqlWTo3KRUJDzR0mLVtaHYnIKF9fX1atWkVERARDhw69tTxPHnNxeO9e+PRTCwMUIofJaBI4DzwIdAH+D8jx351v3DC3Fj71FPj4WB2NuBN16tQhJCSEd999l40bN95a3qsXVKpkkkEq89IIIVKR0SqivTH1fz4DFmqt+zg1Khf49FO4elW6gnKqkJAQatWqxYABAzh37hxgyn5MmgTffw+JcoMQIh0ZnWN4CPAG8CiwXCmVZgXRnCI0FEqUgCZNrI5EZIaPjw+rVq3i8uXLDBw4ELvjQkC3bnD//WbegVSmIxBCJJPR7qBuQCOt9TCgAaZbKMe6cgU++QQ6dwavjFZPEtlOtWrVmDp1Khs2bGDNmjUAeHqaEtNHj5pEL4RIX0aTgC1R8bdYIPY27bO1zZvNXLXSFZTzjRw5kvr16/Pf//6XM2fOANCpE9SsabqGYnP0X6oQzpfRJLBLKfWh4w6hDzG1gXKs0FC45x549FGrIxF3y9PTk5UrVxITE0O/fv1MQSwPmDrVzA+xKs1x7UIIyPiF4VHA25gyEyu01i/eZpNs68IF2LYNunQxU0mKnK9KlSrMnj2brVu38sYbbwCmFlTduqZrKPruqkwLkatl9MJwCaAV0BJoppQq7NSonGjjRtNFIF1BucugQYNo3rw5I0aM4OTJk9hsMG0anD5t5iMWQqQuo9+F1wLhwBjgJPCu0yJystBQqFwZate2OhKRlTw8PFixYgWenp707t2bhIQEWrSAxo1NMrh+3eoIhcieMtwhorVeprU+rLVeAuR3YkxOc/asqS3TtSukMsOkyOHKly/PwoUL+frrr1m4cCE2m7k2EBEBS5daHZ0Q2VNGk8AxpVR3pVQZpVQH4LxSKkApFeDM4LLahx+akaTSFZR79ezZkw4dOjB27FjCw8Np3BhatYJZs8ytwUKIpDKaBO4H+mHm/R0BFAFeB5Y5KS6nCA2FBx6AZFPVilzEZrOxfPly8ufPT8+ePYmLi2PqVIiMhEWLrI5OiOwno3cHNQWeBoKBzlrrpo5HM6dGl4UiIrzYtUvOAtxBqVKlWLp0Kfv372fWrFnUrWumDn3lFVM1Vgjxr4zeHdQZ2A2MA/YqpZ51alROsHVrAcDcGipyv86dO9O1a1cmT57MoUOHmDIFLl+GuXOtjkyI7CWj3UEjgNpa645ALSDIeSE5x2efFaB2bahSxepIhKu8+uqrFCtWjOeee47774+mSxdYsAAc9eaEEGQ8CSRora8CaK2vAFHOCynr/fILHDmSV7qC3EzRokV58803+fHHH5k0aRKTJpkS4i+/bHVkQmQfGU0Cvyil5iql/qOUmgv84sygstqePeDpaeeZZ6yORLhau3bt6Nu3L7Nnz+bixT089xy89hr8+afVkQmRPWQ0CfTDDBJr6fj5vNMicoKnn4ZNm05SvrzVkQgrzJs3j3LlytGzZ09GjbpOXBzMmGF1VEJkDxktpPyJ1jrHTinp6wuVKsVYHYawSIECBXj77bdp1qwZy5cH07fvIpYvhxdfhAoVrI5OCGtl9EzgkqMr6P6cOEhMiKZNmzJ06FAWL15M06bf3Ko0KoS7u20SUEoVAO4FhgFLyYGDxIQAmDlzJgEBAYwe3Y3evaN55x04ftzqqISwVrpJQCn1X+AwUBh4OScOEhPipnz58rFy5Ur++OMPLl0KJk8emDzZ6qiEsNbtzgS6AQqoRw4cGyBEcvXq1WPMmDGEhi6gTZsTrFkDP/1kdVRCWOd2SSBKax2jtY4EfFwRkBDONnHiRGrUqMHOnU/g55fAxIlWRySEde5kbi0pvixyhTx58rBq1SouXjxBhQobWL8evv8+r9VhCWGJ290i+oBSag2cpt8NAAAYrElEQVQmAdx8DoDWultaGymlPIAlQE0gGuintT6RrE1xTD2i6lrrKKVUXkyV0hLAFaCn1loG+AunqFmzJhMnTmT8+D4UL/44I0eWpVUrKFbM6siEcK3bnQk8w793A918fvORno6Ar9a6PqbyaJKyXUqp1sDnQMlEiwcBR7TWjYBVwPgM/g5CZMqYMWOoWzeQ6OgnOH/eg27dID7e6qiEcK10zwS01mGZ3G9DYKtjH3uVUnWSrU8AWgAHk20z2/H8M2BCJt9biAzx8vJi5cqV1K1bl/z5x/LFF7OZPNlMTi+Eu8joiOE7VQC4nOh1vFLKS2sdB6C1/gJAKZXWNleAgqntODw8PFMBRUVFZXrb3EiOx7+WL1/O888/j5/fQ0yd2pUyZX7nsceuWh2WZeRvI6ncfjyclQT+AfwTvfa4mQAyuI0/cCm1RoGBgZkKKDw8PNPb5kZyPP4VGBiIzWajf/8gfHyqMWZMIIcOeXLvvVZHZg3520gqNxyPgwcPprnuTu4OuhPfAG0BlFL1gCN3sg3QBtjpnNCESKlGjRp8+eXH5MnTnatXr9KhQzRROapguhCZ46wksBGIUkrtBuYDw5VSI5RST6SzzVLMHUi7gP6AjOUULvXwww+zfftb5M07gJ9+ykPPnjIzvcj9nNIdpLVOAAYmW3wslXYVEz2/DnR2RjxCZFSdOnX4+uvRNGgwh3XrRlGr1jmCg4tbHZYQTuOsMwEhcqyHHnqInTub4+W1g7Fj/dm8+Q+rQxLCaSQJCJGKOnVq8fnnxfDwiOSppxI4ePCU1SEJ4RSSBIRIQ9Om1Vm58gbx8aVo0OAkP/984vYbCZHDSBIQIh3PPluFcePOER3dnDp1NnBcJiAQuYwkASFuY9q0srRvf4krV0ZRr95LaK2tDkmILCNJQIjbsNkgNLQQVarEcPHiqzRq1J1jx1Lc7CZEjiRJQIgM8PODLVt88fMrxMWLy3nssVa5upSAcB+SBITIoCpVYNUqT+LiHuLKlSk0adKEo0ePWh2WEHdFkoAQd+DJJ2H0aLhxoxfR0c/QpEkTfvzxR6vDEiLTJAkIcYemT4cmTSA6eiF2e3WaNWvGkSMZKY8lRPYjSUCIO+TlBaGhULiwB35+n+HlVZRmzZrxww8/WB2aEHdMkoAQmVCyJHzwAZw540O1avvJk8eXZs2acfjwYatDE+KOSBIQIpMaNIBXXoEvvshP9+6HyJcvH82aNePQoUNWhyZEhkkSEOIuBAXBM8/AnDlFmTlzL/nz56dZs2Z89913VocmRIZIEhDiLths8OabEBAAI0aUITR0JwUKFKBFixbpzuYkRHYhSUCIu+TvDxs2wLVrMGpUeb744qtbieDAgQNWhydEuiQJCJEFAgNhxQrYvRtee60iYWFhFCpUiBYtWrB//36rwxMiTZIEhMgizzwDw4bBwoWwd28FwsLCKFKkCC1atGDfvn1WhydEqiQJCJGFZs82dw317QtXr5YnLCyMYsWK0apVK/bu3Wt1eEKkIElAiCzk7Q3r1pmCc506QaFC5QgLC6N48eK0atWKPXv2WB2iEElIEhAii5UpA2vXws8/mzOCsmXvISwsjFKlStGqVSu++eYbq0MU4hZJAkI4QZMmMHOmGVW8cCGULVuWHTt2UKZMGR5//HF27dpldYhCAJIEhHCaF1+Ejh3Nz127/k0EZcuW5fHHH2fnzp1WhyiEJAEhnMVmg3fegXvvNXcORURAmTJl2LFjB+XKlaNNmzaEhYVZHaZwc5IEhHCiggVh/Xq4dAm6doW4OChdujQ7duygfPnytG3blq+++srqMIUbkyQghJNVrw7Ll0NYGIwbZ5aVKlWKHTt2ULFiRdq2bcv27dutDVK4LUkCQrjAs8/C4MGm6uiGDWZZyZIl2bFjB5UqVaJ9+/b873//szZI4ZYkCQjhIvPmQd260KuXuX0UoESJEuzYsYP77ruPdu3asWDBAhISEiyNU7gXSQJCuEiePPDhh+DjYwaSXbtmlhcvXpzt27fTokULhg8fTtOmTTl58qS1wQq34ZQkoJTyUEotU0rtUUp9pZS6L9n655VSB5RSe5VS7R3LiiilIh3tv1JKBTkjNiGsVK4cvP8+/PQTDBgAdrtZXqxYMT7++GNWrFjBoUOHqFGjBkuWLJGzAuF0zjoT6Aj4aq3rA8HA3JsrlFKlgKFAA6A1MFMplQd4CHhfa93E8VjopNiEsFTLljBlCqxeDUuX/rvcZrPRu3dvfvzxRxo0aMALL7xAq1atOH36tHXBilzPWUmgIbAVQGu9F6iTaF1d4ButdbTW+jJwAqgB1AYeUkqFKaU+UEqVdlJsQlhu3Dho185UHU1eV65cuXJs3bqVZcuWsXfvXqpVq8Zbb72F/eZpgxBZyMtJ+y0AXE70Ol4p5aW1jktl3RWgIHAMOKi1/lIp1R1YDDydfMfh4eGZCigqKirT2+ZGcjySsuJ4jB/vweHD9/LkkzY+/PAURYrEJ1nfuHFjNm7cyPjx4+nXrx8rV65kypQplCxZ0qlxyd9GUrn9eDgrCfwD+Cd67eFIAKmt8wcuAfuA645lG4Epqe04MDAwUwGFh4dnetvcSI5HUlYdj82boX59mDQpgK1bwdMz6frAwECaN2/Oa6+9xpgxY3jyySdZtGgRzz77LDabzSkxyd9GUrnheKQ31amzuoO+AdoCKKXqAUcSrfsWaKSU8lVKFQQCgR+BN4FOjjbNAZmgVeR6tWrBkiXw5ZcwcWLqbTw8PBgyZAiHDx+matWq9OjRgyeffJKIiAjXBityJWclgY1AlFJqNzAfGK6UGqGUekJrHQEsAnYC24EQrXUU5gLyIKXUV8BAQO4OEm6hTx/o1w+mT4cePeDcudTbValSha+//po5c+awdetWqlWrxtq1a10brMh97HZ7jnkcOHDAnllHjx7N9La5kRyPpKw+HtHRdvv48Xa7t7fdXqSI3b5ihd2ekJB2+6NHj9offvhhO2Dv3Lmz/e+//86yWKw+FtlNbjgejs/OVD9XZbCYENmAjw9MnQqHDkHVqubsoGlT0Dr19oGBgezevZsZM2bw0UcfUa1aNTZu3OjaoEWuIElAiGykalVTaO6NN+DwYahRAyZNgujolG29vLwYO3YsBw8epGzZsjz11FM8++yzXLhwweVxi5xLkoAQ2YyHh7lGcOwYPP00TJ5skkFaFaerV6/Ovn37mDRpEmvXrqVatWps2bLFpTGLnEuSgBDZVMmSZlTx1q0QG2u6h/r0gfPnU7b19vZm4sSJ7Nu3j6JFi9K+fXv69OnD5cuXUzYWIhFJAkJkc61bw48/QnAwvPsu3H+/+ZnaAOKHHnqIAwcOMHbsWFauXEm1atX4/PPPXR+0yDEkCQiRA+TLZyau/+47qFLF3ErasiUcP56ybZ48eZgxYwZ79uwhf/78tG7dmoEDB3LlyhXXBy6yPUkCQuQg1aubSeuXLoUDB8zradMgJiZl27p16/Ldd98xcuRIli9fTo0aNdixY4frgxbZmiQBIXIYDw8YOBDCw+E//4EJE+DBB01ySC5v3rzMmTOHnTt34uXlRbNmzRg6dCjXbk5mINyeJAEhcqjSpWHtWtiyBa5fh0aN4Pnn4eLFlG0bNGjAoUOHGDJkCIsXL+bBBx/km2++cX3QItuRJCBEDte2rZmkZtQoePttc+H4/fdTXjj28/Nj0aJF7Nixg7i4OBo1asSoUaO4ceOGNYGLbEGSgBC5gJ+fmcT+wAGoUAG6dYM2bSC1WSqbNGnCDz/8QP/+/Zk7dy4PPfQQ3377reuDFtmCJAEhcpEHH4Q9e2DxYti9Gx54AGbNMuMMEvP392fZsmVs27aNq1evUr9+fUJCQohObWiyyNUkCQiRy3h6wn//C0ePmq6isWPhoYdMckiuVatWHDlyhJ49ezJjxgzq1KnDkSNHUjYUuZYkASFyqXvugfXrYdMmuHQJGjSAwYPN88QKFSrEihUr+Pjjj4mMjKRLly7Uq1eP119/nUvJG4tcR5KAELncE0+Ys4KgIHj9dQgMhHXrUl44bt++PUePHuXFF1/k6tWrDBw4kNKlS9OtWzc+//xz4uPjU38DkaNJEhDCDfj7w/z58O23UKYMdOkC7dvDr78mbVe4cGF69+7NkSNH2L9/P3369OGzzz6jdevWVKxYkZCQEI6nNkxZ5FiSBIRwI7Vrw759JiGEhZkLx3PmQFxc0nY2m406derw2muv8ddff7F27VqqV6/OrFmzCAgIoGHDhrz11lv8888/1vwiIstIEhDCzXh5wbBhpouoeXN48UWoU8ecJaTG19eXZ555hk8//ZTTp08za9YsIiMj6devH6VLl6ZHjx5s376dhIQE1/4iIktIEhDCTZUvby4ar19v5jWuVw+GDoWrV9P+WChbtixjxowhPDycPXv28Oyzz7Jp0yaaN29OpUqVmDhxIqdOnXLhbyHuliQBIdyYzQZPPWXqEL3wArz6KjRrdh/t28Ps2abrKPkYA7Od7dYdRBEREaxevZqAgACmTp1KpUqVaNKkCStXrpQaRTmAJAEhBAUKmAFm+/ZBmzb/cOIEjBljzg4KFTJlq6dOha+/hqiopNvmzZv31h1Ev/76K9OmTePMmTP06tWLUqVK0adPH3bu3Ik9tQkQhOUkCQghbnn4YZg8OYJjxyAiAj74wMxm9vffMHEiPPYYFCwIjRub6qVffAFXr/67ffny5QkJCeHnn39m586dPPPMM3zwwQc0btyYKlWqMG3aNE6fPm3dLyhSkCQghEhVyZJmjuPFi82k95GR5hrCkCHmbGDmTGjVCgoXNmcMo0fDJ5+YwWg2m+3WHUQRERGsXLmScuXKMWHCBCpWrEjLli1ZvXo1169ft/rXdHuSBIQQGVKkiBl4NmeOuZPo4kUz//Ho0eDtDQsXQocOpl2tWmZw2vr1cP26Hz169GDHjh2cPHmSl156iRMnTvDss89SunRp+vfvz549e6S7yCKSBIQQmeLvb+Y/nj4ddu40ZwA7dsCkSVC0KLzxhjmTKFECqlY1E+Hs3Xsvzz8/iV9++YXt27fTsWNHVq9ezaOPPkpgYCCzZs3izJkzVv9qbsXL6gCEELlD3rzQpIl5gJny8uBBczE5LAzWrDFlKwAqVfKgceOmNG3alJEjX2P//nW8887bjB07lpCQEFq1akWbNm0oU6YMZcqUoXTp0pQuXRpfX1+rfr1cS5KAEMIpfHygfn3zGDMG4uPNtYWwMJMYNm+Gd94ByE/Zsn147LE+TJnyN3/+Gconn8xh69atKfZZuHBhSpcufSsxJE4QiZ/ny5fP1b9ujiVJQAjhEp6epqT1Qw/B8OGQkGDGJ9xMCtu3w5o1JYChFC8+hEaNYvHzu4a392VstgvY7X8THX2G69d/4+zZXwgPD+fs2fXExaW8uFywYMHbJosyZcrg5+fn+gORzTglCSilPIAlQE0gGuintT6RaP3zwAAgDpimtf5EKVUMWAPkBf4Eemut5dYBIXIpDw9Tu+iBB0yJa7sdTpy42X1k4/hxH06c8CEysjCXLlVMcz/+/gkUKhSLn98N8uS5gqfnRez2v4mNjeC3337nxx9PcvHit8TG/glEApcAu2Nb/9smi/Pnz3P16lXy5s2Lp6enKw6NSznrTKAj4Ku1rq+UqgfMBf4DoJQqBQwF6gC+wC6l1BfAS8AarfU7SqlgTJKY76T4hBDZjM0GVaqYR9++SdfFxsL58+Y21chIU+bi5vPISA/OnctDZGQeIiMLce5cOc6dSzmo7SYPjwTy54/G1/ca3t4XuXbtPMeORXDo0BmuXTtFbOx3mGQRCZwDrgExQCw+Pjb8/PLi5+dHvnz5UjzuZHlabb28XNtB46x3awhsBdBa71VK1Um0ri7wjdY6GohWSp0Aaji2meFo85njuSQBIQTe3lCqlHlkhN0O16+nljDg3DkPIiPzOh7FOHeuCpGR8M8/posqPTExEBsbzz//xOPhEYfNFovNFgvEAjHY7dHY7TEkJESTkBCF3R51a93NRGKSS+LXSZ97eibg42PHx8eDPHls5Mljo2BBD9au7UHVqgGZOXzpclYSKABcTvQ6XinlpbWOS2XdFaBgsuU3lwkhxB2z2cDPzzwqVMjYNgkJ5jbXxEnj3Dk4eTKCokVLERNzMwl4EhPjSUyMj+M1idYlfR4dbSc6OoGoqARiYuxER9sTtbMRG2sjLs5GfLwHCQnmjv34eLhxwzxu+v132L//T6pWzfpj5awk8A/gn+i1hyMBpLbOH9NJd3P5jUTLUggPD89UQFFRUZneNjeS45GUHI9/ybEwA96KFIGAAKhdOwpf34tOf8/4eIiLM0khNvbmA2JjbXh6wj33xBIefvn2O7pDzkoC3wAdgHWOawKJZ67+FpiulPIF8gCBwI+ObdoC7wBtgJ2p7TgwMDBTAYWHh2d629xIjkdScjz+JcciqdxwPA4ePJjmOmclgY1AS6XUbsAG9FZKjQBOaK03K6UWYT7kPYAQrXWUUmoasNJx51Ak0M1JsQkhhHBwShLQWicAA5MtPpZo/RvAG8m2OQs87ox4hBBCpE5qBwkhhBuTJCCEEG5MkoAQQrgxSQJCCOHGJAkIIYQbs+Wk2XwOHjyYc4IVQohspHbt2rbUlueoJCCEECJrSXeQEEK4MUkCQgjhxnL9zGK3m+DGnSilvIEVQEVM3aZpWuvNlgaVDSilSgAHgZZa62O3a5+bKaXGAk8APsASrfVbFodkCcf/lZWY/yvxwPO59W/DHc4Ebk1wAwRjJrhxV88C57XWjTBF+l61OB7LOf6zv46pXuvWlFJNgEeBBsBjQDlLA7JWW8BLa/0oMAWYbnE8TuMOSSDJBDeYGc3c1QfAhESv49Jq6EbmAMswU5q6u9aYir8bgY+BT6wNx1I/A16OnoQCmBlfciV3SAKpTnBjVTBW0lpf1VpfUUr5Ax8C462OyUpKqV7AOa31NqtjySaKYb4kdcYUgFytlEr1tkI3cBXTFXQMU+xykaXROJE7JIH0JrhxO0qpcsAO4F2t9Rqr47FYH0zJ86+AB4FVjjmw3dV5YJvWOkZrrYEooLjFMVllOOZYBGCuJ650zIGS67jDN+L0JrhxK0qpksDnwH+11v+zOh6raa0b33zuSAQDtdYR1kVkuV1AkFJqHlAa8MMkBnd0kX+7gC4A3oCndeE4jzskgRQT3Fgcj5XGAYWBCUqpm9cG2mit3f6iqACt9SdKqcaY2f88gBe01vEWh2WV+cAKpdROzJ1S47TW1yyOySlkxLAQQrgxd7gmIIQQIg2SBIQQwo1JEhBCCDcmSUAIIdyYJAEhhHBj7nCLqMgFHHVtPgKqa61/dyybBRzTWr+TyX1WBEK11vWyKMzE+/YEPsXca99Ba33Rsfwd4CHMvedeQCQwXGt96g73XwR4XGu9xrHPUK311qz7DYS7kDMBkZPEAG/nkFIGpYFiWuuGNxNAIqO11k201g0xBQ3XZWL/NTDVPoW4K3ImIHKS7TgGMZGoAmryb/RKqb1AV6AXcB+mJk4RTEnxTkAA0BOIAIorpTYDJYAtWuupjtIaywFfTOmE/pjRoh9jRtB+qrWenej9uwPDMKXKjzvaLweqKKVe11oPSOsX0lrvVErFKqXuc2yf5H211r8rpWZiavr4A+Fa695ACFBTKdXfsasBSqnRQEFgEPADJrkUBPJiEs9XGTvMwp3ImYDIaQYBw5VSVTLY/obW+nFgA9BWa90BmIVJEgD5gecw5ZPbKKVqYiqLLtJaN3U8n+VoWwpolSwBFAUmA80c3+wvAQOAwcDR9BJAImcxiSrF+yqlCgAXtdYtMWWe6ymlymJKG2/XWi937OOg1roZsBiT/Co74u0AdAPyZfB4CTcjSUDkKFrr85hv3e+Q9t9v4u6i7xw/LwFHHc8vYr5tAxzWWl92lEf4FnOWUB0Y56gn9BLmLAHglNY6Jtl7VQJ+0lpfcbz+GnjgDn+tCsAfabzvDaCEUup9zLwH+TF1bJI76PgZAeTTWv8EvAa8jzkDkv/rIlXyhyFyHK31x4DGfOMF03VSQinlqZQqBNybqPnt6qIEKqXyO8qLPwL8hCkfPEZr3QTzrf5DR9uEVLY/BVRVSvk5Xj+GqUWfIUqplsB1rfUfabxvG6Cc1vr/MLWf8mKSXAJJ///ak+23OuCvtW6H6fpanNGYhHuRawIipxoGNAfQWkcopb4A9gMnHI+MugCsxZRMXqu1PqqUGgUsdZQOzgsEpbWx1jpSKTUR2KGUSnC8dzCmKyYts5VSwZhpC68AXRzLU3vfU5iCf3sx1wxOAmWAX4DqSqlhabzHcWCiUqoH5oL6S7c5DsJNSQE5IYRwY9IdJIQQbkySgBBCuDFJAkII4cYkCQghhBuTJCCEEG5MkoAQQrgxSQJCCOHGJAkIIYQb+3+j7KAjZ1ggmQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "lamda = np.sum(freq*np.arange(freq.size)) / np.sum(freq)\n",
    "pmf = []\n",
    "for i in range(freq.size):\n",
    "    pmf.append(np.exp(-lamda) * (lamda**i) / mt.factorial(i))\n",
    "\n",
    "plt.plot(np.arange(freq.size), freq/np.sum(freq), color='black')\n",
    "plt.plot(np.arange(freq.size), pmf, color='blue')\n",
    "plt.xlabel(\"Number of Deaths\")\n",
    "plt.ylabel(\"Proportion\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Problem 2\n",
    "\n",
    "Fit the $\\alpha$, $\\mu_1$ and $\\mu_2$ parameters for the mixture of two Poissons model proposed in Question 11. You must write your own fitting code (cannot use `scipy optimize` or R `optim`). You can write any sort of optimizer you like such as gradient ascent, Newton's method (or any of its varients), or an EM approach. Discuss how you select starting parameters and monitor convergence. You can confirm your results by also using Python or R built in optimizers."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The answer to problem 1 used standard notation for the data i.e. $x_1$ was the number of deaths on day 1. Now we let $n_i = $ number of days with $i$ deaths. So\n",
    "$$n_i = \\sum_j \\mathbb{1}_{X_j=i}$$\n",
    "\n",
    "The likelihood is\n",
    "$$f(\\vec{n}|\\alpha,\\mu_1,\\mu_2) = \\prod_{i=0}^9 \\left(\\underbrace{\\alpha e^{-\\mu_1} \\frac{\\mu_1^i}{i!} + (1-\\alpha)e^{-\\mu_2}\\frac{\\mu_2^i}{i!}}_{\\equiv g_i}\\right)^{n_i}$$\n",
    "\n",
    "So the log likelihood is\n",
    "$$ \\log \\, f = \\sum_{i=0}^9 n_i \\log(g_i) $$\n",
    "\n",
    "The gradient is\n",
    "$$ \\frac{\\partial log f}{\\partial \\alpha} = \\sum_{i=0}^9 \\frac{n_i}{g_i} \\frac{\\partial g_i}{\\partial \\alpha}$$\n",
    "$$ \\frac{\\partial log f}{\\partial \\mu_j} = \\sum_{i=0}^9 \\frac{n_i}{g_i} \\frac{\\partial g_i}{\\partial \\mu_j}$$\n",
    "for $j=1,2$ and where\n",
    "$$\\frac{\\partial g_i}{\\partial \\alpha} = e^{-\\mu_1} \\frac{\\mu_1^i}{i!} - e^{-\\mu_2}\\frac{\\mu_2^i}{i!}$$\n",
    "$$\\frac{\\partial g_i}{\\partial \\mu_1} = -\\alpha e^{-\\mu_1} \\frac{\\mu_1^i}{i!} + \\alpha e^{-\\mu_1} \\frac{\\mu_1^{i-1}\\mathbb{1}_{i > 0}}{(i-1)!}$$\n",
    "$$\\frac{\\partial g_i}{\\partial \\mu_1} = -(1-\\alpha) e^{-\\mu_2} \\frac{\\mu_2^i}{i!} + (1-\\alpha) e^{-\\mu_2} \\frac{\\mu_2^{i-1}\\mathbb{1}_{i > 0}}{(i-1)!}.$$\n",
    "Note the indicator functions. These are important for the second terms in the $\\partial g_i / \\partial \\mu_j$ which are $0$ for $i=0$ (because we are differentiating a constant)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now write functions to implement the log likelihood and gradient."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "import scipy\n",
    "\n",
    "def g(alpha,mu1,mu2,N):\n",
    "    iis = np.arange(N)\n",
    "    iisf = scipy.special.factorial(iis)\n",
    "    return alpha*np.exp(-mu1)*np.power(mu1,iis)/iisf + (1-alpha)*np.exp(-mu2)*np.power(mu2,iis)/iisf\n",
    "\n",
    "## params = alpha, mu1, mu2 1-d array of parameters\n",
    "## n = 1-d array of counts of the number of deaths\n",
    "def log_likelihood(params,n):\n",
    "    return np.sum(n*np.log(g(params[0],params[1],params[2],n.size)))\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The gradient is a bit more complicated."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def dgda(mu1,mu2,N):\n",
    "    iis = np.arange(N)\n",
    "    iisf = scipy.special.factorial(iis)\n",
    "    return np.exp(-mu1)*np.power(mu1,iis)/iisf - np.exp(-mu2)*np.power(mu2,iis)/iisf\n",
    "\n",
    "## omits the leading alpha for mu1 and leading (1-alpha) for mu2\n",
    "def dgdmu(mu,N):\n",
    "    iis = np.arange(N)\n",
    "    iisf = scipy.special.factorial(iis)\n",
    "    temp = np.exp(-mu)*np.power(mu,iis[:(N-1)])/iisf[:(N-1)]\n",
    "    return -np.exp(-mu)*np.power(mu,iis)/iisf + np.concatenate((np.zeros(1),temp))\n",
    "\n",
    "## params = alpha, mu1, mu2 1-d array of parameters\n",
    "## n = 1-d array of counts of the number of deaths\n",
    "def gradient_log_likelihood(params,n):\n",
    "    alpha = params[0]\n",
    "    mu1 = params[1]\n",
    "    mu2 = params[2]\n",
    "    N = n.size\n",
    "    nog = n / g(alpha,mu1,mu2,N)\n",
    "    c1 = dgda(mu1,mu2,N)\n",
    "    c2 = alpha*dgdmu(mu1,N)\n",
    "    c3 = (1-alpha)*dgdmu(mu2,N)\n",
    "    return np.array((np.sum(nog*c1),np.sum(nog*c2),np.sum(nog*c3)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now implement gradient ascent with step halving."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0.3,  1. ,  2.5])"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "## gradient ascent with step halving\n",
    "Niter = 1000\n",
    "param_initial = np.array((0.3,1.0,2.5))\n",
    "param = np.zeros((Niter+1)*3).reshape(Niter+1,3)\n",
    "param[0,:] = param_initial\n",
    "param[0,:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\Kanchan\\Anaconda3\\lib\\site-packages\\ipykernel_launcher.py:11: RuntimeWarning: invalid value encountered in log\n",
      "  # This is added back by InteractiveShellApp.init_path()\n"
     ]
    }
   ],
   "source": [
    "for ii in np.arange(Niter):\n",
    "    grad = gradient_log_likelihood(param[ii,:],freq)\n",
    "    half = True\n",
    "    while half:\n",
    "        param_temp = param[ii,:] + grad\n",
    "        ## the new likelihood could produce nan because\n",
    "        ## we take the log of a negative number\n",
    "        ## if so, step halve\n",
    "        newlik = log_likelihood(param_temp,freq)\n",
    "        if np.isnan(newlik) or (newlik < log_likelihood(param[ii,:],freq)):\n",
    "            grad = grad / 2\n",
    "        else:\n",
    "            half = False\n",
    "    param[ii+1,:] = param_temp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.3       ,  1.        ,  2.5       ],\n",
       "       [ 0.22706808,  1.01877166,  2.52908317],\n",
       "       [ 0.25187606,  1.01496951,  2.52085513],\n",
       "       ..., \n",
       "       [ 0.3598852 ,  1.25609473,  2.6634041 ],\n",
       "       [ 0.3598852 ,  1.25609473,  2.6634041 ],\n",
       "       [ 0.3598852 ,  1.25609473,  2.6634041 ]])"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "## looks correct\n",
    "param"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Functionalize Gradient Code\n",
    "\n",
    "We now turn gradient code into a function so we can easily plot parameter path from several starting points."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## arguments\n",
    "##        loglike : log likelihood function\n",
    "##       dloglike : gradient (derivative of log likelihood)\n",
    "##  param_initial : initial parameter values for algorithm\n",
    "##              X : data\n",
    "##          Niter : number of iterations\n",
    "##\n",
    "##\n",
    "##    value\n",
    "##       param    :  2-D array with rows giving the parameter values at\n",
    "##                   each iteration\n",
    "##\n",
    "##    note: the function loglike and dloglike must have arguments\n",
    "##          (parameters,data)\n",
    "##          loglike must return a scalar and dloglike must return a 1-D\n",
    "##          np array of length p\n",
    "def gradient_ascent(loglike,dloglike,param_initial,X,Niter):\n",
    "    p = param_initial.size\n",
    "    param = np.zeros((Niter+1)*p).reshape(Niter+1,p)\n",
    "    param[0,:] = param_initial\n",
    "    for ii in np.arange(Niter):\n",
    "        grad = dloglike(param[ii,:],X)\n",
    "        half = True\n",
    "        while half:\n",
    "            param_temp = param[ii,:] + grad\n",
    "            ## the new likelihood could produce nan because\n",
    "            ## we take the log of a negative number\n",
    "            ## if so, step halve\n",
    "            newlik = loglike(param_temp,X)\n",
    "            if np.isnan(newlik) or (newlik < log_likelihood(param[ii,:],X)):\n",
    "                grad = grad / 2\n",
    "            else:\n",
    "                half = False\n",
    "        param[ii+1,:] = param_temp\n",
    "    return param"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We call function with 3 different starting values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\Kanchan\\Anaconda3\\lib\\site-packages\\ipykernel_launcher.py:11: RuntimeWarning: invalid value encountered in log\n",
      "  # This is added back by InteractiveShellApp.init_path()\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array([[ 0.3       ,  1.        ,  2.5       ],\n",
       "       [ 0.22706808,  1.01877166,  2.52908317],\n",
       "       [ 0.25187606,  1.01496951,  2.52085513],\n",
       "       ..., \n",
       "       [ 0.3598852 ,  1.25609473,  2.6634041 ],\n",
       "       [ 0.3598852 ,  1.25609473,  2.6634041 ],\n",
       "       [ 0.3598852 ,  1.25609473,  2.6634041 ]])"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "## call function\n",
    "Niter = 1000\n",
    "\n",
    "param_initial = np.array((0.3,1.0,2.5))\n",
    "param1 = gradient_ascent(log_likelihood,gradient_log_likelihood,param_initial,freq,Niter)\n",
    "param1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\Kanchan\\Anaconda3\\lib\\site-packages\\ipykernel_launcher.py:11: RuntimeWarning: invalid value encountered in log\n",
      "  # This is added back by InteractiveShellApp.init_path()\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array([[ 0.5       ,  0.5       ,  4.1       ],\n",
       "       [ 0.24923513,  0.90173495,  3.95861443],\n",
       "       [ 0.56916043,  0.90473252,  3.86359071],\n",
       "       ..., \n",
       "       [ 0.3598856 ,  1.25609547,  2.66340461],\n",
       "       [ 0.3598856 ,  1.25609547,  2.66340461],\n",
       "       [ 0.3598856 ,  1.25609547,  2.66340461]])"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "param_initial = np.array((0.5,0.5,4.1))\n",
    "param2 = gradient_ascent(log_likelihood,gradient_log_likelihood,param_initial,freq,Niter)\n",
    "param2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\Kanchan\\Anaconda3\\lib\\site-packages\\ipykernel_launcher.py:11: RuntimeWarning: divide by zero encountered in log\n",
      "  # This is added back by InteractiveShellApp.init_path()\n",
      "C:\\Users\\Kanchan\\Anaconda3\\lib\\site-packages\\ipykernel_launcher.py:11: RuntimeWarning: invalid value encountered in log\n",
      "  # This is added back by InteractiveShellApp.init_path()\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array([[ 0.8       ,  0.1       ,  1.        ],\n",
       "       [-0.07149436,  0.48151832,  1.33857954],\n",
       "       [-0.22943851,  0.47434128,  1.65158431],\n",
       "       ..., \n",
       "       [ 0.35988535,  1.25609502,  2.6634043 ],\n",
       "       [ 0.35988535,  1.25609502,  2.6634043 ],\n",
       "       [ 0.35988535,  1.25609502,  2.6634043 ]])"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "param_initial = np.array((0.8,0.1,1.0))\n",
    "param3 = gradient_ascent(log_likelihood,gradient_log_likelihood,param_initial,freq,Niter)\n",
    "param3"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plot Parameter Paths\n",
    "\n",
    "We plot the $\\mu_1$, $\\mu_2$ parameter paths. You could also plot these paths for $\\alpha$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "plt.style.use('ggplot')\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAECCAYAAAD9z2x7AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAIABJREFUeJzt3Xl8VOXZ//HPLElmJnsYhhA2CQgo\nSTGAEQJSwYBFrAiiVMsPbbEqCBKBQpCyG+F5QNAHH8S6tETQUiEo8EgtiICgFReCoIDITkggIQtZ\nZsksvz+GDAlZSGC2JNf79eJF5szkzDdBr3PmPve5boXD4XAghBCiyVP6OoAQQgjvkIIvhBDNhBR8\nIYRoJqTgCyFEMyEFXwghmgkp+EII0UxIwRdCiGZC7esA13P+/HlfR0Cv15OXl+frGFX4YyaQXA3l\nj7n8MRNIrvqKiYmp9Tk5wxdCiGZCCr4QQjQTUvCFEKKZkIIvhBDNhBR8IYRoJqTgN2MZGVoSEw20\nbduaxEQDGRlaX0cSQniQ30/LFJ6RkaFl+vRwjEbnMT8rS8306eEAjBxp9GU0IYSHSMFvoMzMTObP\nn0+HDh1QKBSYzWaSk5MZOXJkrd+zefNmhg4dyqFDh9i0aRNz5syp8z1MJhPTpk1j+vTptG/f3t0/\nAgCLF4e6in0Fo1HJjBnh/PKLmjZtbLRta6NNG+cfrVaWTRCisZOCfwMSEhJcRdtisTB27FiGDBlC\nSEhIja9fu3Yt9913X732ffToUZYtW0Zubq7b8tbk/HlVjdvLyhSsWBGC3a6osj0qquoBoOJPxbYW\nLewoFDXuUgjhJxp1wdfmfIgue51b91nWejTG6Efq/Xqj0YhKpUKlUpGZmUl6ejrgPEufOXMmP/zw\nA/n5+SxYsIBRo0Zx7tw5ZsyYQUFBAUlJSTz55JNV9mexWFi4cCEvv/yyO3+samJibGRlVf/nb9PG\nxpdfXuTCBRXnzqnIyrr6d1aWiuPH1ezaFURZWdVPBxqNg5gYGx07KjEYwqsdFFq3thEY6NEfSQhx\nHY264PvK/v37SUlJQalUolKpmDRpElqtllOnTvHiiy+i1+tZs2YNu3btYsyYMbz33nvMmTOHn376\nyVXQ7XY7o0ePrlbw4+PjvfIzpKYWVxnDB9Bq7aSmFqNW4yrWNXE4oLBQceUgoK5yULhwIYiDBzVc\nvFj1E4RC4aBVKzsxMVc/FbRtayUm5upBISzMIZ8ShPAgjxb8oqIiUlNT+ctf/kKbNm1c27/99ls2\nbNiAUqlk4MCBJCcn39D+jdGPNOhs3F0qD+lUptfrWbFiBVqtlry8POLi4qq9pmPHjgReOdVVqWoe\nVvGGiguzixeHcv68ipgYG6mpxfW6YKtQQGSkg8hIK3Fx1irPVfQVMZkgO/vqJwPnQcF5cPjhhwD+\n9S8NFkvV6h4SYqdtW9s1BwXn3zExVqKj7fjwVyZEo+exgm+1WvnrX//qKm6Vt69evZpFixah0WiY\nPXs2vXv3JiIiwlNRvGbp0qW8//776HQ6Fi1aRMX68EqlErvdDoDCj05hR440emxGjkYDHTva6Nix\n5k8Jdjvk5SkrHQxUnD9/9ZPC998HUlhYddhIpXLQurWt1oNCmzY2dDq5uCxEbTxW8N977z0GDx7M\nRx99VGV7VlYW0dHRrgucXbt25fDhw/Tt29dTUbxm8ODBTJgwgZCQEKKiorh06RLgHKZJTU3liSee\n8HFC/6FUgsFgx2Cwk5BQXuNrSksV13xCuPr1vn2BfPyxCput6gE0MrL6QaBNGxtxcQqCg5Xo9XJx\nWTRfCkfFaagb7dy5k0uXLvHwww8zb948/vSnP7mGdI4cOcLWrVt54YUXAFi3bh16vZ577723xn1Z\nLBZ3x2swtVqN1Wq9/gu9yB8zgXdzWa2QnQ1nzig4e1bB2bNw+rTz6zNnnNtLSqpW96AgB+3aQbt2\nDjp0cNCunYP27bnyt4O2bSEoyCvxAf/8d/THTCC56uvaUZXKPHKG//nnnwNw8OBBTp06xeuvv86M\nGTOIiIhAq9ViMplcrzUajQQHB9e6L3/oM+1v/a7BPzOB93NptdC1q/PPtRwOKCpyfkq4fDmSI0dK\nycpSuz4pbN2q4sKF6v8LGAzVp51WXGBu29ZGeLj7Li7747+jP2YCyVVfdfXD90jBnz9/vuvrijP8\nijH6Nm3akJ2dTUlJCRqNhsOHD/Pggw96IoZo5hQKiIhwEBFhRa930LdvWbXXmM1XLy5XvY6g5qef\nAti+XYPJVLW6Dx5s4u9/z/fWjyGE23htWuaePXswmUwkJyczduxY0tLSsNvtDBw4kKioKG/FEKKK\noCC45RYbt9xS+xTUS5eUrk8Ff/tbMN99F+DllEK4h8cL/rx58wCqTMvs3bs3vXv39vRbC3HTFArQ\n6+3o9XbuuKOco0fVfPVVEOXlECB1XzQy0i1TiAZo2dI5vTYvT/7XEY2P/FcrRAMYDBUFX+4AE41P\nk26tEHJmJZbQHlgi+7m2BRbsJbD4ACXtJ9zQPnNychg3bhy33nqra1tCQgL9+vVj7969DZ5rn5KS\nwpQpU2rsivm///u/tGvXTi5q+xG93jnWf/GinCuJxqdJF3xLaA8if3qWgttXYYnsR2DBXtfjm9Gh\nQwdeffXVats7d+58U/utUFhYyKJFizh37hyjR492yz6Fe1Sc4c+cGU7XrlbatrXRrp1z2ma7ds6v\npXOo8FdNu+BH9qPg9lVE/vQsZTFj0Z1PdxV/d8vMzHT1uh8zZgxxcXGcPXuWyMhI5s+fj8lkYsmS\nJZSWllJUVMSwYcMYPnx4jfsyGo088cQT7Nu3z+05xc1p29bG9OmXOXQogLNna24BodHYXcVfDgjC\nnzTpgg/Ool8WM5bQ069S3CHFLcX+9OnTpKSkuB7PmjWryvPZ2dksW7YMg8HAxIkTOXr0KGq1mkGD\nBjFgwADy8vJISUmpteC3bt2a1q1bS8H3QwoFTJ5cUmVbcbGCc+dUnD3rbBDn/FtVrwNC584qWrYM\nkQOC8IomX/ADC/aiO59OcYcUdOfTMUck3XTRr2lIJysry/V1eHg4BoMBAIPBgMViwWAwsH79er74\n4gt0Op1f3Yotbk5oqIPbbrNy221WwFzt+boOCJmZSvLzw6q8Xqu1V/l00K6dtdLXNqKi5IAgbkyT\nLviVx+wtkf0wRyRVeexN69ato3v37gwfPpz9+/fz9ddfe/X9he/UdUDQ6/WcOHGJc+cqDgLX/4RQ\n0wEhMdFCr141N6ETokLTLvjFB6oU94ox/cDiA14v+ElJSSxfvpzt27cTFhaGSqXyi8ZwwvfCwhzc\nfruV22+v+RPC5cuKeh0QRo4s4y9/uUyrVnYv/wSisfBIt0x3On/+vK8j+F1zJPDPTCC5GsoduQoL\nFbz1VggrV4YQGOhg2rRi/vCHUtQ3eDrXlH9XnuBvuepqniaTiYVo5CIiHPz5z8Xs2HGRxEQL8+aF\n85vftOTrr2URYVGVFHwhmoiOHW2kp+fzzjv5XL6sYORIPc8/HyE3iQkX+S9BiCZEoYDf/MbErl25\nPP98MZs3axkwwMA77wQjE8NEk75oK0RzpdU6mDGjmFGjypg9O5w5c8L5xz903HefiagoO5GRdtff\nFV/rdO5b2EX4Jyn4QjRhnTrZWLs2n61bNaSlhfHqqyE4HDVX9cBAB1FRdvR6BWFhLaocDGr6OzLS\nTliYHCQaEyn4QjRxCgXcf7+J++83YbNBUZGS/HwFBQUqCgoU5OcrKShQuv4uKdFy4QL8/LOa/Hwl\nhYXKaovFV1CrHUREVP+00LNnOb/7XZkcDPxMky/42gsZhJ5YjMp8HltQDMWxqRhbjbzh/eXk5LBg\nwQJWrlxZ4/P5+fmsXr2aF154gQMHDhASEkKnTp2YM2cOCxYsaPA+CwsLmThxIu+++26dixMLUR8q\nFURF2XEuMlfzKl96fQB5eZdcj+12570AlQ8Klb+uvO3UKTXffqvk/feDKSpS8Oyzpd75wUS9NOmC\nr72QQfjR6SjtRgDU5izCj04HuKmiX5eoqCheeOEFALZu3cqgQYPo1KlTrcW+Lvv27eOtt96ioKDA\n3TGFqDelsmJtYBsdO9Z8kKjMbocJEyJZuDCcVq3sjBhh9EJKUR9NuuCHnljsKvYVlHYjoScWu6Xg\np6Sk0LlzZ06ePElZWRlz584FYMGCBUyePJl9+/Zx7NgxOnTowPjx48nIyCAzM5P09HQATCYTM2fO\nJKCWtfKUSiVLly7lmWeeuemsQniLUgmvvVZAXp6SF16IQK+3cffdcle5P2jS0zJV5prv0q1t+43o\n1q0br7zyCr169WLHjh2u7V27diUxMZFnnnmGVq1aubafOnWKF198kWXLlpGUlMSuXbtq3Xfv3r0J\nDw93W1YhvCUoCN55J59Onaw89VQUa9bo2LYtiG++CeCXX1RcuqSUaaI+0KTP8G1BMajNWTVud5eK\nla8MBgP5+fnXfb1er2fFihVotVry8vKIi4tzWxYh/El4uIP33rvEyJF6ZsyIqPE1YWHOC70REdf+\n7ai2vUcP5zUIceOadMEvjk2tMoYPYFdqKY5Nddt7KOqYhqBQKLDbqzayWrp0Ke+//z46nY5Fixbh\n562MhLgpMTF2du26SHa2ioIC54yfwkLlla8Vrm0Vf586paawUElRkaLa9FG12sGcOcH88Y+lMvvn\nBjXpgl8xTu/OWToNcdttt/HWW2/RunVr17bBgwczYcIEQkJCiIqK4tKlS3XsQYjGLygIbrnFxi23\nXP+CbwXn9FFFlQPEP/8ZyZw54ezfH8B//3cROp2cLDWUdMusB3/rhgf+mQkkV0P5Yy5/zAQQFaVn\n7lwTS5aE0q2blcceKyM62ub6YzDYqWX+g0f52++rrm6ZHjvDt9vtrFq1iuzsbJRKJePHjyc6Otr1\n/JYtW9ixYwdhYc7Vfp5++uk6gwohmjel0rm8ZI8e5UyeHMGcOVUnNCgUDlq2tFc6CFz9unXrq1+H\nhjbfu4M9VvC//fZbABYuXMiPP/5Ieno606dPdz1/4sQJJk6cSGxsrKciCCGaoHvuMbN//wUKCpRk\nZyvJyVFV+uN8fPasmn37VNVWCwPQ6ezXHAxstGplp317K/fea0bZhOcueqzgJyYm0qtXLwByc3Or\nTS88efIkGzdupLCwkJ49ezJixAhPRRFCNDFKJbRoYadFCztxcbXP7zQa4cKFqgeE7Oyrj/ftC+TC\nBRXl5c5T/nXr8ujfv+neM+DxMfzXX3+db775hilTptCjRw/X9g8//JD77rsPnU7HkiVLGDJkiOsA\nUZk/LAOoVqv9btFxf8wEkquh/DGXP2YCz+U6cQJ69w4gPt7B9u3WBl8H8LffV10tWLxy0bawsNB1\ns5FGo8HhcGA0GtHpdAB8+umnFBcXM2rUqGrfKxdta+aPmUByNZQ/5vLHTOC5XC+8EME//6lDq7Wj\n11+d+19xP0CLFnaGDy+jc+faeg/51+/LJxdtd+/ezaVLlxgxYgSBgYEoFAqUVwbHjEYjU6dOZfny\n5QQFBXHo0CEGDhzoqShCCFGrp58uoX17q6sBXMU00DNnrt4TsGJFCE89VcrkycWEhvr1xMY6eewM\n32QysXLlSoqKirBarTz00EOYzWZMJhPJycns3r2brVu3olariY+P59FHH61xP/52hp+Zmcn8+fPp\n0KEDCoUCs9lMcnIyI0fWPrd/8+bNDB06lEOHDrFp0ybmzJlT62s/++wz1q9fj0qlIjY2lpSUFNeB\nsrZM/kRyNYw/5vLHTOC7XLm5ShYvDuUf/wjGYLCRmnqZhIRytFoHOp2Dtm2jKC3N85uLvT45w9do\nNEyZMqXW5wcMGMCAAQM89fYelZCQ4CraFouFsWPHMmTIEEJCQmp8/dq1a7nvvvuuu1+z2cy7777L\nO++8g0ajYeHChXz11Vf069fPrfmFEPXXsqWdV14pYswY5+phU6ZE1vCqGDQaO6NGGVm0qMhviv+1\nGvWdth9+qGXdOp1b9zl6dBmPPFL/dq5GoxGVSoVKpaqxE+YPP/xAfn4+CxYsYNSoUZw7d44ZM2ZQ\nUFBAUlISTz75pGtfAQEBrFixAo1GA4DNZpMe+EL4iYSEcjZtyuOrrwK5dEmJ0ajAaFQAIVy6ZOTk\nSRVr1gQTGupg9OgyAgMdBAU5CAx0EBzsICjI1z9BIy/4vrJ//37XUItKpWLSpElotVpXJ0y9Xs+a\nNWvYtWsXY8aM4b333mPOnDn89NNPWCwWFi5ciN1uZ/To0VUKvlKpJMq5MgUZGRkYjUZ69+7to59S\nCHEtpRL69as6c1Cv15GXV4zD4VxL+I03Qnjjjaqf9rVaO9u35zaovYQnNOqC/8gjxgadjbtL5SGd\nyurTCbNjx46us3ZVDa3/7HY7b775JufOnWP+/Pl1NmcTQvgPhQIWLy7igQdMFBYqMJudf4xGBS+/\nHMZTT0XRtatz7D883MHkycWEhXn3AnCjLvj+prZOmEql0tU183oFfNmyZQQEBLBw4cIaL9YKIfyX\nSgW//rW52najUcHmzVoyMwMxmRTk5Kg4c0bFXXdZ6NXLQkJCuVfyScF3o9o6YcbHx5OamsoTTzxR\n5/f//PPPfPLJJ8THx7sueD/88MPcfffdHs8uhPCcyZNLmDy5xPV43LhItm7V8MknWlQqB1OnFvPb\n3xqJjfXskI90y6wHf5ym5o+ZQHI1lD/m8sdM0PRy2Wxw5oyKP/whimPHAoiOtrF1ay4tW9pvqrlb\nXdMyZcxACCF8QKWCjh1t7NyZy6ZNuVy4oCQhIZqEhFZXZv+4nxR8IYTwsV69ysnIcA4B5+aqePzx\nKI+8jxR8IYTwA4mJFt56y7ku9jffeOb+Gyn4QgjhJzIytAAMHWryyP6l4AshhB84eDCArVudBf+V\nVwo98h5S8IUQwse+/jqQmTOdi0R9/vlFj92Q1aTn4a9cGUKPHpYqt0Lv3RvIgQOBTJhQUsd31i4n\nJ4dx48Zx6623urYlJCTQr18/9u7de9259tdKSUlhypQptG/f3rXtl19+4X/+539QKpUEBAQwc+ZM\nV8sFIUTT8q9/aRg3zvn/9+zZRXTp4rnFVJp0we/Rw8Kzz0ayalUB/fpZ2Ls30PX4ZnTo0IFXX321\n2vbOnTvf1H4rrFixgueff57OnTuzadMmPvjgA5577jm37FsI4R9sNli3TkdaWhjdupXz8stF3HWX\nZ1f4a9IFv18/C6tWFfDss5GMHVtGerrOVfzdLTMz09XrfsyYMcTFxXH27FkiIyOZP38+JpOJJUuW\nUFpaSlFREcOGDWP48OE17mvOnDm0aNECkI6ZQjRFGRlaJk1ytllu29bKqlUF3Hqr55dJbNIFH5xF\nf+zYMl59NZSUlGK3FPvTp0+TkpLiejxr1qwqz2dnZ7Ns2TIMBgMTJ07k6NGjqNVqBg0axIABA8jL\nyyMlJaXWgl9R7A8dOsRHH33Ea6+9dtOZhRC+l5+vYNmyUP72N2c3zRdeKGbq1OKburO2IZp8wd+7\nN5D0dB0pKcWkp+tISjLfdNGvaUgnKyvL9XV4eDgGgwEAg8GAxWLBYDCwfv16vvjiC3Q63XUXPd6x\nYwdr165l0aJFRERE3FReIYTvnTun4q67WgHQrVu5187qK2vSBb/ymH2/fhaSksxVHnvTunXr6N69\nO8OHD2f//v18/fXXtb5227ZtbN68meXLlxMWFubFlEIId3M4YM6cMN5992qP/H//O5cauqN7XJMu\n+AcOBFYp7hVj+gcOBHq94CclJbF8+XK2b99OWFgYKpUKi6V6BpvNxooVKzAYDK6e+z169OAPf/iD\nV/MKIW7eiRMq7r67levxxInFzJhR7LMlEKVbZj34Y5c+f8wEkquh/DGXP2aCxpXL4YAxY6LYudO5\nXGmnTuV89lkuAQGez+OTRcyFEKI5+u67AMaObUFhofM0/pNPcunRwzsLnFyP3GkrhBBuYDTCm28G\n8+CDLSksVPL888UcPpztN8Ue5AxfCCFuSnk5pKaGs2mTlqIiJXfcYWH69OIalzr0NSn4Qghxg374\nIYChQwMB582R69fn0bevdyeENITHCr7dbmfVqlVkZ2ejVCoZP3480dHRrue//fZbNmzYgFKpZODA\ngSQnJ3sqihBCuJXZDLGxVS+OfvddDtHRdh8lqh+PjeF/++23ACxcuJBHH32U9PR013NWq5XVq1cz\na9Ys5s+fz2effUZhoWfagQohhLuUlip4881gkpKuTrUcN87GuXPn/b7YgwcLfmJiIs888wwAubm5\nhIeHu57LysoiOjqakJAQ1Go1Xbt25fDhwx7JkZGhJTHRQNu2rUlMNLgWGLhROTk5TJgwodbn8/Pz\nWb58OQAHDhzg+PHjAK459Q3Z54cffsj48eMZP348q1evvqncQogb53DAe+/p6NKlNQsWhNOpk5UP\nPsjj3LnzrFxp81prhJvl0TF8lUrF66+/zjfffMOUKVNc241GIzqdzvVYq9VSVlZW4z70ev0Nv/8H\nHyiZMUNFWZnzXyMrS82MGRGEhoby2GP1Pxqr1WpXDrPZTEBAQK259Ho9aWlpgLM9wtChQ9Hr9axc\nubLW/de0z7Nnz7Jz507+8Y9/oFAoeOKJJ3jggQfo2rVrtUz+RHI1jD/m8sdM4JtcNhukpKhIT1di\nMjnryLPP2njtNQUQ5rNcN8rjF20nTpxIYWEhL774IsuWLUOj0aDVajGZri7hZTQaCQ4OrvH7b+ZG\ni1mzDK5iX6GsTMGsWTB4cP33W/nGioKCAsrLy10N0Dp37szJkycpKytj7ty5ACxYsIDJkyfzxRdf\ncOjQISIjIxk/fjwZGRlkZma6hrdMJhMzZ84kICDAtc8KAQEBvPzyyxQUFLheW1ZW5npNY7oJxR9I\nrvrzx0zg/VyV+9QDLF1ayMiRZQQFQeUY/vb7quvGK48N6ezevZuNGzcCEBgYiEKhQHnlfuI2bdqQ\nnZ1NSUkJVquVw4cP06VLF7dnOH++5mYVtW2/Ed26deOVV16hV69e7Nixw7W9a9eurmGtVq2ujved\nOnXKdfBLSkpi165dNe5XrVYTHh6Ow+HgjTfeoHPnzrRr185tuYUQ1VmtsGWLhuHD9VWK/T//mcdj\njzmLfWPmsTP8xMREVq5cydy5c7FarTz55JPs27cPk8lEcnIyY8eOJS0tDbvdzsCBAz2yolNMjI2s\nrOo/YkyMzW3vUbHylcFgID8//7qv1+v1rFixAq1WS15eHnFxcbW+1mKx8F//9V/odLoq7ZiFEO5V\nXKzggw90vPNOMOfOqenQwcqCBUWMHl1GSIhfd59pEI8VfI1GU2Xc/lq9e/emd+/ennp7AFJTi5k+\nPRyj8eoHGa3WTmpqsdveQ1HH1RqFQoHdXvVawdKlS3n//ffR6XQsWrSI2loZORwOZs2aRc+ePXns\nscfcllcIcZXRCEOGGDhxwlkK77rLzPz5lxk82OSTbpae1qRvvBo50gjA4sWhnD+vIibGRmpqsWu7\np91222289dZbtG7d2rVt8ODBTJgwgZCQEKKiorh06VKN37tnzx4OHDhAeXm5q5Xyn/70J7p37+6V\n7EI0ZZV701dISyvkySdrnjzSVEi3zHrwt4sy4J+ZQHI1lD/m8sdM4J5cixaF8vrroVW2rV+fR58+\nlhueWulvvy/plimEaLby8pT06BFdbfsHH1xiwAD/63fjSVLwhRBN0pQpEaxbp6u2ffv2i9x2m3eX\nFvQXUvCFEE3GL7+o2LxZy9//Hkxe3tWrrrGxVj79NBedzq9HsD1OCr4QolE7c0bFpk1aNm3S8uOP\nASgUDhITLfTubWHixBISEvynH72vScEXQjQ6WVlKNm/WsnmzlsxMZ2vinj0tzJtXxAMPGGnd2v8b\nmfmCFHwhRKNgMsHSpaEsX351lk18vIVZsy7z298aadfOfTdUNlVS8BsoMzOT+fPn06FDBxQKBWaz\nmeTkZEaOHFnr92zevJmhQ4dy6NAhNm3aVGfnzF27dvHBBx+gUCh44IEHGDZsmCd+DCEaBaMR3ngj\nhF27NBw/HkBBgfNsvl8/M4sXFxIbK0W+IaTg34CEhARX0bZYLIwdO5YhQ4YQEhJS4+vXrl3Lfffd\nd9392mw23nrrLVatWoVWq+UPf/gD/fv3r9JaWoim7p13gsnI0KJUwsGDAZSXOyfIP/qojREjCujf\n34JSVuO+IY264Gs//BDdunVu3WfZ6NEYH3mk3q83Go2oVCpUKlWNnTB/+OEH8vPzWbBgAaNGjeLc\nuXPMmDGDgoICkpKSePLJJ137UqlUrF69GpVKRUFBAQ6HA6325vr3C+HvzGbIylKxYYOOPXuC+PZb\n51l8YqKZZ54pITbWyh13lNOvXwR5ef67fGBj0KgLvq/s37+flJQUlEolKpWKSZMmodVqXZ0w9Xo9\na9asYdeuXYwZM4b33nuPOXPm8NNPP2GxWFi4cCF2u53Ro0dXKfjgLPq7d+/mtddeo0+fPqiaYkMP\n0eydPKmif/9WNT43e3YRjzxipEULufDqbo264BsfeaRBZ+PuUnlIp7L6dMLs2LEjgYHOM5jaivmA\nAQPo378/ixcv5t///jdDhw517w8ghJeVlSnYty+QffsC2bUryDWzBmDatMu0bWujuFjJ6NFlBAc3\n77nyntSoC76/qa0TplKpdHXNrKu7ZmlpKS+++CJLliwhMDAQrVZb5+uF8IWfRs1G0z+O2JTRrm0n\nXl2Hac8hbl+/EHBebP3kEy0//RTA/v0BfP99IOXlClQqB3Fx5cyadZlhw4x06CAXXb1JCr4b1dYJ\nMz4+ntTUVJ544ok6vz84OJjk5GQmT56MWq0mNjaWwYMHeyO6EPWm6R9Ht6VpHAFiU0Zz4tV1dFua\nxpaHF7Plv0P5z38C2b8/EIvFebLSo4eFp58uoV8/C3feaWn2d7v6Up3dMouLi1m9ejWXLl2ib9++\nDBkyxPXc0qVLmTZtmscDSrfMmvljJpBcDeWPueqT6dxLr3PbGyvYGvkQAwv+xaP8k50MRKl08Ktf\nldOnj4UuXcoZPNhEVJR7Crw//q7A/3LdcLfMv/71r7Rq1YpevXqxfv16zpw5w1NPPQXAxYsX3ZtS\nCOGfzGVo9/+LoK//TeB3magX2U5LAAAexUlEQVQOZRNzwdl87NGCNbwaPA3N3X146+F87r7bTGio\nnMH7qzoLfk5ODlOnTgWcwxLz5s1jy5YtPPDAA14JJ4TwMrsdDu8h7NP3Cfx+H+pDJ1EeK4ErsyEd\n4UrscS0517EDEf85wjdxwxjz49/oHR9D7P2j69638Lk6C77dbsdms6FSqQgJCWHatGnMnj2bdu3a\nycVEIRo7u53An78k6MBuAvd/jfrH46i+LwAgEHAEgCM2GOv97THddS/GpAcpj+3Nif/5kG5L0zj4\n53ncnjKaI1fG8CvG9IX/qrPgJyQk8PLLL/P//t//45ZbbiE6OprJkyfzyiuveCufEMIdyssJytxO\n4M/7CMr8isAfDuL4ERRXRl8cSqB9ELZuUfDovVzuMQBjwm8gqHo/edOeQxyZNstV3GNTRnPkynak\n4Pu1Ogv+73//ez7//HPKy6+2F42Li2P27Nl88MEHHg8nhGg45aWzBB3YSdCBrwj46UdUv2SjOFbq\nKu6ogHZgu6899ohISh98HGPvYRAcCTgvQhrruAhZMfWystiU0VLsG4E6C75CoWDQoEHVtsfGxjJr\n1iyPhRJC1ENpMZoDnxH449cEHjmI6thZlKcLUORdndvuCAZHbCjWUd2wRrXDNPi3GH91LwRH+DC4\n8JV6zcOfOnVqjWP2S5cudXsgdwpZuRJLjx5Y+vVzbQvcu5fAAwcomTDhhvaZk5PDuHHjuPXWW13b\nEhIS6NevH3v37r3uXPtrpaSkMGXKFNq3b+/adurUKV555RUcDgedOnXi+eeflxYLzZiyKJfAH/cS\nePRrAk4fRXXkNKqfL8HF8qtDMmqgdQC2+GisXTtjuf0OzAn3UH5LL1DKfzvCqV4Ff9y4ca6vrVYr\ne/fupVWrmvtg+BNLjx5EPvssBatWYenXj8C9e12Pb0aHDh149dVXq23v3LnzTe23wttvv81TTz1F\njx49WLx4MV9++SV33323W/Yt/JciP995AfX41wT8eAj1ydMozxWiyLFBpZmOjpZq7LdGYr2vA+Vx\n8Zi79cHc/degCa1950JQz4J/++23V3kcHx/PX/7yl1p7wFutVt544w1yc3MpLy/n4Ycfpnfv3q7n\nt2zZwo4dOwgLCwPg6aefrvNmgRtl6dePglWriHz2WcrGjkWXnu4q/u6WmZnp6nU/ZswY4uLiOHv2\nLJGRkcyfPx+TycSSJUsoLS2lqKiIYcOGMXz48Br3NX/+fFQqFeXl5eTn5xMZGen2vMIHHA6Uly6h\nOnUK3Xfr0P64GUVeMa3PB6DIKoeySi8NAqIDsXduifU3t1De5XbMcX2xdOmPIyTMZz+CaNxuqLVC\ncXExBQUFtT7/xRdfEBoayqRJkyguLmb69OlVCv6JEyeYOHEisbGxN/L2DWLp14+ysWMJffVVilNS\n3FLsT58+TUpKiuvxtdczsrOzWbZsGQaDgYkTJ3L06FHUajWDBg1iwIAB5OXlkZKSUmvBV6lU5OTk\nMG3aNIKDg2nXrt1NZxZeZDYTcPQo6sOHUR8/jio7G/XJk6hPnEBZVHT1dQogDIgph75ASyjtMwbz\nrQMw3z4AR6CcsQv3avAYvsPhIC8vr84eL3379qVPnz6ux9eOP588eZKNGzdSWFhIz549GTFixI1k\nr5fAvXvRpadTnJKCLj0dc1LSTRf9moZ0srKyXF+Hh4djMBgAMBgMWCwWDAYD69ev54svvkCn02G1\nWut8j+joaNasWcP//d//sXLlSmbOnHlTmYUblZejOn8e9enTqE6fRnXmDOpffnEW9jNnqhZ1wNa6\nNdaOHTE++CDWW24hxPK/qPT50AIIuvo6a1Abivr+l3d/FtGs1HsMPz8/n9LSUtq3b09wcDDKOpac\n0Wg0gHNxkGXLlvG73/2uyvNJSUncd9996HQ6lixZwnfffUevXr1q3Jder6/vz1KNYudO1BMmYP3g\nA4LuuQf70KG0+P3vsa5di+Oee+q9H7Va7cphNpsJCAiolis8PJygoCD0ej1KpdL1fFBQEOHh4Wza\ntIm77rqL0aNHs2/fPr755hv0ej0BAQFERERU2d+kSZOYNm0aHTp0oFWrVmi12mrvVzmTP2kSucxm\nOHcOxcWLKI4cQXH0KGRno/j5ZxQnT6K45tOtQ63G0aULtGuH4667sLZti6NLFxzx8dCpE6hUKICA\nK3+U61+q8W1V5vN+8btrEv+GXuSvuWpSr4L/zTff8Omnn7ra9TocDhQKBW+//Xat35OXl8fSpUsZ\nMmQI/fv3d213OBwMGzYMnc55Q0fPnj05efJkrQX/ZpoShezejWXlSixxcZCXB3FxBK5cSeDu3ZTU\n0Ku+NpWbIxUUFFBeXl4tV1FREWazmby8POx2u+t5s9lMUVERCQkJLF++nI8++oiwsDAUCgXnz5+n\nvLycwsLCKvt7+OGHmTFjBgEBAQQFBfHnP/+52vv5W8OmCn6fy25HWVCA6tQplIWFzqGWkydR5uej\nzM93nqnn5FT5XkdQEPaoKCy33oo9ORlbu3bYoqOxduiArUMHbDExUNssqhqGPg1BMajNWdW224Ji\n/OJ35/f/hn7G33LVdT20zm6ZFZ577jkWL15MaGj9xhQLCwuZP38+f/zjH4mPj6/yXFlZGVOnTmX5\n8uUEBQWxbNkyBg4cSM+ePWvcl3TLrJk/ZgLf51IYjaiOH0eVl4fy4kXnEEtODtrcXGwnT6LKzUV5\n+XKV73EEBWEzGLDr9Vg7dcJ6yy3YW7XCZjBg7dwZW7t2tRf0G6C9kEH40eko7UbXNrtSS1HX/8bY\nquaJEN7k63/D2kiu+rnhbpkVWrduTXBwcL3fcOPGjZSUlLBhwwY2bNgAwL333ovZbCY5OZnHHnuM\n+fPno1ariY+Pr7XYC1Ejux3ViRMEHjzoHEs/fBj12bOozp1DefEiCnvVpfFsej1ER2O99VYsSUlY\nb7kFW5s2zqLepg221q3x5qrYFUU99MRiVObz2IJiKI5N9YtiL5q2ep3hf/fdd3z00Ud0794dtfrq\nMWLUqFEeDQdyhl8bf8wEns2l2bKFiBkzUBYWVtlujYnBHhPjHGJp3975d+vWzqJuMOAIDm6Wv68b\n5Y+ZQHLV102f4X/00UfodDrKysqu/2Ih3M1uJ/K559Bu2gRA+W23Ufboo1h69cJ66604wmReuhD1\nUa+CbzabWbiwesMkIbwhptJ9CDmZmdhbtvRhGiEar3oNXMbExHD69GlPZxGiGkVpqevr86dPS7EX\n4ibU6ww/Ly+P1NRUDAYDAQEBru3+3jxNNH5Be/ZcfaC+oRvDhRBX1OsM//HHH2f27Nk888wz/PGP\nf3T9aQy0GRkYEhNp3bYthsREtBkZN7W/nJwcJtTRaTM/P5/ly5cDcODAAY4fPw7AnDlzGrzPjRs3\n8uyzzzJ+/Hi++uqrm8rdWNmjogCw3nKLb4MI0QTcUPO0xkKbkUH49Okojc75zuqsLMKnTwfAWEvj\nt5sVFRXFCy+8AMDWrVsZNGgQnTp1YsGCBQ3aT1FRER9//DFvv/02FouFJ598kj59+jS7pSUVV9oU\nqE+dAocDmtnPL4Q7NenPyKGLF7uKfQWl0Ujo4sVuKfgpKSl07tyZkydPUlZWxty5cwFYsGABkydP\nZt++fRw7dowOHTowfvx4MjIyyMzMJD09HQCTycTMmTOrDJNVCA8P55133nE1UgsJCWl2xR7AnJyM\nPSICe2SkFHshbpL37jbxAVUtc/hr234junXrxiuvvEKvXr3YsWOHa3vXrl1JTEzkmWeeqbJ2wKlT\np3jxxRdZtmwZSUlJ7Nq1q/b8KhUbN27kueee49e//rXbMjcqZrOrBYLub3/zdRohGrUmXfBttdyA\nUNv2G1Gx8lVFV8zr0ev1rFixgsWLF5OZmXndrpkjRoxg/fr1HDhwgP3797slc6MSFMSl1asBUNTj\n9yuEqF2TLvjFqanYtdoq2+xaLcWpqW57j7qGWRQKBfZrbvNfunQpM2bMIDU1lRYtWlDbjc5nzpxh\nzpw5OBwO1Go1gYGBzXJIB3AO5wDWjh19nESIxq1Jj+FXjNOHLl6M6vx5bDExFKemeuyC7bVuu+02\n3nrrLVq3bu3aNnjwYCZMmEBISAhRUVFcunSpxu9t3749nTp14rnnnkOhUJCYmMgdd9zhldz+RnVl\nrYGQVauw3HknDlkBTIgbUq9eOr4kvXRq5o+ZwDO5FJcvEz57Nrr168n9+GPKK62e5stc7uCPufwx\nE0iu+qqrl06THtIRTYMjLAzTkCHOr6+soyCEaDgp+KJRUFxp3Oe45pqMEKL+pOCLRkFx5X4KR1DQ\ndV4phKiNFHzRKFg7dQIg+IMPfJxEiMZLCr5oFCz9+mHt0AHVqVO+jiJEoyUFX/i3SjemOQID5eYr\nIW6CFHzht7Tr1hHdtSu6v//d2TgtIACk4Atxw5r0jVeicQt+7z0UVisRs2ah+fxzFCYTivJyX8cS\notGSgi/8kurkSQL37+fyrFk4goIIe+klFBYLtuhoX0cTotGSgi/8kvajj3AoFJQ99BD2mBjMffoQ\nOWkStkqdR4UQDSMFX/gfhwPtxo1Y+vTBfuU2cWv37uR+9plzLF8IcUM8UvCtVitvvPEGubm5lJeX\n8/DDD9O7Uv+Tb7/9lg0bNqBUKhk4cCDJycmeiCEaqYBDhwg4fpzCZ56p+oRCIYugCHETPFLwv/ji\nC0JDQ5k0aRLFxcVMnz7dVfCtViurV69m0aJFaDQaZs+eTe/evYmIiPBEFNEIaTMycAQEYLz/fl9H\nEaJJ8ci0zL59+zJ69GjXY5VK5fo6KyuL6OhoQkJCUKvVdO3alcOHD3sihvAmmw1FcbFb9qPdtAnT\nwIHSBlkIN/NIwddoNGi1WoxGI8uWLeN3v/ud6zmj0YiuUsdDrVZL2ZXGWKLxCnnjDaJvvx31qFEE\n7dwJ1yz8Ul+BX32FKicH44gR7g0ohPDcRdu8vDyWLl3KkCFD6N+/v2u7VqvFZDK5HhuNRoKDg2vd\nj16v91TEelOr1X6RozJ/y6QqKgKlEsXXX9Ni82YcsbHY/vQn7E88AS1a1H8///oXjpAQQn73O0Lc\n2ArZ335fFfwxlz9mAsnlDh4p+IWFhaSlpfHHP/6R+Pj4Ks+1adOG7OxsSkpK0Gg0HD58mAcffLDW\nffnDwgL+tsAB+F+m8JISNGFh2I4fp+S99whevZqgmTNxzJuH8be/pXTsWMp79qz7oqvZTPSGDRh/\n8xsKy8rAjZ/8/O33VcEfc/ljJpBc9VXXAigeKfgbN26kpKSEDRs2sGHDBgDuvfdezGYzycnJjB07\nlrS0NOx2OwMHDiQqKsoTMYQ32WygVkNgIKbhwzENH476yBGC09PRbtiAbv16yrt3p/SJJzCOGFHj\nQiaaHTtQXr4swzlCeIgscVgP/nYEB//LFD51KpqdO7GdPl0tl6KkBG1GBsHp6QQcPow9NJSyUaMo\nGzsWa5curtdFPv00gV9/zYXvvnMePNzI335fFfwxlz9mAslVX7LEofA4hc2Go5Yi7QgJoWzsWHK3\nbSP3o48wDR5M8Nq1GAYOpMWoUWg2bUKZn49m+3aMDz7o9mIvhHCS/7OEe1QM6dRFoaD8zjspvPNO\nLs+bh27dOnTvvUfU+PHYNRoUZjPGhx7yTl4hmiEp+MItFFYrKOv/gdHeogUlEyZQ8uyzBO3cSXB6\nOthszgu7QgiPkIIv3KOOIZ06KZWYBw3CPGiQ+zMJIaqQMXzhHjYbVLqjWgjhf6TgC7dQ2Gw4pOAL\n4dek4Av3qM9FWyGET0nBF+4hQzpC+D0p+MItFFarDOkI4eek4Av3kDN8IfyeFHzhHlLwhfB7UvCF\nWyis1hubhy+E8Bop+OLmmc0oTCY5wxfCz8kpmWgQRVkZ6p9+ci40fvAggQcPov75ZxTl5ZR37Spn\nEEL4MSn4olaKy5edhf1KcQ84dAj1L7+guLJ8oS0ykvJf/QrTPfdQHheH+de/pv5rWwkhvE0KvgBA\neenS1cJeUdxPnXI9b4uOpjw+HtOwYZTHx2OJi8MeE1P3ClZCCL8iBb+5cThQ5uQ4i/qPP7oKvLrS\nQjPW9u0pj4ujbPRoyuPjKY+Lw96ypQ9DCyHcQQp+M6L65Rf0o0ejyskBwKFQYO3UCctdd1EaF+cs\n7t2744iI8HFSIYQnSMFvRtRZWahycih5+mmM99+P9fbbcQQH+zqWEMJLpOA3I/YrZ+7mvn0pv/NO\nH6cRQnibzKJrRuzh4QAoi4p8nEQI4QtS8JsRKfhCNG9S8JsRR1gYDoUCZWGhr6MIIXxACn5zolLh\nCAtDIWf4QjRLHr1oe+zYMdauXcu8efOqbN+yZQs7duwgLCwMgKeffpqYmBhPRhFX2MPD5QxfiGbK\nYwX/448/Zvfu3Wg0mmrPnThxgokTJxIbG+uptxe1sEdESMEXopny2JBOq1atmDZtWo3PnTx5ko0b\nNzJ79mw2btzoqQiiBo7wcLloK0Qz5bGC36dPH1S1tMtNSkriT3/6E3PnzuXIkSN89913noohrmEP\nD0chZ/hCNEtev/HK4XAwbNgwdDodAD179uTkyZP06tWrxtfr9XpvxquRWq32ixyV3WgmVXQ0yn37\nPPbz+OPvCiRXQ/hjJpBc7uD1gm80Gpk6dSrLly8nKCiIQ4cOMXDgwFpfn5eX58V0NdPr9X6Ro7Ib\nzRSq0RCSn09ebq5HOl364+8KJFdD+GMmkFz1VdcEGK8V/D179mAymUhOTuaxxx5j/vz5qNVq4uPj\n6dmzp7diNHv2iAgUViuKsjLpoyNEM+PRgm8wGEhLSwOgf//+ru0DBgxgwIABnnxrUQtHxd22hYXY\npOAL0azIjVfNTEV7BblwK0TzIwW/manomClTM4VofqTgNzP2SkM6QojmRQp+M+OQM3whmi0p+M1M\nxZCOjOEL0fzIilfNhdlM0NdfE7RtGwDKy5d9HEgI4W1S8Jsw5aVLBH32GZpt2wjatQtlaSkOjQZT\ncjLGBx/0dTwhhJdJwW9KHA7UR4+i2bYNzbZtBHz/PQqHA1t0NMaHHsI0eDCW/v1xaLW+TiqE8AEp\n+I2d2UzQf/5D0LZtaLZvR332LACWX/2K4ilTMCcnUx4f75E2CkKIxkUKfiOkvHQJ5SefELlxo2uo\nxq7RYLn7bkomTcJ0773Yo6N9HVMI4Wek4DcGDgfqI0fQbN9eZahGER2NccQITMnJMlQjhLguKfj+\nqrahmh49KJ4yBe0jj5Dbtq0M1Qgh6k0Kvh9R5uU5Z9Vs337doRqNXg9+1JJVCOH/pOD7UsVQzZWz\n+CqzakaMwDR4MOZ+/UCGaoQQbiAF39vMZoK++oqgK+Px6nPngCtDNVOnYho8GGv37jJUI4RwOyn4\nXlDbUI15wABKJk/GNGiQzKoRQnicFHxPqDxUs20bAfv3y1CNEMLnpOC7y5WhGs22bQRt3y5DNUII\nvyMF/ya4hmoqetWUlVUdqrn3XuytWvk6phBCAFLwG8bhQH348NUboCqGalq3xvjww86hmqQkGaoR\nQvglKfjXYzaj+Pe/Cd+woepQzR13yFCNEKJRkYJfg5qGalRarXOoJiXFOatGhmqEEI2MFHy4OlRT\nMasmM/PqUM2oUQSOHEluXJwM1QghGrXmW/DNZoK+/NI5N37bNtRZWQBYEhKqDdXopY2BEKIJ8GjB\nP3bsGGvXrmXevHlVtn/77bds2LABpVLJwIEDSU5Odtt7ajMyCF28GNX589hiYihOTcU4ciQAytxc\ngnbsqDqrpmKo5oUXnLNqDAa3ZRFCCH/isYL/8ccfs3v3bjQaTZXtVquV1atXs2jRIjQaDbNnz6Z3\n795EXFlc+2ZoMzIInz4dpdEIgDori4hp0wj69FPUWVnVhmpcs2quySiEEE2R0lM7btWqFdOmTau2\nPSsri+joaEJCQlCr1XTt2pXDhw+75T1DFy92FfsKCrMZ3ZYtABRPm8bFTz/lwjffULRoEeZBg6TY\nCyGaDY+d4ffp04eLFy9W2240GtHpdK7HWq2WsrKyWvej1+vr/Z6q8+dr3O5QKOA//0ED3Eh5V6vV\nDcrhDf6YCSRXQ/ljLn/MBJLLHbx+0Var1WIymVyPjUYjwcHBtb4+rwEXSw0xMa6Lr5XZYmIatJ9r\n6fX6m/p+T/DHTCC5Gsofc/ljJpBc9RUTE1Prcx4b0qlNmzZtyM7OpqSkBKvVyuHDh+nSpYtb9l2c\nmor9mqmTdq2W4tRUt+xfCCEaM6+d4e/ZsweTyURycjJjx44lLS0Nu93OwIEDiYqKcst7VMzGqW2W\njhBCNGcKh8Ph8HWIupyvZVzem/ztIxv4ZyaQXA3lj7n8MRNIrvryqyEdIYQQviEFXwghmgkp+EII\n0UxIwRdCiGZCCr4QQjQTfj9LRwghhHvIGb4QQjQTUvCFEKKZkIIvhBDNRPNd8eoadrudt99+m9On\nTxMQEMCzzz5LdHS06/ktW7bw5ZdfApCQkMAjjzziF7kqXrN48WJ69+7NkCFD/CLX/v37Wb9+PQAd\nO3Zk3LhxKDy80Pv1Mm3atIm9e/eiVCoZMWIEiYmJHs1zLV8sCHQzufbs2cMnn3yCUqmkffv2PPXU\nUyiV3jlHrC1ThTfffJOQkBB+//vfeyXP9XL98ssvpKen43A4iIiIYNKkSQQGBno1W33IGf4V33zz\nDeXl5aSlpfH444+Tnp7ueu7ChQvs2bOHl156iZdeeokffviB06dP+zxXhX/84x+UlJR4JU99chmN\nRtasWcOMGTNIS0ujZcuWFBcX+zRTaWkpW7duJS0tjVmzZvH3v//d43kq+/jjj1m1ahXl5eVVtlcs\nCDRr1izmz5/PZ599RmFhoc9zWSwW1q1bx9y5c3nppZcoKyvj+++/92mmCtu2bePMmTNeyVJZbbkc\nDgdvvvkmEyZMYOHChdxxxx1+1WqhMin4Vxw5coQ77rgDgC5dunD8+HHXcy1atODFF19EqVSiVCqx\nWq0EBAT4PBfAf/7zH5RKpes13lJXrqNHj9KuXTvS09OZM2cO4eHhhIWF+TRTUFAQLVu2xGQyYTab\nPf5p41q+WBDoZnKp1WoWLlxIUFAQ4Pz05K3/5mvLBPDzzz9z7NgxBg8e7JUsldWWKzs7m9DQUP7v\n//6PuXPnUlJSUmc/G1+Sgn/FtQuzKJVKbDYb4PyPPywsDIfDQXp6Oh07dvTaP2hduc6cOcOePXt4\n9NFHvZKlvrmKi4v58ccfGTNmDC+++CKffPKJV5rg1ZUJnAfuKVOmMGPGDIYOHerxPJX16dMHlUpV\nbXtDFwTyVi6lUuladnTr1q2YTCZ+9atf+TRTQUEBH374IePGjfNKjmvVluvy5cscPXqU++67j9mz\nZ3Po0CEOHjzog4TXJ2P4V2i1WoyVlkd0OBxV/nEtFgtvvPEGWq2Wp556yi9y7d69m/z8fBYsWEBu\nbi5qtRqDweCVs/26coWGhtKpUydXwbjttts4deqUxw+SdWXKzMyksLCQ119/HYC0tDS6detG586d\nPZrpehq6IJA32e121qxZQ3Z2NlOnTvX6p6JrffXVV1y+fJlFixZRWFiI2WymTZs23HPPPT7NFRoa\nSnR0NG3btgWgR48enDhxgvj4eJ/mqokU/Cu6du3Kd999R1JSEj///DPt27d3PedwOFiyZAndu3fn\noYce8ptcY8aMcX39z3/+k4iICK8N7dSVKzY2lrNnz3L58mWCg4M5duwY9957r08zBQcHExgYSEBA\nAAqFguDgYEpLSz2e6XoqLwik0Wg4fPgwDz74oK9jAfDXv/6VgIAA/vznP3vtYm1d7r//fu6//34A\ndu7cSVZWls+LPTiHekwmEzk5OURHR3PkyBEGDRrk61g1koJ/RWJiIj/88AN/+ctfcDgcTJgwgS1b\nthAdHY3dbuenn36ivLyczMxMAB5//HG3rdR1o7l69+7t8fe/0VyPP/44aWlpAPTt27dK8fVVpoMH\nDzJr1iwUCgXdunXz2hBFTbyxINDN5IqNjeXzzz+nW7duLFiwAHAWXG/PbKqcyduzl66ncq7x48fz\n2muvAc7rRz179vRxuppJawUhhGgmfP85TQghhFdIwRdCiGZCCr4QQjQTUvCFEKKZkIIvhBDNhBR8\nIdystLSUadOmVWuDIYSvyTx8Idzo+++/Z/Xq1eTm5vo6ihDVSMEX4ooff/yR999/H71ez/nz59Fo\nNAwfPpytW7dy/vx57rrrLu68807effddXnnlFdf3VH68detWJk2axPLly335owhRIyn4QlRy/Phx\nnnrqKTp27MjLL7/MRx99xLx58ygrK+OZZ565bu+dWbNmeSmpEA0nBV+ISgwGAx07dgScPVJ0Op2r\nW6pOp/P6ugNCuJNctBWikmt7vl/bDrd9+/ZU7kZitVq9kksId5CCL0QDBAcHk5eXR1FREQ6Hg717\n9/o6khD1JkM6QjSAQqEgOTmZ1NRUIiMj6dmzp0y/FI2GdMsUQohmQoZ0hBCimZCCL4QQzYQUfCGE\naCak4AshRDMhBV8IIZoJKfhCCNFMSMEXQohmQgq+EEI0E/8fqjPCjcqAl30AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x154ef7b89e8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#plt.plot(param1[:,1],param1[:,2])\n",
    "\n",
    "## sequence 1\n",
    "plt.plot(param1[:,1],param1[:,2],color='orange',label=\"Path 1\")\n",
    "plt.plot(param1[Niter,1],param1[Niter,2],'x',color='orange',label=\"Final 1\")\n",
    "plt.plot(param1[0,1],param1[0,2],'o',color='orange',label=\"Initial 1\")\n",
    "## sequence 2\n",
    "plt.plot(param2[:,1],param2[:,2],color='blue',label=\"Path 2\")\n",
    "plt.plot(param2[Niter,1],param2[Niter,2],'x',color='blue',label=\"Final 2\")\n",
    "plt.plot(param2[0,1],param2[0,2],'o',color='blue',label=\"Initial 2\")\n",
    "## sequence 3\n",
    "plt.plot(param3[:,1],param3[:,2],color='red',label=\"Path 3\")\n",
    "plt.plot(param3[Niter,1],param3[Niter,2],'x',color='red',label=\"Final 3\")\n",
    "plt.plot(param3[0,1],param3[0,2],'o',color='red',label=\"Initial 3\")\n",
    "plt.xlabel(\"mu1\")\n",
    "plt.ylabel(\"mu2\")\n",
    "leg = plt.legend(loc='upper left')\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}