Stat_654_Homework_3.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "# Srujan Jabbireddy - Homework 3\n",
    "Due 11:55pm on Wed 4th March\n",
    "\n",
    "We study Table 13.2 on p240 of Lange and Chapter 13 Question 11 (p240) and Chapter 14 Example 14.10.1 (p264).\n",
    "\n",
    "The table contains counts of the number of days that a given number of deaths occurred. The data are from The London Times in the years 1910-1912."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy.special\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+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAdD0lEQVR4nO3df1CT9x0H8HcCBMRQqNVVezQBcazUbkNwW29KO725m157dRUlmKIbtJydnbZa6o9cKaMG8Cg77zyVYn+oNFNaZ3/ArrYynVZ662quuRYb9IoQdDpqK3YLMULMsz8cj6YCDcqTB/N9v+68y/N98n3y4Sv3vPk+T57n0UiSJIGIiISkVbsAIiJSD0OAiEhgDAEiIoExBIiIBMYQICISGEOAiEhgkWoXMBR2u13tEoiIbkqZmZn9tt9UIQAM/IMEw+l0Ii0tbRiruXlxLAJxPK7gWAQKh/EY7A9oHg4iIhIYQ4CISGCKHA7y+/0oKSnBsWPHoNPpsG7dOhiNRnl9Q0MDtm/fjoiICKSmpqKkpARarRZz585FXFwcACAxMRHl5eVKlEdERP+nSAg0Njaip6cHdXV1cDgcqKiowJYtWwAAXq8XGzZsQH19PUaNGoUVK1bgwIEDmD59OgCgtrZWiZKIiKgfihwOstvtyMrKAgCkp6ejublZXqfT6bBr1y6MGjUKAODz+RAdHY2WlhZcuHAB+fn5WLRoERwOhxKlERHRVRQJAbfbDb1eLy9HRETA5/Nd/kCtFmPHjgVw+a9+j8eDadOmISYmBgUFBXj55Zfxxz/+EU8//bTch4aPzWZDUlISJk+ejKSkJNhsNrVLIiIVKXI4SK/Xo7u7W172+/2IjIwMWK6srERbWxs2btwIjUaD5ORkGI1G+XVCQgLOnj2LCRMmBGzb6XRed11er/eG+t/sGhoaUFxcDK/XCwBwuVx49NFHcfr0aTzwwAMqV6cu0X83rsaxCBTu46FICGRkZODAgQOYM2cOHA4HUlNTA9YXFxdDp9Nh8+bN0GovT0Z2796N48ePo6SkBJ2dnXC73Rg3btw1276R7+uGw/d9b8Ts2bPlAOjj9XqxadMmFBUVqVTVyCD678bVOBaBwmE8Qn6dwKxZs6DT6WAymVBeXo41a9agvr4edXV1OHr0qLzDX7x4MfLy8rBv3z5kZ2fjv//9L3Jzc/HUU0+hrKwsYPZws+s7DKPValU7DNPR0TGkdiIKf4rsZbVaLUpLSwPaUlJS5NctLS399quqqlKiHNXZbDYUFhbC4/EAuHwYprCwEABgNptDVofBYIDL5eq3nYjExIvFQsBiscgB0Mfj8cBisYS0DqvVitjY2IC22NhYWK3WkNZBRCMHQyAERsphGLPZjJqaGvkEvNFoRE1NTUhnI0Q0sjAEQmCgwy1qHIYxm81ob2/H0aNH0d7ezgAgEhxDIAR4GIaIRiqGQAjwMAwRjVTh8x3MEc5sNnOnT0QjDmcCREQCYwgQEQmMIUBEJDCGABGRwBgCREQCYwgQEQmMIUBEJDCGABGRwBgCREQCYwgQEQmMIUBEJDCGABGRwBgCREQCYwgQEQmMIUBEJDCGABGRwBgCREQCYwgQEQmMIUBEJDCGABGRwBgCREQCYwgQEQmMIUBEJDCGABGRwBgCREQCYwgQEQmMIUBEJLBIJTbq9/tRUlKCY8eOQafTYd26dTAajfL6hoYGbN++HREREUhNTUVJSQkADNqHiIiGnyIzgcbGRvT09KCurg4rV65ERUWFvM7r9WLDhg3YsWMHdu3aBbfbjQMHDgzah4iIlKFICNjtdmRlZQEA0tPT0dzcLK/T6XTYtWsXRo0aBQDw+XyIjo4etA8RESlDkcNBbrcber1eXo6IiIDP50NkZCS0Wi3Gjh0LAKitrYXH48G0adPw7rvvDtjnak6n87rr8nq9N9Q/nHAsAnE8ruBYBAr38VAkBPR6Pbq7u+Vlv98fsDP3+/2orKxEW1sbNm7cCI1G8519+qSlpV13XU6n84b6hxOORSCOxxUci0DhMB52u33AdYocDsrIyMChQ4cAAA6HA6mpqQHri4uLcfHiRWzevFk+LPRdfYiIaPgpMhOYNWsWmpqaYDKZIEkSysrKUF9fD4/Hg3vuuQe7d+/G1KlTsXjxYgDAokWL+u1DRETKUiQEtFotSktLA9pSUlLk1y0tLf32+3YfIiJSFi8WIyISGEOAiEhgDAEiIoExBIiIBMYQICISGEOAiEhgDAEiIoExBIiIBMYQICISGEOAiEhgDAEiIoExBIiIBMYQICISGEOAiEhgDAEiIoExBIiIBMYQICISGEOAVGGz2ZCUlAStVoukpCTYbDa1SyISkiKPlyQajM1mQ2FhITweDwDA5XKhsLAQAGA2m9UsjUg4nAlQyFksFjkA+ng8HlgsFpUqIhIXQ4BCrqOjY0jtRKQchgCFnMFgGFI7ESmHIUAhZ7VaERsbG9AWGxsLq9WqUkVE4mIIUMiZzWbU1NTAaDRCo9HAaDSipqaGJ4WJVMBvB5EqzGYzd/pEIwBnAkREAmMIEBEJjCFARCQwhgARkcAYAkREAmMIEBEJjCFARCQwRa4T8Pv9KCkpwbFjx6DT6bBu3ToYjcaA91y4cAG/+93vYLVakZKSAgCYO3cu4uLiAACJiYkoLy9XojwiIvo/RUKgsbERPT09qKurg8PhQEVFBbZs2SKv/+yzz/Dcc8+hs7NTbrt48SIAoLa2VomSiIioH4ocDrLb7cjKygIApKeno7m5OWB9T08PNm3ahIkTJ8ptLS0tuHDhAvLz87Fo0SI4HI5hq6fvASaTJ0/mA0yIiK4S1EygqakJr776Knp6euS2HTt2DPh+t9sNvV4vL0dERMDn8yEy8vLHZWZmXtMnJiYGBQUFmD9/Ptrb2/HYY49h7969cp8+TqczmJJlDQ0NKC4uhtfrBXD5ASaPPvooTp8+jQceeGBI2wonXq93yGMZzjgeV3AsAoX7eAQVAuXl5Vi7di3Gjx8f1Eb1ej26u7vlZb/ff83O/NuSk5PlG4olJycjISEBZ8+exYQJEwLel5aWFlQNfWbPni0HQB+v14tNmzahqKhoSNsKJ06nc8hjGc44HldwLAKFw3jY7fYB1wV1OGjChAn4+c9/jokTJ8r/BpORkYFDhw4BABwOB1JTU7/zM3bv3o2KigoAQGdnJ9xuN8aNGxdMeYPiA0yIiAYW1EzgtttuQ3FxMe6++25oNBoAQE5OzoDvnzVrFpqammAymSBJEsrKylBfXw+PxzNgv+zsbKxZswa5ubnQaDQoKyv7ztlDMAwGA1wuV7/tRESiC2ovm5iYCAD46quvgtqoVqtFaWlpQFvf10CvdvU3gXQ6HaqqqoLa/lBYrdaAh5oDfIAJEVGfoA4HPfHEE7jnnnsQHR2Nu+66C0888YTSdQ0bPsCEiGhgQYVAVVUV9uzZg6ioKLz11ltYv3690nUNK7PZjPb2dhw9ehTt7e0MACKi/wvqcNDHH3+MXbt2AQAWL16MBQsWKFoUERGFRlAzAZ/PB7/fDwCQJEk+OUxERDe3oGYCc+bMQW5uLn784x/j008/xZw5c5Sui4iIQiCoEMjPz8f06dNx4sQJZGdnB/W9fyIiGvkGDYE33ngD8+fPR1VVlXwI6PPPPwcArFixQvnqiIhIUYOGQN9tIr59hTDPCRARhYdBTwz33Qn0s88+w29+8xv534cffhiS4oiISFmDzgRsNhu2bNmC8+fP4/3335fb+7v6l4iIbj6DhoDZbIbZbEZ1dTWWLFkSqpqIiChEgrpO4ODBg0rXQUREKgjqK6Lx8fHYvn07kpOTodVezo3p06crWhgRESkvqBC49dZb0dLSgpaWFrmNIUBEdPML+slix48fxxdffIHk5OSb/ik7RER0WVAhUFtbi4aGBvzoRz/CK6+8gtmzZ6OgoEDp2oiISGFBhUBDQwNsNhsiIyPR29sLk8nEECAiCgNBfTtIkiT5UY9RUVGIiopStCgiIgqNoGYCmZmZWLZsGTIzM2G32zFlyhSl6yIiohAIKgRWrVqFv//97zhx4gTmzZuH+++/X+m6iIgoBIIKga+//hqHDx9GW1sbzp49i/T0dMTHxytdGxERKSyocwJPPvkkUlJSUFRUhMTERDzzzDNK10VERCEQ1EwAAHJzcwEAd911F/bu3atYQUREFDpBzQQmTpyId955B52dndi/fz8SEhLQ1taGtrY2pesjIiIFBTUTOHHiBE6cOIHdu3dDkiQAQHFxMTQaDXbs2KFogUREpJygrxju6urCyZMnkZiYiDFjxihdFxERhUBQh4PeffddmEwmVFdXIycnB2+//bbSdRERUQgENRPYtm0b9uzZg9GjR8PtdmPx4sV46KGHlK6NiIgUFtRMQKPRYPTo0QAAvV6P6OhoRYsiIqLQCGomYDAYUFFRgalTp+LIkSMwGAxK10VERCEQ1Exg3bp1uPPOO/Hhhx/izjvvxPPPP690XUREFAJBzQSWLFmCV155RelaiELOZrPBYrGgo6MDBoMBVqsVZrNZ7bKIQiaomUBcXBwaGxvR2toa1EVifr8fxcXFyMnJQV5eHlwu1zXvuXDhAkwmE1pbW4PuQzScbDYbCgsL4XK5IEkSXC4XCgsLYbPZ1C6NKGS+cybgdrtx6tQpbN++XW77rovEGhsb0dPTg7q6OjgcDlRUVGDLli3y+s8++wzPPfccOjs7g+5DNNwsFgs8Hk9Am8fjgcVi4WyAhDFoCLz22mt45ZVXEBERgeXLl+O+++4LaqN2ux1ZWVkAgPT0dDQ3Nwes7+npwaZNmwJuRPddfYiGW0dHx5DaicLRoCHQ0NCAvXv3wu1245lnngk6BNxuN/R6vbwcEREBn88nP50sMzNzyH36OJ3OoGroj9frvaH+4YRjAYwfPx5nzpzpt13kseHvRqBwH49BQ0Cn00Gn02HMmDHo7e0NeqN6vR7d3d3yst/vv2Znfr190tLSgq7j25xO5w31DyccC6CyshKFhYUBh4RiY2NRWVkp9NjwdyNQOIyH3W4fcF1QJ4YByDeOC0ZGRgYOHToEAHA4HEhNTVWkD9GNMJvNqKmpgdFohEajgdFoRE1NDc8HkFAG/fP8iy++wMqVKyFJkvy6T1VV1YD9Zs2ahaamJphMJkiShLKyMtTX18Pj8SAnJyfoPkRKM5vNMJvNYfHXHtH1GDQENmzYIL82mUxBb1Sr1aK0tDSgLSUl5Zr31dbWDtqHiIiUNWgI/PSnPw1VHUREpIKgzwkQEVH4YQgQEQmMIUBEJDCGABGRwBgCREQCYwgQEQmMIUBEJDCGABGRwBgCREQCYwgQEQmMIUBEJDCGABGRwBgCREQCYwgQEQmMIUBEJDCGABGRwBgCREQCYwgQEQmMIUBEJDCGABGRwBgCREQCYwgQEQmMIUBEJDCGABGRwBgCREQCYwgQEQmMIUA0AthsNiQlJUGr1SIpKQk2m03tkkgQkWoXQCQ6m82GwsJCeDweAIDL5UJhYSEAwGw2q1kaCYAzASKVWSwWOQD6eDweWCwWlSoikTAEiFTW0dExpHai4cQQIFKZwWAYUjvRcFLknIDf70dJSQmOHTsGnU6HdevWwWg0yuv379+PTZs2ITIyEvPmzcOCBQsAAHPnzkVcXBwAIDExEeXl5UqURzSiWK3WgHMCABAbGwur1apiVSQKRUKgsbERPT09qKurg8PhQEVFBbZs2QIA6O3tRXl5OXbv3o1Ro0YhNzcXM2bMwC233AIAqK2tVaIkohGr7+SvxWJBR0cHDAYDrFYrTwpTSCgSAna7HVlZWQCA9PR0NDc3y+taW1thMBgQHx8PAMjMzMSRI0dwxx134MKFC8jPz4fP58OKFSuQnp6uRHlEI47ZbOZOn1ShSAi43W7o9Xp5OSIiAj6fD5GRkXC73fIhHwAYPXo03G43YmJiUFBQgPnz56O9vR2PPfYY9u7di8jIwBKdTud11+X1em+ofzjhWATieFzBsQgU7uOhSAjo9Xp0d3fLy36/X96Zf3tdd3c34uLikJycDKPRCI1Gg+TkZCQkJODs2bOYMGFCwLbT0tKuuy6n03lD/cMJxyIQx+MKjkWgcBgPu90+4DpFvh2UkZGBQ4cOAQAcDgdSU1PldSkpKXC5XDh//jx6enpw5MgRTJkyBbt370ZFRQUAoLOzE263G+PGjVOiPCIi+j9FZgKzZs1CU1MTTCYTJElCWVkZ6uvr4fF4kJOTg9WrV6OgoACSJGHevHm4/fbbkZ2djTVr1iA3NxcajQZlZWXXHAoiIqLhpcheVqvVorS0NKAtJSVFfj1z5kzMnDkzYL1Op0NVVZUS5RAR0QB4sRgRkcAYAkREAmMIEBEJjCFARCQwhgARkcAYAkREAmMIEBEJjCFARCQwhgARkcAYAkREAmMIEBEJjCFARCQwhgARAQBsNhuSkpIwefJkJCUlwWazqV0ShQDv1UxEsNlsAQ+7d7lcKCwsBAA+9jLMcSZARLBYLHIA9PF4PLBYLCpVRKHCECAidHR0DKmdwgdDgIhgMBiG1E7hgyFARLBarYiNjQ1oi42NhdVqVakiChWGABHBbDajpqYGRqMRGo0GRqMRNTU1PCksAH47iIgAXA4Cs9kMp9OJtLQ0tcuhEOFMgIhIYAwBIiKBMQSIiATGECAiEhhDgIhIYAwBIiKBMQSIaETpu5upVqvl3UxDgNcJENGIwbuZhh5nAkQ0YvBupqHHECCiEYN3Mw09hgARjRgj6W6mojxpjSFARCPGSLmbad+5CZfLBUmS5HMTagSB0ifKFQkBv9+P4uJi5OTkIC8vDy6XK2D9/v37MW/ePOTk5OD1118Pqg8Rhb+RcjfTkXJuIiRhJCngvffek1atWiVJkiR98skn0pIlS+R1PT090i9/+Uvp/Pnz0sWLF6WHH35Y+vLLLwft0+fIkSM3VNfnn39+Q/3DCcciEMfjCo6FJGk0GgnANf80Gk1I6zAajf3WYTQah7SdwfadiswE7HY7srKyAADp6elobm6W17W2tsJgMCA+Ph46nQ6ZmZk4cuTIoH2IiEJppJybCMWJckWuE3C73dDr9fJyREQEfD4fIiMj4Xa7ERcXJ68bPXo03G73oH2u5nQ6r7sur9d7Q/3DCcciEMfjCo4FsHTpUhQXF8Pr9cptMTExWLp0aUjHZvz48Thz5ky/7cNVhyIhoNfr0d3dLS/7/X55Z/7tdd3d3YiLixu0z9Vu5GEXfFjGFRyLQByPKzgWl/czd9xxBywWCzo6OmAwGGC1WkN+bqKysjLg4jng8onyysrKIf0f2e32AdcpcjgoIyMDhw4dAgA4HA6kpqbK61JSUuByuXD+/Hn09PTgyJEjmDJlyqB9iIhCzWw2o729HUePHkV7e7sqVyyH4kS5IjOBWbNmoampCSaTCZIkoaysDPX19fB4PMjJycHq1atRUFAASZIwb9483H777f32ISISXd9jP5WiSAhotVqUlpYGtKWkpMivZ86ciZkzZ35nHyIiUhYvFiMiEhhDgIhIYAwBIiKBMQSIiASmkSRJUruIYA32XVciIhpYZmZmv+03VQgQEdHw4uEgIiKBMQSIiAQW9iHA5xQE6u3tRVFRERYuXIjs7Gz87W9/U7sk1X399de4//770draqnYpqnvxxReRk5ODhx9+GG+88Yba5aiqt7cXK1euhMlkwsKFC8P29yPsQ6CxsRE9PT2oq6vDypUrUVFRoXZJqnrnnXeQkJCAP//5z9i6dSuef/55tUtSVW9vL4qLixETE6N2Kar76KOP8Mknn2Dnzp2ora3Fv//9b7VLUtXBgwfh8/mwa9cuLF26FBs2bFC7JEWEfQjwOQWBfv3rX2P58uXyckREhIrVqG/9+vUwmUz43ve+p3Ypqjt8+DBSU1OxdOlSLFmyBL/4xS/ULklVycnJuHTpEvx+P9xud793NQ4H4flTXSXY5xSIYvTo0QAuj8uyZcvw5JNPqlyRevbs2YMxY8YgKysLNTU1apejuq6uLpw+fRrV1dU4deoUHn/8cezduxcajUbt0lQRGxuLf/3rX5g9eza6urpQXV2tdkmKCPuZQLDPKRDJmTNnsGjRIjz00EN48MEH1S5HNX/5y1/w4YcfIi8vD06nE6tWrcLZs2fVLks1CQkJmD59OnQ6HSZOnIjo6GicO3dO7bJUs23bNkyfPh3vvfce3n77baxevRoXL15Uu6xhF/YhwOcUBPrqq6+Qn5+PoqIiZGdnq12Oqmw2G1577TXU1tYiLS0N69evx7hx49QuSzWZmZn44IMPIEkSOjs7ceHCBSQkJKhdlmpuueUW+SmI8fHx8Pl8uHTpkspVDb+w/5OYzykIVF1djf/85z/YvHkzNm/eDADYunUrT4wSZsyYgY8//hjZ2dmQJAnFxcVCnzP67W9/i7Vr12LhwoXo7e3FU089hdjYWLXLGna8YpiISGBhfziIiIgGxhAgIhIYQ4CISGAMASIigTEEiIgExhCgm8JHH32EqVOn4syZM3LbCy+8gD179lz3Nk+dOoUFCxYMR3nXuHTpEgoKCpCbm4tvvvlGbl+9ejUefPBB5OXlITc3F7///e9x8uTJIW///PnzqK+vl7fZdy0M0VAxBOimERUVhTVr1uBm+Fbz2bNn0dXVhZ07dyI+Pj5gXVFREWpra7Fz507k5+df1607jh07hv379w9XuSSwsL9YjMLHvffeC7/fD5vNhkceeURuP3XqFFasWIHXX38dALBgwQL86U9/wptvvgmXy4Wuri588803WLhwId5//320tbVh/fr1GDt2LM6dO4clS5bg3LlzuP/++7F06VKcOXMGzz77LC5evIjo6Gg8//zzuHTpEh5//HEkJCTgvvvuw2OPPSZ//jvvvIPt27dDp9MhKSkJpaWlePbZZ9He3o7i4mKUlpYO+DNNnToVUVFRcLlc0Ol013zuhAkTUFVVhebmZnR3dyMlJQXl5eWorq5GS0sL6urqAAB1dXV46aWX4Ha7UVJSgh/84AdYvnw53G43vF4vioqK8LOf/Uyh/xm6mXEmQDeVkpISbNu2De3t7UG9PyYmBi+//DJ+9atf4eDBg6iurkZhYSH++te/AgA8Hg8qKyuxc+dOfPDBB2hpacH69euRl5eH2tpaFBQU4IUXXgBw+a/7l19+OSAAurq6sHHjRmzfvh07d+5EXFwc6urq8Nxzz2HSpEmDBkCf2267DV1dXf1+rtvtxi233IJXX30Vu3btgsPhQGdnJ5YsWYJ7770XOTk5AIDJkydjx44deOSRR7Bnzx50dHTgq6++QnV1NaqqquD1eoc40iQKzgTopnLrrbdi7dq1WL16NTIyMvp9z9WHi+6++24AQFxcHCZNmgTg8n1g+m4Edtddd8n3h/nhD3+ItrY2HD9+HC+++CJeeuklSJKEqKgoAEBiYiJ0Ol3AZ508eRKTJk2S71T7k5/8BIcPHx7SbZhPnz6N8ePH9/u5fTdxW7FiBWJjY+HxeNDb23vNNiZPngwAGDt2LLxeL77//e/DbDZjxYoV8Pl8yMvLC7oeEgtDgG46M2fOxL59+/Dmm2+iqKgI0dHR+Prrr3Hp0iV0d3fj1KlT8nu/6zbIra2t6O7uRnR0ND799FPk5ORg4sSJyM/PR0ZGBlpbW/Hxxx8DALTaayfOiYmJaG1thcfjQWxsLP75z38iOTk56J+lqakJMTExGD9+fL+fe+jQIZw5cwYbNmzAuXPnsG/fPkiSBK1WC7/fP+DPeezYMXR3d6OmpgZffvklTCYTZsyYEXRdJA6GAN2ULBYL/vGPfwAAxo0bh2nTpiE7OxsGgwFGozHo7cTHx+Opp57CuXPnMGfOHEyaNAmrVq1CSUkJLl68CK/XC4vFMmD/MWPG4A9/+AMWLVoErVYLg8GAp59+etBbUldWVmLr1q3QarUYPXq0/MSq/j43MTERmzdvxoIFC6DT6XDnnXfiyy+/hMFgwPHjx7Ft27Z+PyMpKQmbNm3CW2+9haioKCxbtizoMSGx8AZyREQC44lhIiKBMQSIiATGECAiEhhDgIhIYAwBIiKBMQSIiATGECAiEhhDgIhIYP8D9+lgnjzJze4AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "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\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Problem 1\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  𝜆  (this should 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$$ f(x|\\lambda) = \\frac{e^{-\\lambda}\\lambda^x}{x!}$$\n",
    "\n",
    "The MLE maximizes the log likelihood. \n",
    "\n",
    "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",
    "\n",
    "Differentienting, setting equal to $0$, and solving for $\\lambda$, we obtain$$ \\widehat{\\lambda} = \\bar{X} = \\frac{1}{n}\\sum_{i=1}^n x_i $$In other words, the MLE is the sample mean.\n",
    "\n",
    "We now compute the sample mean for this data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([  0, 267, 542, 555, 444, 305, 162,  56,  24,   9])"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.arange(freq.size)*freq"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2.156934306569343"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# mean number of deaths = total deaths / number of days\n",
    "death_mle = np.sum(np.arange(freq.size)*freq)/(np.sum(freq))\n",
    "death_mle"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYEAAAEBCAYAAACe6Rn8AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3dfVzUdb7//8cMMCiikBeZrYJEa1qZhta6a9Ypj23Z6dSeCFAW22TlZxfbheZmkUQqiGuau66XlWk0KeUxC3dzS2uz7GTJLRI30EIEXS/CC+w3jFzOfP+YGJ0EHIxhkM/z/g/M+zPvmde8xXnO+/OZz/tjcjqdTkRExJDM/i5ARET8RyEgImJgCgEREQNTCIiIGJhCQETEwBQCIiIGFujvAloiLy/P3yWIiFyQhg0b1mj7BRUC0PQL8UZhYSGDBg1qxWouXBoLTxqP0zQWnjrCeDT3AVq7g0REDEwhICJiYAoBEREDUwiIiBiYQkBExMAUAgZjLbDSf2F/rnrjKvov7I+1wOrvkkTEjy64r4jK+bMWWEnJTcFeaweg9GQpKbkpACQOTvRnaSLiJ5oJGEjqllR3ADSw19pJ3ZLqp4qk3SmxEr1rNLxuhg39oeSnzxT379/PH/7wB5KSkkhISCA9PR2bzXbW/QoLC/nrX//a5ONs3bqVnJycFj//LbfcQnV1tUdbUlISsbGx7pqeeOIJTpw40ezj5OTkUFtb2+Lnb+80EzCQspNlLWoXgymxwucpWOp/+KBgL4XPXTNFos5vplhVVcWDDz7I7NmzGTJkCABvvfUWU6dOZfny5R73HTRoULMnZd14443nVUNT5s6dS3R0NADvvPMOaWlpLFq0qMn7L1++nLvvvrtVa2gPFAIGEhEWQenJ0kbbRfgqFeo9Z4rU213t5xkC//znP7nuuuvcAQDwm9/8hjVr1rB//34WL15MRUUFFRUVJCcn8/e//50XXniBN998E6vVSlhYGEFBQYwdOxaAvXv3kpCQwNSpU7nkkkvYv38/gwcP5rnnnuPw4cOkp6dTXV1NRUUFDz30EP/5n//pVZ3//d//zcKFC6muruarr75yz0iqqqpISUlh165dlJeX8/jjj7No0SLS0tI4fPgwJ06c4MYbb+Sxxx47r/FpD7Q7yEAyRmcQEhTi0RYSFELG6Aw/VSTtir2JGWFT7V7Yv38/ERFnf8jo27cvBw8eBGDEiBGsXbuWbt26AXD8+HFeeukl1qxZw8qVKzl16tRZ/fft20dGRgZvvvkmW7dupby8nL1793L//ffzyiuvMGPGDKzWlu3K6tatG99//z3ffPMN8+bN49VXX+WWW27h008/5d5776VXr1688MILHDp0iKFDh/Lyyy+zZs0a1qxZcx4j035oJmAgDQd/U7ekUnayjIiwCDJGZ+igsLiERLh2ATXWfp569+7Nzp07z2rft28fl156KQBRUVEe28rKyoiOjqZz584AXHvttWf1j4iIIDQ0FIBevXpRXV1Nr169WLp0KevWrcNkMlFXV+d1nU6nk6NHj9KjRw969+5NRkYGISEhHDlyhH79+nncNzw8nIKCAj777DNCQ0Opqanx+nnaI80EDCZxcCL7HtvHv+L+xb7H9ikA5LQhGRDgOVMkIMTVfp5Gjx7Np59+6hEEb775Jt27d3e/uZpMJo8+ERER7N27l6qqKhwOR6Mh8uM+AH/+85+56667mDdvHr/4xS9wOp1e17lu3TpGjBiB2WzmmWeeITMzk6ysLC6++GKP53Q4HKxfv56uXbsyf/58Jk6cSFVVVYueq73RTEBEXH7Y71+zYxqW2sOuGcCQjPM+HgDQpUsXli1bRmZmJhUVFdTX13PFFVewYMGCJvt0796dSZMmMX78eMLDw6muriYwMPCcn+xvu+02MjIyWL58OX369Dnnt32efPJJ92yjd+/ePPvsswDcddddxMXF0a1bN3r27EllZSUAw4cPJyUlhbS0NKZMmUJeXh6dO3cmMjKS7777jt69e7dkaNoP5wVkx44dP6n/119/3UqVXPg0Fp40Hqf5eyxqa2udS5Yscd8eP3688/PPP/dbPf4ej9bQ3HunZgIi0q4EBgZy6tQpfvOb3xAUFMQ111zD8OHD/V1Wh6UQEJF2Z8qUKUyZMsXfZRiCDgyLiBiYQkBExMB8sjvI4XCQnp7O7t27sVgszJ49m8jISPf2jRs3snr1agICAhgwYADp6emYzWbuvvtuunbtCrhOJpkzZ44vyvMLa4FV388XkXbHJyGwefNmampqyMnJIT8/n6ysLJYuXQq4TsNeuHAhubm5dO7cmSlTpvDhhx9yww03AJCdne2LkvxKq3eKSHvlk91BeXl5jBo1CoChQ4eya9cu9zaLxcLatWvd38+tq6sjODiYoqIiTp06xcSJE5kwYQL5+fm+KM0v2tXqnSVW2NCfgV9e1WqrRErHYS2wMnrjaMzPmVvlehPbt2/nl7/8JUlJSSQlJREXF9fsB73y8nLS09N/0nN6a9GiRfz6178mKSmJ8ePHM3HiRL7++msA1q9fz5YtW5rs29SKoo8//jg1NTVMnz6drVu3elVHRUUFubm5AKxYsaLRk+N8ySczAZvN5j6lGyAgIIC6ujoCAwMxm8307NkTcH3qt9vtjBw5kj179pCcnMy9997Lvn37mDRpEps2bSIw0LPEwsLC866rqqrqJ/U/X82t3tmW9XQ7vpE+ZWmYnVWYAOylOD77PYcOHuT77v/VZnW0R/7622hPNpZuJG1HGlX1VYBrxvr7t3/PwX8f5L8iz+/vo7S0lCuvvJInnngCgNraWh588EEGDhzo8R5xpvj4+Db5tygvL+e2227jtttuA+DAgQM8/PDDvPDCC+7VTAsLCxv921i0aBEDBw7EYrF4tKekpFBcXExFRQX79+/36nUUFBTwj3/8g8svv9z94bkt/xZ9EgKhoaHus+zAdYzgzDdzh8PBvHnzKCkpYdGiRZhMJqKiooiMjHT/Hh4eTnl5OX369PF47OaWmj2XwsLCn9T/fDW3emeb1rPhdnBWeTSZnVX8rHwxPxs5re3qaIf89bfRntz+j9vdAdCgqr6KxUWLmXbb+f19fP/993Tr1s09tsePH6dTp05ceeWVlJaWMmvWLAICAggODmbWrFk4HA6mTJnCG2+8wQsvvMBnn32Gw+Hgjjvu4He/+x1Wq5UNGzZgNpuJiYnhySef5MCBA6SmplJXV4fJZOKZZ55h4MCB3HrrrcTExFBSUkKPHj1YtGgRAQEB7tp69epFz5493bUNGjSITZs2UVlZydatW+nZsye//vWvSUlJoXPnztTW1vLcc8+xc+dOTp48yfLly7nvvvt4/vnnCQoKIi4ujr/85S+8++67hIeH8/HHH/Pee+9RX19PRkYGAQEB7tcGEBcXx4IFC3j33XcpKipi586dfPnll4wdO5Zf/vKXPP300+zfv5/6+nruv/9+xo4dS1JSEgMHDuSbb77BZrPx5z//mZ/97Gfn/HfIy8trcptPdgfFxMS4p0L5+fkMGDDAY3taWhrV1dUsWbLEvVto3bp1ZGVlAXDkyBFsNhu9evXyRXltLmPwWEJ+tNRJiMnV3qZ8sEqkdBy+ut7EZ599RlJSEhMmTGDatGnMmDGDLl268Mwzz5CWlsZrr73GuHHj3P//G2zYsIHnn38eq9VKp06dANdumtTUVHJycujXrx91dXX86U9/IikpCavVSmpqKk8//TTgWsH00UcfJScnh+PHj1NQUHDOWnv06OGx3MTOnTsJCQnhxRdf5JlnnsFms3msKApQXV3N66+/fta1BmJiYli9ejWTJk1i3rx5TT7n5MmTGTFiBPHx8e62nJwcLrroItauXcsrr7zCwoULOX78OADXXHMNq1atYuTIkfztb38752s6F5/MBMaMGcO2bdtISEjA6XSSmZlJbm4udrudq6++mnXr1jF8+HDuu+8+ACZMmEBsbCxPPfUU48aNw2QykZmZedauoAtV4v//d7gYUo9BWR1EBEJGjx/a25IPVomUjsNX15sYMWKE+w3zTN999537U/h1113H/PnzPbYvWLCABQsWcPToUfdukjlz5rBy5Uqef/55hg4ditPppLi4mOuuuw5wfZo/fPgwABdddJF7T0KfPn3OurpYYw4ePMitt97K3r17AdeFbLZv386DDz5IYGAgDzzwwFl9frwKaoOGs5yvvfZa/vSnP5213dnMonPFxcX86le/Alx7VqKjo9m/fz8AV155JQCXXHIJR48ePedrOhefvMuazWZmzpzp0dZwBR+AoqKiRvv9+I+gw7CXkdgNErud3d6mhmS4rhR15oVDfuIqkdJxZIzO8PgWG/j2ehMXX3wxRUVFDBw4kC+++IL+/fu7t9XU1LBp0yYWLFiA0+nkjjvu4I477uCNN97gueeeIzg4mOTkZL788kuio6PZsWMHo0ePprCw0H3MsbGVRpuzZ88evv32W4YOHcr//d//Aa4D2927d2flypV8+eWXLFiwgOzsbPeKouB6v2vMzp07iYmJYceOHfz85z8nODiYY8eOUV9fT2VlJQcOHHD3b3isBg2vacyYMdhsNvbs2UPfvn1b9Hq81TE+ard37eUTeMNqkF+l4rSXYWqFVSKl42j4uvK0TdM4bD/s8/NZZs+ezaxZs3A6nQQEBJCZmeneZrFYCAsL46677iIsLIyRI0dy6aWXcsUVVxAbG8tFF11E7969GTJkCJdeeikzZsxg5cqV1NXVkZHhfWitWrWKv//975jNZgIDA/nLX/7isQdi4MCBLFiwgM2bN2M2m3nooYeA0yuKNtxuzFdffcWECRPcezZ69erFyJEjiY2NJSIiwn3uVEREBHv27GHVqlXuvnFxccyYMYNx48ZRXV3Nww8/TI8ePbx+XS1hcjY3J2ln8vLyGDZs2Hn399vBvx+u3XrWJ/DrV/jtDVgHQj1pPE7TWHjqCOPR3Hunlo1oC1GJrjf8kEjA5PrpxwAQEWmg3UFtJSpRb/oi0u5oJiAiYmAKARERA1MIiIgYmEJARMTAFAIiIgamEBARMTCFgIiIgSkEREQMTCEgImJgCgEREQNTCIiIGJhCQETEwBQCIiIGphAQETEwhYCIiIEpBEREDEwhICJiYAoBEREDUwiIiBiYQkBExMAUAiIiBqYQEBExMIWAiIiBKQRERAws0BcP6nA4SE9PZ/fu3VgsFmbPnk1kZKR7+8aNG1m9ejUBAQEMGDCA9PR0gGb7iIhI6/PJTGDz5s3U1NSQk5PD1KlTycrKcm+rqqpi4cKFvPrqq6xduxabzcaHH37YbB8REfENn4RAXl4eo0aNAmDo0KHs2rXLvc1isbB27Vo6d+4MQF1dHcHBwc32ERER3/BJCNhsNkJDQ923AwICqKurcz2h2UzPnj0ByM7Oxm63M3LkyGb7iIiIb/jkmEBoaCiVlZXu2w6Hg8DAQI/b8+bNo6SkhEWLFmEymc7Zp0FhYeF511VVVfWT+nckGgtPGo/TNBaeOvp4+CQEYmJi+PDDDxk7diz5+fkMGDDAY3taWhoWi4UlS5ZgNpu96tNg0KBB511XYWHhT+rfkfh7LKwFVlK3pFJ2soyIsAgyRmeQODjRb/X4ezzaE42Fp44wHnl5eU1u80kIjBkzhm3btpGQkIDT6SQzM5Pc3FzsdjtXX30169atY/jw4dx3330ATJgwodE+0jFZC6yk5KZgr7UDUHqylJTcFAC/BoGIEfkkBMxmMzNnzvRoi46Odv9eVFTUaL8f95GOKXVLqjsAGthr7aRuSVUIiLQxnSwmba7sZFmL2kXEdxQC0uYiwiJa1C4ivqMQkDaXMTqDkKAQj7aQoBAyRmf4qSIR41IISJtLHJzIijtXEBkWiQkTkWGRrLhzhY4HiPiBTw4Mi5xL4uBEvemLtAOaCYiIGJhCQETEwBQCIiIGphAQETEwhYCIiIEpBEREDEwhICJiYAoBEREDUwiIiBiYQkBExMAUAiIiBqYQEP8oscKG/vC62fWzxOrvikQMSQvISdsrscLnKVD/w9XF7KWu2wBRWlROpC1pJiBt76vU0wHQoN7uaheRNqUQkLZnb+Iykk21i4jPKASk7YU0cRnJptpFxGcUAtL2hmRAgOflJQkIcbWLSJtSCEjbi0qE61dASCRgcv28foUOCov4gVffDtq2bRuvvPIKNTU17rZXX33VZ0WJAUQl6k1fpB3wKgTmzJnD008/zSWXXOLrekREpA15FQJ9+vThV7/6la9rERGRNuZVCPTo0YO0tDSuvPJKTCYTAPHx8T4tTEREfM+rEOjbty8AR48e9WkxIiLStrz6dtDDDz/M1VdfTXBwMAMHDuThhx/2dV0iItIGvAqB+fPns379eoKCgtiwYQNz585t9v4Oh4O0tDTi4+NJSkqitLT0rPucOnWKhIQEiouL3W133303SUlJJCUl8dRTT7XwpYiISEt5tTvoiy++YO3atQDcd999xMXFNXv/zZs3U1NTQ05ODvn5+WRlZbF06VL39oKCAp599lmOHDnibquurgYgOzu7xS9CRETOj1czgbq6OhwOBwBOp9N9cLgpeXl5jBo1CoChQ4eya9cuj+01NTUsXryYyy67zN1WVFTEqVOnmDhxIhMmTCA/P79FL6RZPyxbPPDLq7RssYjIGbyaCYwdO5Zx48YxZMgQdu7cydixY5u9v81mIzQ01H07ICCAuro6AgNdTzds2LCz+nTq1Ink5GTuvfde9u3bx6RJk9i0aZO7T4PCwkJvSnbrdnwjfcrSMDurMAHYS3F89nsOHTzI993/q0WP1ZFUVVW1eCw7Mo3HaRoLTx19PLwKgYkTJ3LDDTewd+9eYmNjGTBgQLP3Dw0NpbKy0n3b4XCc9Wb+Y1FRUURGRmIymYiKiiI8PJzy8nL69Onjcb9BgwZ5U/JpG24HZ5VHk9lZxc/KF/OzkdNa9lgdSGFhYcvHsgPTeJymsfDUEcYjLy+vyW3N7g568803AdeB4Y0bN/L111+zceNGFixY0OwTxsTEsHXrVgDy8/PPGRoA69atIysrC4AjR45gs9no1avXOfudk5YtFhFpUrMfzxuWiThz3z1wzmMCY8aMYdu2bSQkJOB0OsnMzCQ3Nxe73d7kSWaxsbE89dRTjBs3DpPJRGZm5jlnD14JiXBduaqxdhERg2v2Xbbh4G5BQQFpaWnu9j/+8Y/cfffdTfYzm83MnDnToy06Ovqs+535TSCLxcL8+fO9q7olhmR4XsoQtGyxiMgPmg0Bq9XK0qVLqaio4L333nO3N/aG3m41rFT5VSpOexmmkAhXAGgFSxGR5kMgMTGRxMREli1bxuTJk9uqptb3w7LFRR3gAI+ISGvy6jyBjz76yNd1iIiIH3h15DUsLIzVq1cTFRWF2ezKjRtuuMGnhYmIiO95FQIXXXQRRUVFFBUVudsUAiIiFz6vryy2Z88evv32W6KiorRfXUSkg/AqBLKzs9m4cSPXXHMNK1eu5Pbbbyc5OdnXtYmIiI95FQIbN27EarUSGBhIbW0tCQkJCgERkQ7Aq28HOZ1O99m7QUFBBAUF+bQoERFpG17NBIYNG8YjjzzCsGHDyMvL49prr/V1XSIi0ga8CoEnn3ySf/7zn+zdu5d77rmHm266ydd1iYhIG/AqBI4dO8Ynn3xCSUkJ5eXlDB06lLCwMF/XJiIiPubVMYHHHnuM6Ohopk2bRt++ffnjH//o67pERKQNeL1W87hx4wAYOHAgmzZt8llBIiLSdryaCVx22WW88847HDlyhA8++IDw8HBKSkooKSnxdX0iIuJDXs0E9u7dy969e1m3bh1OpxOAtLQ0TCYTr776qk8LFBER3/H6jOETJ06wf/9++vbtS/fu3X1dl4iItAGvdge9++67JCQksGzZMuLj43n77bd9XZeIiLQBr0Jg1apVrF+/niVLlvDWW29dcLuArAVW+i/sz1VvXEX/hf2xFlj9XZKISLvg1e4gk8lEly5dAAgNDSU4ONinRbUma4GVlNwU7LWuawyXniwlJTcFgMTBusSkiBibVyEQERFBVlYWw4cPZ8eOHURERPi6rlaTuiXVHQAN7LV2UrekKgRExPC82h00e/Zs+vXrx6effkq/fv2YNWuWr+tqNWUny1rULiJiJF7NBCZPnszKlSt9XYtPRIRFUHqytNF2ERGj82om0LVrVzZv3kxxcfEFd5JYxugMQoJCPNpCgkLIGJ3hp4pERNqPc84EbDYbBw4cYPXq1e62C+kksYb9/qlbUik7WUZEWAQZozN0PEBEhHOEwGuvvcbKlSsJCAjg0Ucf5cYbb2yrulpV4uBEEgcnUlhYqOsji4icodndQRs3bmTTpk3k5ORcMJ/8RVpC55CI0TU7E7BYLFgsFrp3705tbW1b1STSJnQOiYiXB4YB98JxIh1Fc+eQiBhFszOBb7/9lqlTp+J0Ot2/N5g/f36T/RwOB+np6ezevRuLxcLs2bOJjIz0uM+pU6e4//77ycjIIDo62qs+Iq2prJGvDjfXLtIRNRsCCxcudP+ekJDg9YNu3ryZmpoacnJyyM/PJysri6VLl7q3FxQU8Oyzz3LkyBGv+4i0toigAEpr6xttFzGKZkPg+uuvP68HzcvLY9SoUQAMHTqUXbt2eWyvqalh8eLFHpepPFcfkdaW0b2elO/AfsaezhCTq13EKLy+vGRL2Gw2QkND3bcDAgKoq6sjMND1dMOGDWtxnwaFhYXnXVdVVdVP6t+RaCzg3h59gEOkHoOyOogIhIwernYjj43+Njx19PHwSQiEhoZSWVnpvu1wOM56Mz/fPj/le/46T+A0jQXQaR6Jn6eQ2O2Mg8MBITB8HoOijDs2+tvw1BHGIy8vr8ltXn87qCViYmLYunUrAPn5+QwYMMAnfUR+kqhEuH4FhETixAQhka7bUfp6qBiHT2YCY8aMYdu2bSQkJOB0OsnMzCQ3Nxe73U58fLzXfUR8LioRohIp6gCf9kTOh09CwGw2M3PmTI+26Ojos+6XnZ3dbB8REfEtn+wOEhGRC4NCQETEwBQCIiIGphAQETEwhYCIiIEpBEREDEwhICJiYAoBEREDUwiIiBiYQkBExMAUAiIiBqYQEBExMIWAiIiBKQRERAxMISAiYmAKARERA1MIiIgYmEJARMTAFAIiIgamEBARMTCFgIiIgSkEREQMTCEgImJgCgEREQNTCIiIGJhCQKQ9KLHChv7wutn1s8Tq74rEIAL9XYCI4ZVY4fMUqLe7bttLXbcBohL9V5cYgmYCIv72VerpAGhQb3e1i/iYQkDE3+xlLWsXaUU+2R3kcDhIT09n9+7dWCwWZs+eTWRkpHv7Bx98wOLFiwkMDOSee+4hLi4OgLvvvpuuXbsC0LdvX+bMmeOL8kTal5AI1y6gxtpFfMwnIbB582ZqamrIyckhPz+frKwsli5dCkBtbS1z5sxh3bp1dO7cmXHjxnHzzTfTrVs3ALKzs31Rkkj7NSTD85gAQECIq13Ex3yyOygvL49Ro0YBMHToUHbt2uXeVlxcTEREBGFhYVgsFoYNG8aOHTsoKiri1KlTTJw4kQkTJpCfn++L0kTan6hEuH4FhEQCJtfP61fooLC0CZ/MBGw2G6Ghoe7bAQEB1NXVERgYiM1mc+/yAejSpQs2m41OnTqRnJzMvffey759+5g0aRKbNm0iMNCzxMLCwvOuq6qq6if170g0Fp78Px4xcMW7p29WAX6qx/9j0b509PHwSQiEhoZSWVnpvu1wONxv5j/eVllZSdeuXYmKiiIyMhKTyURUVBTh4eGUl5fTp08fj8ceNGjQeddVWFj4k/p3JBoLT/4eD2uBldQtqZSdLCMiLIKM0RkkDvbPTMDfY9HedITxyMvLa3KbT3YHxcTEsHXrVgDy8/MZMGCAe1t0dDSlpaVUVFRQU1PDjh07uPbaa1m3bh1ZWVkAHDlyBJvNRq9evXxRnki7Yi2wkpKbQunJUpw4KT1ZSkpuCtYCnTAmvueTmcCYMWPYtm0bCQkJOJ1OMjMzyc3NxW63Ex8fz/Tp00lOTsbpdHLPPffQu3dvYmNjeeqppxg3bhwmk4nMzMyzdgWJdESpW1Kx13qeJ2CvtZO6JdVvswExDp+8y5rNZmbOnOnRFh0d7f79lltu4ZZbbvHYbrFYmD9/vi/KEWnXyk42fj5AU+0irUkni4n4WURY4+cDNNUu0poUAiJ+ljE6g5CgEI+2kKAQMkbrPAHxPYWAiJ8lDk5kxZ0riAyLxISJyLBIVty5QscDpE3oyKtIO5A4OFFv+uIXmgmIiBiYQkBExMAUAiIiBqYQEBExMIWAiIiBKQRERAxMISAiYmAKARFxKbHChv4M/PIq2NDfdVs6PJ0sJiKuN/wfLnFpAtc1jz9PcW3TFc46NM0ERAS+SvW8xjG4bn+V6p96pM0oBEQE7E0sW91Uu3QYCgERgZAmlq1uql06DIWAiMCQDAjwXM6agBBXu3RoCgERgahErD3vo/++AMzfQP99AVh73qeDwgagEBAR18Xut6+mtLYeJ1BaW0/K9tW62L0BKAREpNmL3UvHphAQEV3s3sAUAiKii90bmEJARHSxewNTCIhI+7rY/Q9rGPG6WWsYtQGtHSQiwOmL3RcWFjJo0CD/FHHGGkaA1jBqA5oJiEj78VUq1hN2+pfgOl+hBKwntIaRL2kmICLthvVwKSnfgd3pul1aBynfAZSieYBvaCYgIu1G6vEAdwA0sDtd7W3OINdX8EkIOBwO0tLSiI+PJykpidLSUo/tH3zwAffccw/x8fG88cYbXvURkY6vrLa+Re0+U2LF+v799P9XKQHfOOn/r1Ks79/vlyCwbnmQ/pmBmJ8z0T8zEOuWB1v18X0SAps3b6ampoacnBymTp1KVlaWe1ttbS1z5sxh5cqVZGdnk5OTQ3l5ebN9RMQYIsIiW9TuK9YPHyXlcC2ldbiW0aiDlMO1WD98tG3r2PIgKduWei7nsW1pqwaBT0IgLy+PUaNGATB06FB27drl3lZcXExERARhYWFYLBaGDRvGjh07mu0jIsbQXs5XSD10rPHdUoeOtW0d21c0Xsf2Fa32HD4JAZvNRq9i68gAAAjySURBVGhoqPt2QEAAdXV17m1du3Z1b+vSpQs2m63ZPiJiDO3lfIWyJt56mmr3WR1tsHvMJ98OCg0NpbKy0n3b4XAQGBjY6LbKykq6du3abJ8zFRYWnnddVVVVP6l/R6Kx8KTxOM3fYxETGMO7v37Xo62t67m0Uzj/rqpotL0ta+kXaKasztFoe2vV4ZMQiImJ4cMPP2Ts2LHk5+czYMAA97bo6GhKS0upqKggJCSEHTt2kJycjMlkarLPmX7KSSx+PQmmndFYeNJ4nKaxgLl1fyXl7YnY62vcbSEBFuaO/Wubjk3miP+PlG1LPXYJhZhc7S2pIy8vr8ltPgmBMWPGsG3bNhISEnA6nWRmZpKbm4vdbic+Pp7p06eTnJyM0+nknnvuoXfv3o32ERHxh4bdT6lbUik7WUZEWAQZozPafLdU4uglrjq2r6Cstp6IoAAyfpHibm8NJqfT6Tz33dqHvLw8hg0bdt799QnnNI2FJ43HaRoLTx1hPJp779TJYiIiBqYQEBExMIWAiIiBKQRERAxMISAiYmAX3LeDRESk5Zr6dtAFFQIiItK6tDtIRMTAFAIiIgbW4UNAF6vxVFtby7Rp0xg/fjyxsbFs2bLF3yX53bFjx7jpppsoLi72dyl+t3z5cuLj4/mf//kf3nzzTX+X41e1tbVMnTqVhIQExo8f32H/Pjp8COhiNZ7eeecdwsPDef3113nxxReZNWuWv0vyq9raWtLS0ujUqZO/S/G77du38+WXX7JmzRqys7M5fPiwv0vyq48++oi6ujrWrl3LQw89xMKFC/1dkk90+BDQxWo83XbbbTz66OmrIwUE+OHare3I3LlzSUhI4OKLL/Z3KX73ySefMGDAAB566CEmT57Mf/zHf/i7JL+Kioqivr4eh8OBzWZrdGn7jqBjvqozNHWxmo76D3ouXbp0AVzj8sgjj/DYY4/5uSL/Wb9+Pd27d2fUqFGsWNF6V2q6UJ04cYKDBw+ybNkyDhw4wAMPPMCmTZswmUz+Ls0vQkJC+Pe//83tt9/OiRMnWLZsmb9L8okOPxPw9mI1RnLo0CEmTJjAXXfdxZ133unvcvzmf//3f/n0009JSkqisLCQJ598kvLycn+X5Tfh4eHccMMNWCwWLrvsMoKDgzl+/Li/y/KbVatWccMNN/CPf/yDt99+m+nTp1NdXe3vslpdhw+BmJgYtm7dCtDsxWqM4ujRo0ycOJFp06YRGxvr73L8ymq18tprr5Gdnc2gQYOYO3cuvXr18ndZfjNs2DA+/vhjnE4nR44c4dSpU4SHh/u7LL/p1q2b+1K4YWFh1NXVUV/fepd1bC86/EdiXazG07Jly/j+++9ZsmQJS5a4Lkzx4osv6sCocPPNN/PFF18QGxuL0+kkLS3N0MeMfve73/H0008zfvx4amtrefzxxwkJCfF3Wa1OZwyLiBhYh98dJCIiTVMIiIgYmEJARMTAFAIiIgamEBARMTCFgFwQtm/fzvDhwzl06JC77fnnn2f9+vXn/ZgHDhwgLi6uNco7S319PcnJyYwbN46TJ0+626dPn86dd95JUlIS48aN48EHH2T//v0tfvyKigpyc3Pdj9lwLoxISykE5IIRFBTEU089xYXwreby8nJOnDjBmjVrCAsL89g2bdo0srOzWbNmDRMnTjyvpTt2797NBx980FrlioF1+JPFpOMYMWIEDocDq9XKb3/7W3f7gQMHmDJlCm+88QYAcXFxLFiwgLfeeovS0lJOnDjByZMnGT9+PO+99x4lJSXMnTuXnj17cvz4cSZPnszx48e56aabeOihhzh06BAzZsygurqa4OBgZs2aRX19PQ888ADh4eHceOONTJo0yf3877zzDqtXr8ZisdC/f39mzpzJjBkz2LdvH2lpacycObPJ1zR8+HCCgoIoLS3FYrGc9bx9+vRh/vz57Nq1i8rKSqKjo5kzZw7Lli2jqKiInJwcAHJycnjppZew2Wykp6dzxRVX8Oijj2Kz2aiqqmLatGn84he/8NG/jFzINBOQC0p6ejqrVq1i3759Xt2/U6dOvPzyy9x666189NFHLFu2jJSUFP72t78BYLfbmTdvHmvWrOHjjz+mqKiIuXPnkpSURHZ2NsnJyTz//POA69P9yy+/7BEAJ06cYNGiRaxevZo1a9bQtWtXcnJyePbZZ7n88subDYAGPXr04MSJE40+r81mo1u3brzyyiusXbuW/Px8jhw5wuTJkxkxYgTx8fEAXHXVVbz66qv89re/Zf369ZSVlXH06FGWLVvG/PnzqaqqauFIi1FoJiAXlIsuuoinn36a6dOnExMT0+h9ztxddOWVVwLQtWtXLr/8csC1DkzDQmADBw50rw8zePBgSkpK2LNnD8uXL+ell17C6XQSFBQEQN++fbFYLB7PtX//fi6//HL3SrXXXXcdn3zySYuWYT548CCXXHJJo8/bsIjblClTCAkJwW63U1tbe9ZjXHXVVQD07NmTqqoqfv7zn5OYmMiUKVOoq6sjKSnJ63rEWBQCcsG55ZZbeP/993nrrbeYNm0awcHBHDt2jPr6eiorKzlw4ID7vudaBrm4uJjKykqCg4PZuXMn8fHxXHbZZUycOJGYmBiKi4v54osvADCbz5449+3bl+LiYux2OyEhIXz++edERUV5/Vq2bdtGp06duOSSSxp93q1bt3Lo0CEWLlzI8ePHef/993E6nZjNZhwOR5Ovc/fu3VRWVrJixQq+++47EhISuPnmm72uS4xDISAXpNTUVD777DMAevXqxciRI4mNjSUiIoLIyEivHycsLIzHH3+c48ePM3bsWC6//HKefPJJ0tPTqa6upqqqitTU1Cb7d+/enT/84Q9MmDABs9lMREQETzzxRLNLUs+bN48XX3wRs9lMly5d3Fesaux5+/bty5IlS4iLi8NisdCvXz++++47IiIi2LNnD6tWrWr0Ofr378/ixYvZsGEDQUFBPPLII16PiRiLFpATETEwHRgWETEwhYCIiIEpBEREDEwhICJiYAoBEREDUwiIiBiYQkBExMAUAiIiBvb/ABtSd5nh2yKGAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from scipy.stats import poisson\n",
    "plt.style.use('seaborn-whitegrid')\n",
    "\n",
    "\n",
    "## proportion of different number of deaths\n",
    "plt.plot(np.arange(freq.size), freq/np.sum(freq), 'o', \n",
    "         color='orange',label = 'Original Data')\n",
    "plt.plot(np.arange(freq.size), \n",
    "         poisson.pmf(np.arange(freq.size),death_mle),'o',\n",
    "         color='green',label='Poisson Distribution')\n",
    "\n",
    "plt.legend()\n",
    "plt.xlabel(\"Number of Deaths\")\n",
    "plt.ylabel(\"Proportion\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Problem 2\n",
    "Fit the  𝛼 ,  𝜇1  and  𝜇2  parameters for the mixture of two Poissons model proposed in Question 11 using a gradient ascent algorithm. You must write your own fitting code (you cannot use scipy optimize or R optim). Discuss how you select starting parameters and monitor convergence. Include a plot to illustrate the convergence for different starting values of the parameters."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In problem 1 we used standard notation for the data i.e. $x_1$ was the number of deaths on day 1. <br>\n",
    "Now we let $n_i =$ number of days with $i$ deaths. \n",
    "So$$n_i = \\sum_j \\mathbb{1}_{X_j=i}$$\n",
    "\n",
    "The likelihood is$$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",
    "\n",
    "\n",
    "So the log likelihood is$$ \\log \\, f = \\sum_{i=0}^9 n_i \\log(g_i) $$\n",
    "\n",
    "The gradient is$$ \\frac{\\partial log f}{\\partial \\alpha} = \\sum_{i=0}^9 \\frac{n_i}{g_i} (e^{-\\mu_1} \\frac{\\mu_1^i}{i!} - e^{-\\mu_2}\\frac{\\mu_2^i} {i!})$$and $$ \\frac{\\partial log f}{\\partial \\mu_j} = \\sum_{i=0}^9 \\frac{n_i}{g_i} \\frac{\\partial g_i}{\\partial \\mu_j}$$ for $j=1,2$ and where$$ \\frac{\\partial log f}{\\partial \\mu_1} = -\\frac{n_0}{g_0} (\\alpha e^{-\\mu_1}) + \\sum_{i=1}^9 \\frac{n_i}{g_i} \\frac{\\alpha e^{-\\mu_1} (i- \\mu_1) \\mu_1^{i-1}}{i!}$$$$ \\frac{\\partial log f}{\\partial \\mu_2} = -\\frac{n_0}{g_0} ((1-\\alpha) e^{-\\mu_2}) + \\sum_{i=1}^9 \\frac{n_i}{g_i} \\frac{(1-\\alpha) e^{-\\mu_2} (i- \\mu_2) \\mu_2^{i-1}}{i!}$$\n",
    "\n",
    "\n",
    "\n",
    "<br>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": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "#defining g function\n",
    "def g(alpha,mu1,mu2,i):\n",
    "    i_fact = scipy.special.factorial(i)\n",
    "    return (alpha*np.exp(-mu1)*np.power(mu1,i)/i_fact) + ((1-alpha)*np.exp(-mu2)*np.power(mu2,i)/i_fact)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Definig LOG-LIKELIHOOD function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "def log_likelihood(params,n):\n",
    "    like=g(params[0],params[1],params[2],n.size)\n",
    "    if np.min(like) < 0:\n",
    "        return -np.inf \n",
    "    else:\n",
    "        z =0\n",
    "        for i in range(0,10):\n",
    "            p = n[i]*np.log(g(params[0],params[1],params[2],i))\n",
    "            z = z+p\n",
    "        return z "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Calculating Loglikehoood - Gradient Ascent derivates\n",
    "\n",
    "Alpha-gradient log"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "def lograd_alpha(params,n):\n",
    "    alpha = params[0]\n",
    "    mu1 = params[1]\n",
    "    mu2 = params[2]\n",
    "    z = 0\n",
    "    for i in range(0,10):\n",
    "        i_fact = scipy.special.factorial(i)\n",
    "        p = (n[i]/(g(alpha,mu1,mu2,i)*i_fact))*((np.exp(-mu1)*np.power(mu1,i)) - (np.exp(-mu2)*np.power(mu2,i)))\n",
    "        z = z+p\n",
    "    return z\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Defining MU gradients"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "#mu1 gradient\n",
    "def lograd_mu1(params,n):\n",
    "    alpha = params[0]\n",
    "    mu1 = params[1]\n",
    "    mu2 = params[2]\n",
    "    \n",
    "    z = 0\n",
    "    const = -n[0]*alpha*np.exp(-mu1)/(g(alpha,mu1,mu2,0))\n",
    "    for i in range(1,10):\n",
    "        i_fact = scipy.special.factorial(i)\n",
    "        p = (n[i]/g(alpha,mu1,mu2,i))*(alpha/i_fact)*(np.exp(-mu1)*np.power(mu1,i-1)*(i-mu1))\n",
    "        z=z+p\n",
    "    return (z+const)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "#mu2 gradient\n",
    "def lograd_mu2(params,n):\n",
    "    alpha = params[0]\n",
    "    mu1 = params[1]\n",
    "    mu2 = params[2]\n",
    "    z = 0\n",
    "    const = -n[0]*(1-alpha)*np.exp(-mu2)/(g(alpha,mu1,mu2,0))\n",
    "    for i in range(1,10):\n",
    "        i_fact = scipy.special.factorial(i)\n",
    "        p = (n[i]/g(alpha,mu1,mu2,i))*((1-alpha)/i_fact)*(np.exp(-mu2)*np.power(mu2,i-1)*(i-mu2))\n",
    "        z=z+p\n",
    "    return (z+const)  \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Functioning the gradient derivates"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "def gradient_ascent(param_initial):\n",
    "    #running for 1000 iterations for the algorithm to converge\n",
    "    #(as we dont know where will it converge)\n",
    "    Niter = 1000\n",
    "    #creating matrix of no. of iterations size\n",
    "    param = np.zeros((Niter+1)*3).reshape(Niter+1,3)\n",
    "    param[0,:] = param_initial\n",
    "\n",
    "    for ii in np.arange(Niter):\n",
    "        lgrad_alpha = lograd_alpha(param[ii,:],freq)\n",
    "        lgrad_mu1 = lograd_mu1(param[ii,:],freq)\n",
    "        lgrad_mu2 = lograd_mu2(param[ii,:],freq)\n",
    "        lgrad = np.array((lgrad_alpha, lgrad_mu1, lgrad_mu2))\n",
    "        half = True \n",
    "         ## step half until new point increases likelihood    \n",
    "        while half:\n",
    "            param_temp = param[ii,:] + lgrad\n",
    "            newlik = log_likelihood(param_temp,freq)\n",
    "            ## for temp far from maximum, logit_loglik returns\n",
    "            ## Nan values. we replace these with -Inf log likelihoods\n",
    "            if np.isnan(newlik): \n",
    "                newlik = -np.inf\n",
    "            if (newlik < log_likelihood(param[ii,:],freq)):\n",
    "                lgrad = lgrad / 2\n",
    "            else:\n",
    "                half = False\n",
    "        param[ii+1,:] = param_temp\n",
    "    return param\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We choose starting values of parameter as follows(taken from textbook):\n",
    "Since in 1st question when we use single poisson distribution we get mean value of 2.156. \n",
    "Therefore, for the mixture of poisson distribution one value should be below and one value should be above this value. \n",
    "Hence here we use mu1=1 and mu2=2.5 \n",
    "For alpha, we choose initial value of 0.3 since we believe that the contribution of parameter with high mean will be more as\n",
    "observed from the graph in the 1st question."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\jsruj\\Anaconda3\\envs\\Stat 654\\lib\\site-packages\\ipykernel_launcher.py:8: RuntimeWarning: invalid value encountered in log\n",
      "  \n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array([0.35988489, 1.25609417, 2.66340372])"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "param_initial_1 = np.array((0.3,1.0,2.5))\n",
    "param1 = gradient_ascent(param_initial_1)\n",
    "param1[1000,:]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Taking other alpha, mu1,mu2 values and checking the converges "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\jsruj\\Anaconda3\\envs\\Stat 654\\lib\\site-packages\\ipykernel_launcher.py:8: RuntimeWarning: invalid value encountered in log\n",
      "  \n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array([0.35988588, 1.25609595, 2.66340494])"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "param_initial_2 = np.array((0.5,0.5,4.1))\n",
    "param2 = gradient_ascent(param_initial_2)\n",
    "param2[1000,:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\jsruj\\Anaconda3\\envs\\Stat 654\\lib\\site-packages\\ipykernel_launcher.py:8: RuntimeWarning: divide by zero encountered in log\n",
      "  \n",
      "C:\\Users\\jsruj\\Anaconda3\\envs\\Stat 654\\lib\\site-packages\\ipykernel_launcher.py:8: RuntimeWarning: invalid value encountered in log\n",
      "  \n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array([0.35988518, 1.25609471, 2.66340409])"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "param_initial_3 = np.array((0.8,0.1,1.0))\n",
    "param3 = gradient_ascent(param_initial_3)\n",
    "param3[1000,:]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot Parameter Paths\n",
    "\n",
    "Plotting the $\\mu_1$, $\\mu_2$ parameter paths. We could also plot these paths for $\\alpha$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "plt.style.use('ggplot')\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAECCAYAAAD9z2x7AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3deVyU9fr/8deszAwgiwSIAWnmknDMJTI1S0P9euy0ly2GbafUtEhNSY/7MT1HU1uOx1OdfkfLOn5za/X0taxc8qRmWpaZucuSIosgszDL74+REYRBlllu4Ho+Hjxg7hnueQPDdd/zuT/3datcLpcLIYQQzZ462AGEEEIEhhR8IYRoIaTgCyFECyEFXwghWggp+EII0UJIwRdCiBZCCr4QQrQQ2mAHuJScnJxgRyAmJob8/Pxgx6hGibmUmAkkV30pMZcSM4EycyUkJNS4XPbwhRCihZCCL4QQLYQUfCGEaCGk4AshRAshBV8IIVoIKfgt2Nq1RtLSYrn88jakpcWydq0x2JGEEH6k+GmZwj/WrjUyaVIEZrN7m5+drWXSpAgA7rzTHMxoQgg/kYJfT3v27GHWrFkkJyejUqmwWq2kp6dz5513ev2eDz/8kKFDh7Jv3z4++OADpk+fXutzWCwWJk6cyKRJk0hKSvL1jwDA/PnhnmJfwWxWM2lSBAcOaGnb1kHbtg4uv9z9OSxMLpsgRFMnBb8Bunfv7inaNpuNjIwMBg8eTFhYWI2PX7lyJUOGDKnTug8cOMCiRYs4ffq0z/LWJCdHU+Nys1nFsmVh2O2qKssjIpyVNgJ2z9cVG4XLLnOilgFCIRStSRd8Y957mHJX+XSdZW2GY46/p86PN5vNaDQaNBoNe/bsYcWKFYB7L/3555/n+++/p6CggNmzZ3P33Xdz8uRJJk+eTGFhIX369OHhhx+usj6bzcacOXN44YUXfPljVZOQ4CA7u/qfv21bB9u3n+LUKTXZ2ZrzH1qyszWcPOn++OYbPWfPVq3uOp2LhAQH7dqpueyySM87g7ZtL2wcDAa//khCiEto0gU/WL777jsyMzNRq9VoNBrGjRuH0Wjk6NGjTJkyhZiYGN5++22++uorRowYwVtvvcX06dP56aefPAXd6XQyfPjwagU/NTU1ID9DVlZJlTF8AKPRSVZWCRoNtGnjpE0bJ716ldf4/WfPqiptEC58/Pabga1bQ/jtNzVOZ9V3CTExjmrvDCp/jopyolLV+HRCCB/wa8EvLi4mKyuLP/3pT7Rt29azfNeuXaxZswa1Ws2AAQNIT09v0PrN8ffUa2/cVyoP6VQWExPDK6+8gtFoJD8/n5SUlGqPadeuHXq9HgCNpuZhlUCoODA7f344OTkaEhIcZGWV1PmAbatWLlq1stOli73K8oq+IuXlkJfnfkdw8Ubhl1+0bNoUgsVS9V2C0eissjFISKi6UYiPd6DT+ebnF6Il8lvBt9vtvPbaa57iVnn58uXLmTdvHgaDgWnTptGrVy8iIyP9FSVgFi5cyDvvvIPJZGLevHlUXB9erVbjdDoBUCloF/bOO81+m5Gj00FiooPEREeN97tcUFio9gwVVf6ck6Phxx915OdX3SCq1S7i4pyeoaKaNgrh4XJwWQhv/Fbw33rrLQYNGsT69eurLM/OziY+Pt5zgLNTp07s37+f66+/3l9RAmbQoEGMGTOGsLAwoqOjOXPmDOAepsnKymLkyJFBTqgcKhVERzuJjnaSmlrzsJHZ7D64fPFxhOxsDXv36vnkEw3l5VU3oK1aOWsYNrLTtauK0FA1cXFycFm0XCpXxW6oD3355ZecOXOGu+66i5kzZ/LHP/7RM6Tz888/s2HDBp599lkAVq1aRUxMDDfffHON67LZbL6OV29arRa73X7pBwaYEnMFMpPTCb/9BsePqzhxQsWxY3DihOr8bffXhYVVNwg6nYu2bSE52UVioovEREhKcn/t/gwmU0DiA8r8G4IycykxEygz18UjKxX8sof/xRdfAPDDDz9w9OhRXn31VSZPnkxkZCRGoxGLxeJ5rNlsJjQ01Ou6lNBnWon9rkGZuQKdSaeDK690f9SkpERFTo6Gs2ej2L//nOc4wsmTGjZt0pCXp6l2cLl166rHES4+wBwd7buDy0r8G4IycykxEygzl7d++H4p+LNmzfJ8XbGHXzFG37ZtW3JzcyktLcVgMLB//35uvfVWf8QQgvBwF5062YmJcXHttWXV7rfb3QeXLz6WkJOj4eBBLV98EVLtBLUuXcr57DP/nichhD8EbFrm1q1bsVgspKenk5GRwdy5c3E6nQwYMIDo6OhAxRCiCq0WLr/cved+3XXV73e5oKhI5TmOsGaNkY8/NnLunIrQUDlALJoWvxf8mTNnAlSZltmrVy969erl76cWotFUKoiKchEVZSclxU5xsYqPPzZy+rSa0NCaZyAJoVQyX0GIeoiNdU+vPX06eOdQCNFQUvCFqIfLLnPv1Z8+Lf86oulp1q0Vwo4vxRbeDVtUX88yfeE29CV7KU0a06B15uXl8dhjj3HVVVd5lnXv3p2+ffuybdu2es+1z8zMZPz48TV2xfzb3/5GYmKiHNRWkMsuc+/hnzolBV80Pc264NvCuxH10ygKr16GLaov+sJtntuNkZyczJIlS6ot79ChQ6PWW6GoqIh58+Zx8uRJhg8f7pN1Ct9o3dqJWu1iyZJwPv3UQGKi+4Cv+7OdxEQHsbFycpdQpuZd8KP6Unj1MqJ+GkVZQgamnBWe4u9re/bs8fS6HzFiBCkpKZw4cYKoqChmzZqFxWJhwYIFnDt3juLiYoYNG8Ztt91W47rMZjMjR45kx44dPs8pGkejgblzi/n66xBOntTw6afVW0Do9S7atnWQmGiXDYJQlGZd8MFd9MsSMgg/toSS5EyfFPtjx46RmZnpuT116tQq9+fm5rJo0SJiY2MZO3YsBw4cQKvVMnDgQPr3709+fj6ZmZleC36bNm1o06aNFHyFysgoIyPjwpz+sjL3tM0TJ9wfJ09qOHFCW6cNQocOGi67LEw2CCIgmn3B1xduw5SzgpLkTEw5K7BG9ml00a9pSCc7O9vzdUREBLGxsQDExsZis9mIjY1l9erVbNmyBZPJpLhTsUXDmUwurrrKzlVX1fw3rW2D8NFHak6dalXl8SEhFzYIFe8OKm8Q5GIzoqGadcGvPGZvi+qLNbJPlduBtGrVKrp27cptt93Gd999xzfffBPQ5xfBU9sGISYmhuPHz5zfCFRsELSeDcO+fTrOnKn6DqGmDUK3buXccINVricgatW8C37J3irFvWJMX1+yN+AFv0+fPixevJjPPvuMVq1aodFoFNEYTgSfyeSiY0c7HTt6f4dQlw3CTTdZmDOnmPbt5YQwUTO/dMv0pZycnGBHUGRzJFBmLiVmguadq6xMxbvvmliwIByrVcXo0aWMG1eK0djwf20l/r6UmAmUmctb8zQZCRSiiTOZXDz22Dk2bz7FH/5g5qWXwrnppsv49FMDyt6dE4EmBV+IZiI21snLLxexZk0+YWEuHn00moyMaI4elTYQwk0KvhDNTO/eNv7zn9PMmFHMjh16Bg6MZeHCcMz+uZqlaEKa9UFbIVoqnQ6eeOIct95q5s9/bsXixeGsWWPk9tvNnktLRkVd+BwV5SQ83CWzfJo5KfhCNGPx8U5efbWIBx4oY+bMCF59NazaFb4qaLUuz0YgNlZDeHiUZ2NQ8VF5AxEd7SQiwiXnBDQhUvCFaAH69LHxf/93GqcTiotVFBaqKSxUU1BQ9XPFR0mJhkOHtJ777PaaNxJqtYvIyIs3Bi46diznj388h0YOHyhKsy/4xt/WEn54PhprDo6QBEraZ2GOu7PB68vLy2P27NksXbq0xvsLCgpYvnw5zz77LHv37iUsLIwrr7yS6dOnM3v27Hqvs6ioiLFjx/Lmm296vTCxEHWlVldc0MUBeJ+vX3mqocsFpaWqGjcOFy87cULL99+rWbXKxPHjWubOLZZhIgVp1gXf+NtaIg5MQu10H63SWrOJODAJoFFFvzbR0dE8++yzAGzYsIGBAwdy5ZVXei32tdmxYwevv/46hYWFvo4pRJ2pVO5rA4eHO0hOrttJXX/+cyv+/vcwEhIcjB1b6ueEoq6adcEPPzzfU+wrqJ1mwg/P90nBz8zMpEOHDhw5coSysjJmzJgBwOzZs3nmmWfYsWMHBw8eJDk5mdGjR7N27Vr27NnDihUrALBYLDz//PPodLoa169Wq1m4cCFPPvlko7MKEUhTppwlL0/NvHmtiItzcM89MkVICZp1wddYaz5L19vyhujcuTNjx47ljTfeYNOmTQwcOBCATp06kZaWxsCBA4mLi/M8/ujRo0yZMoWYmBjefvttvvrqK9LT02tct1z3VzRVajUsWlTE6dMaJk6MpKBAzRVXOIiKcnrG/CMjnXjZ1xF+0qwLviMkAa01u8blvlJx5avY2FgKCgou+fiYmBheeeUVjEYj+fn5pKSk+CyLEEqi18MbbxRw772tmT07osbHhIVV3gC4qm0QKn/u2hWMxgD/EM1Msy74Je2zqozhAzjVRkraZ/nsOVS1HJFSqVQ4nc4qyxYuXMg777yDyWRi3rx5KLyVkRCNEh7u4sMP88nJ0VBU5D6wW1R0YZbQhWXuz9nZOoqKVBQVqWucPjp2bDiTJpXI7J8GatYFv2Kc3pezdOqjS5cuvP7667Rp08azbNCgQYwZM4awsDCio6M5c+ZMQLIIESxaLSQlOUhKqnsXT6cTSkpUno1BUZGazz+P5NVXw9m7V8/SpYVERzsvvSJRhXTLrAMldsMDZeZSYiaQXPWlxFwxMTG8+uo5pkyJJCbGwWOPnaNNGwdt2jiJj3cQF+cgGDOXlfi78tYt0297+E6nk2XLlpGbm4tarWb06NHEx8d77v/oo4/YtGkTrVq5r/bzxBNPeA0phBAA991npksXO6NHR9V4XKB1awfx8e4NwIWPqrejolpuCwm/Ffxdu3YBMGfOHH788UdWrFjBpEmTPPcfPnyYsWPH0r59e39FEEI0Q926lbNt2ymKilTk5WkqfajJzb1we8+e6lcLAzAYXMTFVd8gtG3rYNAgS1DeJQSK3wp+WloaPXv2BOD06dNERFTdGh85coR169ZRVFREjx49uOOOO/wVRQjRzKhUFWcM2+nSxfv1oa1WOH1aQ26uutrGIS9Pw969ej79VIPF4t7lf/HFQu67r/meM+D3MfxXX32VnTt3Mn78eLp16+ZZ/t577zFkyBBMJhMLFixg8ODBng1EZUq4DKBWq1XkRceVmEuJmUBy1ZcSc/krU2Eh9OihIyQEduwop1WrS39PIHI1hrc2LAE5aFtUVMSUKVNYtGgRBoMBl8uF2WzGZDIB8Omnn1JSUsLdd99d7XvloK13SsylxEwguepLibn8lWnRojBefLEVISEuYmIc1c4HiIpy8j//Y6Fbt/KA5mqMgB+03bx5M2fOnOGOO+5Ar9ejUqlQn++jajabmTBhAosXLyYkJIR9+/YxYMAAf0URQgiv7r3XjEZDlfMBCgvV5OZqPctefjmc++47R1ZWCZdd1nSng/ptD99isbB06VKKi4ux2+3cfvvtWK1WLBYL6enpbN68mQ0bNqDVaklNTeXee++tcT1K28Pfs2cPs2bNIjk5GZVKhdVqJT09nTvv9D63/8MPP2To0KHs27ePDz74gOnTp3t97Oeff87q1avRaDS0b9+ezMxMz4aytlxKocRMILnqS4m5gpWppETFSy+F88YboRgMLp59toQbb7RiMrkwGl0kJkZTVpavqOsCBHwP32AwMH78eK/39+/fn/79+/vr6f2qe/funqJts9nIyMhg8ODBhIWF1fj4lStXMmTIkEuu12q18uabb/LPf/4Tg8HAnDlz2L59O3379vVpfiFE3YWHu/jTn85y333nmDkzwkubiAQMBhf9+1tYtqyQkJCAx6yTJn2m7XvvGVm1yuTTdQ4fXlavzn5msxmNRoNGo6mxE+b3339PQUEBs2fP5u677+bkyZNMnjyZwsJC+vTpw8MPP+xZl06n45VXXsFgMADgcDikB74QCtGhg4O33ipg924dOTkazGYVZWUq1Oow8vPLOHVKw1tvhTJhgounny5Fr3cREuIiJASMRqci+gA16YIfLN99951nqEWj0TBu3DiMRmONnTBHjBjBW2+9xfTp0/npp5+w2WzMmTMHp9PJ8OHDqxR8tVpNdHQ0AGvXrsVsNkvHTCEURKWCnj3L6dnzwgHcmBgT+fnunv9t2jj4619bsW5d1R1RjcbFmjVnuPba4M46bNIF/557zEHps115SKeyunTCbNeunWevXVNDByin08k//vEPTp48yaxZs2ptziaEUJanny7luutsnDqlxmpVYbWqsNlULF4cRmZmJL/7XTkmk5PQUBejRpWSkBDYA8BNuuArjbdOmGq12tM181IFfNGiReh0OubMmeP1YK0QQplUKujdu/pevMnk5J//DOOHH3RYLCp++03N/v060tMtXH11OTfcEJg9fyn4PuStE2ZqaipZWVmMHDmy1u//5Zdf+OSTT0hNTfUc8L7rrru44YYb/J5dCOE/991nrnIG7/Tprfh//y+Ur792H90dM6aEW2+1kJpa81x/X5FumXWgxClqoMxcSswEkqu+lJhLiZmg4bmcTvfc/0ceiebbb/UYjU7+85982re3N3qKp7dpmTJmIIQQQaBWQ+vWTj74IJ/t238D4MYbY0lMTODsWf8cu5OCL4QQQZaU5GD9+gvvEnr2jKvl0Q0nBV8IIRQgJcXOF1+cAqCsTI0/Btul4AshhEJ8/LH7pMvkZLtfLtIiBV8IIRQgJ0fNwoXu3sxr1/rn4LQUfCGECLKff9YyYUIkACtWnCE+3j8nZDXrefhLl4bRrZuNvn0vnNSwbZuevXv1jBlT2qB15uXl8dhjj3HVVVd5lnXv3p2+ffuybdu2S861v1hmZibjx48nKSnJs+zXX3/l5ZdfRq1Wo9PpeP755z0tF4QQzcvOnTpuv/0yAEaNKmXgQKvfnqtZF/xu3WyMGhXFsmWF9O1rY9s2ved2YyQnJ7NkyZJqyzt06NCo9VZ45ZVXePrpp+nQoQMffPAB7777Lk899ZRP1i2EUAaXCz780MC8ea1o29bOnDlnGTLE4tfnbNYFv29fG8uWFTJqVBQZGWWsWGHyFH9f27Nnj6fX/YgRI0hJSeHEiRNERUUxa9YsLBYLCxYs4Ny5cxQXFzNs2DBuu+22Gtc1ffp0WrduDUjHTCGao6+/1nPvva1xuVRERDh5882CGlsy+FqzLvjgLvoZGWUsWRJOZmaJT4r9sWPHyMzM9NyeOnVqlftzc3NZtGgRsbGxjB07lgMHDqDVahk4cCD9+/cnPz+fzMxMrwW/otjv27eP9evX89JLLzU6sxAi+MxmePnlcF5+ORyAO+4oY/HiInS6wDx/sy/427bpWbHCRGZmCStWmOjTx9rool/TkE52drbn64iICGJjYwGIjY3FZrMRGxvL6tWr2bJlCyaT6ZIXPd60aRMrV65k3rx5REZGNiqvECL4iotVXH11GwDUahdr1+Zz7bX+7Z1zsWY9S6fymP1zz5V4hne2bQv8EMmqVavo2rUrU6dO5aabbqr1sRs3bmT9+vUsXrzYa08MIUTT8a9/mTzFHuDAgbyAF3to5nv4e/fqq4zZV4zp792r98s4fm369OnD4sWL+eyzz2jVqhUajQabrXoGh8PBK6+8QmxsrKfnfrdu3XjkkUcCmlcI0XiFhSpSUi4U+jvvLGPhwqKgXQJRumXWQXPr0udPSswEkqu+lJhLiZnAe66pUyP4179CPbcPHszFZApMuQ34RcyFEKIlOnpUw0MPtebwYXd5/fe/8+nXz+aXVgn11azH8IUQIlDsdli1ykjfvnEcPqzl9tvL2L07jxtuUEaxB9nDF0KIRnG5YOHCcFasMHHmjIZOncp56qlS7ror8NfbvhQp+EII0UAnTmho21YPuGf+vflmAYMHWxSzR38xvxV8p9PJsmXLyM3NRa1WM3r0aOLj4z3379q1izVr1qBWqxkwYADp6en+iiKEED7lcMA997Tmm28uTLf59NNTpKTUfn5NsPltDH/Xrl0AzJkzh3vvvZcVK1Z47rPb7SxfvpypU6cya9YsPv/8c4qKivwVRQghfMJmg3ffNXHjjbGeYn/TTU6OHctRfLEHPxb8tLQ0nnzySQBOnz5NRESE577s7Gzi4+MJCwtDq9XSqVMn9u/f75cca9caSUuL5fLL25CWFsvatcZGrS8vL48xY8Z4vb+goIDFixcDsHfvXg4dOgTgmVNfn3W+9957jB49mtGjR7N8+fJG5RZCNM5XX4XQrl0CEydGEhbm5LXXCjh+PIdPP7WjbSKD436NqdFoePXVV9m5cyfjx4/3LDebzZhMJs9to9FIWVlZjeuIiYlp8PO/+66ayZM1lJW5B9Sys7VMnhxJeHg4999f937TWq3Wk8NqtaLT6bzmiomJYe7cuYC7PcLQoUOJiYlh6dKlXtdf0zpPnDjBl19+yb///W9UKhUjR47klltuoVOnTjXmUgolZgLJVV9KzBWMTC4XLFigZu5cDRaLu4507+5k+3ZQqcKAMEX+rrzx+3Zp7NixFBUVMWXKFBYtWoTBYMBoNGKxXGgDajabCQ0NrfH7G3OixdSpsZ5iX6GsTMXUqTBoUN3XW/nEisLCQsrLyz0N0Dp06MCRI0coKytjxowZAMyePZtnnnmGLVu2sG/fPqKiohg9ejRr165lz549nuEti8XC888/j06n86yzgk6n44UXXqCwsNDz2LKysiqPUeKJKErMBJKrvpSYK9CZ9uzRMWzYZZ7bEyac5YknzhEW5uLMmeDlqgtvJ175bUhn8+bNrFu3DgC9Xo9KpUKtdj9d27Ztyc3NpbS0FLvdzv79++nYsaPPM+TkaOq1vCE6d+7Miy++SM+ePdm0aZNneadOnTzDWnFxF65Af/ToUc/Gr0+fPnz11Vc1rler1RIREYHL5eLvf/87HTp0IDEx0We5hRDVuVzuoZsRI6KrFPu//a2Q8eNLCQtTdGOCS/LbHn5aWhpLly5lxowZ2O12Hn74YXbs2IHFYiE9PZ2MjAzmzp2L0+lkwIABfrmiU0KCg+zs6j9iQoLDZ89RceWr2NhYCgoKLvn4mJgYXnnlFYxGI/n5+aSkpHh9rM1m4y9/+Qsmk6lKO2YhhG9ZLLBunYnXXw/lwAEdsbEOnnvuLBkZZURH++dyg8Hgt4JvMBiqjNtfrFevXvTq1ctfTw9AVlYJkyZFYDZfeCNjNDrJyirx2XOoaplwq1KpcDqrvlgWLlzIO++8g8lkYt68eXhrZeRyuZg6dSo9evTg/vvv91leIcQF5eUwfnwka9e6jylefXU5ixcXcttt5qA1OPOnJnJsuWHuvNN9ptv8+eHk5GhISHCQlVXiWe5vXbp04fXXX6dNmwvd8gYNGsSYMWMICwsjOjqaM5UHAyvZunUre/fupby8nG+++QaAP/7xj3Tt2jUg2YVozi7uYgnw6KOlzJ59VrEnTfmCdMusAyUelAFl5lJiJpBc9aXEXL7I9O9/G5kwIarKsn/96wwDB1rRNPDQnhJ/V9ItUwjRIpnNkJISj8VSdY7K4sWF3HOPuVnv0V9MCr4QollatiyUOXMiqi3/3//ND/gFkJRCCr4QotnIzlbz4YdG3n/fyPffV72U6fff59G6dfOZcdMQUvCFEE3ab7+p+fhjIx98YGDnTvfUmtRUG8OGmbn77jIGD7YGOaFySMEXQjQ5Z86o+fhjAx98YOS//9Xjcqno0qWcSZPOcuutZtq18925Ns2JFHwhRJNgt8M775h47rlIz7IrrywnM7OUW28107Gj8rtVBpsU/Hras2cPs2bNIjk5GZVKhdVqJT09nTvvvNPr93z44YcMHTqUffv28cEHH9TaOfOrr77i3XffRaVSccsttzBs2DB//BhCNAk2G6xebTp/NSkdOTnuYh8V5WDVqjNcfbW9Rc2yaSwp+A3QvXt3T9G22WxkZGQwePBgwsLCanz8ypUrGTJkyCXX63A4eP3111m2bBlGo5FHHnmEfv36VWktLURz9+mnBv7yl3Cio518953OM51y6FAnc+YUcvPNFnS6IIdsopp0wTe+9x6mVat8us6y4cMx33NPnR9vNpvRaDRoNJoaO2F+//33FBQUMHv2bO6++25OnjzJ5MmTKSwspE+fPjz88MOedWk0GpYvX45Go6GwsBCXy4XR2Lj+/UIoXXk55OZq+OQTA19+aWDLFveB19/9zsaIEWV07mynU6dyBg+OID/fcom1ido06YIfLN999x2ZmZmo1Wo0Gg3jxo3DaDR6OmHGxMTw9ttv89VXXzFixAjeeustpk+fzk8//YTNZmPOnDk4nU6GDx9epeCDu+hv3ryZl156id69e6Np6Ol/QijYb7+p6dHDfclTtdqF03lhXCYzs4QHHiijbVs58OprTbrgm++5p157475SeUinsrp0wmzXrh16vXt+sLdi3r9/f/r168f8+fP5v//7P4YOHerbH0CIALNaYfduPd98o2fr1hC2b7/Qmezpp0tJTLRjsagYOtRCXFzLnivvT0264CuNt06YarXa0zWztu6a586dY8qUKSxYsAC9Xo/RaKz18UIEw093T8PQL4X2mcM9yw4vWYVl6z6uXj0HcA/TfPFFCDt36tm3T8eOHXrPWHznzuWMH1/CsGFmOneWmTWBJAXfh7x1wkxNTSUrK4uRI0fW+v2hoaGkp6fzzDPPoNVqad++PYMGDQpEdCHqzNAvhc4L5/Iz0D5zOIeXrKLzwrlsuPcFNr4Uxn//q2fnTr2nLXnnzuU8+GAZ/fpZue46GxERiu7X2KzV2i2zpKSE5cuXc+bMGa6//noGDx7suW/hwoVMnDjR7wGlW6Z3SsylxEwguerrUrmOv7CMq/+2hM+jhtKvcBP38r98yQAAunQpp3dvK5072xk40EJCgm+GaJrq7yoYGtQt87XXXiMuLo6ePXuyevVqjh8/zuOPPw7AqVOnfJ9SCKE8jnJCfvgc49h8OfYAAB3OSURBVPb/oPv2W7T7TpJwwt187I7C/2Wx8TlMA3qzcGARQ4ZYmtUVopqbWgt+Xl4eEyZMANzDEjNnzuSjjz7illtuCUg4IUTgaU/sQ/3RZ7TeugndD7+iOlCM6vw1g1yhKpxXR5Mdn0TEzkPs6vo/PPTTm1zbtQ3t7x9e+4pF0NVa8J1OJw6HA41GQ1hYGBMnTmTatGkkJibKwUQhmjqXC92R7wj5/kv0u7ej+/Eg6p2nUZ2fDanRgCvJgGNgAtZr+2O+/hZsXfpz+JXVdF44lx+em0mXzOH8fH4Mv2JMXyhXrQW/e/fuvPDCCzz00ENcccUVxMfH88wzz/Diiy8GKp8QwhecTvQ/bEF/YDuGvdvR7d0Nh5yozrrvdqmABB2uKyKw3Ho9upv+wKkuN0JoVLVVWbbu4+eJUz3FvX3mcH4+vxwp+IpWa8F/8MEH+eKLLygvL/csS0lJYdq0abz77rt+DyeEqD918WlC9n2Ffs829D9+j+bXHFSHSlBZzs/PUAFtwZEWjyM8mrLbhmO59hackfGedcTExICXA5EVUy8ra585XIp9E1BrwVepVAwcOLDa8vbt2zN16lS/hRJC1IHVTMi+LYT88DW6A9+j/eUo6iNnUP12YW67ywCuK0Kx/8+VOKIvx5z+e6zdB+FsFRvE4CJY6jQPf8KECTWO2S9cuNDngXwpbOlSbN26Yevb17NMv20b+r17KR0zpkHrzMvL47HHHuOqq67yLOvevTt9+/Zl27Ztl5xrf7HMzEzGjx9PUlKSZ9nRo0d58cUXcblcXHnllTz99NPSYqElO1dMyE/b0P+yA/3hn9AcPIbmwGnItXrG211qIF6L86oYym9pT/nVv8N6TX9sHfuCWk63EW51eiU89thjnq/tdjvbtm0jLi7Ob6F8xdatG1GjRlG4bBm2vn3Rb9vmud0YycnJLFmypNryDh06NGq9Fd544w0ef/xxunXrxvz58/n666+54YYbfLJuoVyq0lL0e7YScvi/6PbtRXvoKOqTBZBr9xR2AFekGmfHSBw3dKG869XYuvTGkjIAV1h08MKLJqFOBf/qq6+ucjs1NZU//elPXnvA2+12/v73v3P69GnKy8u566676NWrl+f+jz76iE2bNtGqVSsAnnjiCa8nCjSGrW9fCpctI2rUKMoyMjCtWOEp/r62Z88eT6/7ESNGkJKSwokTJ4iKimLWrFlYLBYWLFjAuXPnKC4uZtiwYdx22201rmvWrFloNBrKy8spKCggKqr6gTPRNKkKC9GeOIFpx0oMP76P+nQJbfL0qE7YoLTSA7XgitPiTI7GPiAJ+1WdsXW9DkuX/rgiYoKWXzRtDXqvV1JSQmFhodf7t2zZQnh4OOPGjaOkpIRJkyZVKfiHDx9m7NixtG/fviFPXy+2vn0py8ggfMkSSjIzfVLsjx07RmZmpuf2xcczcnNzWbRoEbGxsYwdO5YDBw6g1WoZOHAg/fv3Jz8/n8zMTK8FX6PRkJeXx8SJEwkNDSUxMbHRmUUA2e1oDx5E9/PPaH/5BU1uLpqjR9EePozmfLsNj3Ag3gbXAq2h7Lq7sXS4AWvKAFym1sFIL5qxeo/hu1wu8vPza+3xcv3119O7d2/P7YvHn48cOcK6desoKiqiR48e3HHHHQ3JXif6bdswrVhBSWYmphUrsPbp0+iiX9OQTnZ2tufriIgIYmPdB8ViY2Ox2WzExsayevVqtmzZgslkwm6vvWlUfHw8b7/9Nh9//DFLly7l+eefb1Rm4UMOB+q8PLRHj6I9fhzNsWNoDx1Ck5OD5sSJakXd0bo19g4dsAwejKNdO0Kty9DEFkA0YLrwOHtIW4qufymwP4toUeo8hl9QUMC5c+dISkoiNDQUtVrt9fEGgwFwXxxk0aJF3HfffVXu79OnD0OGDMFkMrFgwQK+/fZbevbsWeO6YmIa/vZV9eWXaMeMwf7uu4TcdBPOoUNp/eCD2FeuxHXTTXVej1ar9eSwWq3odLpquSIiIggJCSEmJga1Wu25PyQkhIiICD744AOuu+46hg8fzo4dO9i5cycxMTHodDoiIyOrrG/cuHFMnDiR5ORk4uLiMBqNNf4eKudSCiVmgnrmKi+HnBxUp06hOnAA1YEDkJ2N6tdfUf36KxQWonJeaB/gUqngqqtwJSXhuuYa7ImJ0L49rmuuwdWhA+j1qADd+Q/16nk1Pq3GmqOY350S/45KzATKzVWTOhX8nTt38umnn3ra9bpcLlQqFW+88YbX78nPz2fhwoUMHjyYfv36eZa7XC6GDRuGyeTetenRowdHjhzxWvAb05QobPNmbEuXYktJcc8pTklBv3Qp+s2bKa2hV703lZsjFRYWUl5eXi1XcXExVquV/Px8nE6n536r1UpxcTHdu3dn8eLFrF+/nlatWqFSqcjJyaG8vJyioqIq67vrrruYPHkyOp2OkJAQnnvuuRp/D0ps2qTETFApl8uFqrgY7dGjqIuK3EMtR46gPn0adXGxe0/95ElUlXoKurRanDExlLdrh+PGG3EkJuKIjcXRrh32pCQciYmg9fKvdPZstUWxIQlordnVljtCEhTzu1Pi31GJmUCZubwdE621W2aFp556ivnz5xMeHl6nJysqKmLWrFk8+uijpKamVrmvrKyMCRMmsHjxYkJCQli0aBEDBgygR48eNa5LumV6p8RcQc9ktaI9dAh1fj6aU6fQnDyJJicHY0EBjkOHUJ8+jeai408urRZHbCzO6GjsHTrgaNfOXdDj491F/Yor8OVFVI2/rSXiwCTUTrNnmVNtpLjTXzHH1TwRItCC/nesgRIzgTJzNahbZoU2bdoQGhpa5ydbt24dpaWlrFmzhjVr1gBw8803Y7VaSU9P5/7772fWrFlotVpSU1O9FnshauRyoTl5Et3eve7PP/+M5uRJ91776dOoLjo+4oyMRBUfjyM5Gdu11+JITsZ++eU44+LcRf3yy73voftBRVEPPzwfjTUHR0gCJe2zFFPsRfNVpz38b7/9lvXr19O1a1e0lf4x7r77br+GA9nDr40Sc/kzk/7rr4kcPx7tiRNVljtiY3G0bYv9iivcn9u3xxkX5y7ql12GKyJCkb8rUObfEJSZS4mZQJm5GrWHv379ekwmE2VlZT4NJURdhb30Eq3++lcA7AkJlD30ELa0NMo7dsQVLSccCVEXdSr4VquVOXOqN0wSIhCiMzIwfP45AKe+/BJ7pbYWQoi68z63spKEhASOHTvm7yxC1Kii2Of++KMUeyEaoU57+Pn5+WRlZREbG4uu0mwFpTdPE02fdt8+z9euyMggJhGi6avTHv4DDzzAtGnTePLJJ3n00Uc9H02Bce1aYtPSaHP55cSmpWFcu7ZR68vLy2NMLZ02CwoKWLx4MQB79+7l0KFDAEyfPr3e61y3bh2jRo1i9OjRbN++vVG5mypn60rtBZxyrVQhGqNBzdOaCuPatURMmoTa7J7vrM3OJmLSJADMXhq/NVZ0dDTPPvssABs2bGDgwIFceeWVzJ49u17rKS4u5v333+eNN97AZrPx8MMP07t37xZ3aUl1QUGVr51N5IxGIZSoWTfKDp8/31PsK6jNZsLnz/dJwc/MzKRDhw4cOXKEsrIyZsyYAcDs2bN55pln2LFjBwcPHiQ5OZnRo0ezdu1a9uzZw4oVKwCwWCw8//zzVYbJKkRERPDPf/7T00gtLCysxRV7AHvXrth69kT/7bdV9/aFEPVWpyGdpkrjZQ6/t+UN0blzZ1588UV69uzJpk2bPMs7depEWloaTz75ZJVrBxw9epQpU6awaNEi+vTpw1dffeU9v0bDunXreOqpp7jxxht9lrmp0Zw8CUArmSkmRKM064Lv8HLygbflDVFx5auKrpiXEhMTwyuvvML8+fPZs2fPJbtm3nHHHaxevZq9e/fy3Xff+SRzU3NmxQpcWq27qZkQosGadcEvycrCaTRWWeY0GinJyvLZc9Q2zKJSqXBedKBx4cKFTJ48maysLFq3bo23E52PHz/O9OnTcblcaLVa9Hp9ixzSAXDGxqKy23G0axfsKEI0ac16DL9inD58/nw0OTk4EhIoycry2wHbi3Xp0oXXX3+dNm3aeJYNGjSIMWPGEBYWRnR0NGcuviDGeUlJSVx55ZU89dRTqFQq0tLSuOaaawKSW2k05681YHrnHSwDB+K44orgBhKiiapTL51gkl463ikxl18yWa2EL1xI+NKlFP7tb5hvv10ZuXxActWdEjOBMnN566XTrId0RDMREoLlD38AwGkyXeLBQghvpOCLJkF1vnGf6/zV1IQQ9ScFXzQJqorzKaTgC9FgUvBFk2BPSsKl02F6551gRxGiyZKCL5oEx5VXYktLQ3P0aLCjCNFkScEXylbpxDRXSAgqOflKiAaTgi8US79lC206dCB8wQKw23HpdKis1mDHEqLJatYnXommLXTlSnC5CF+yhJAtW0ClkvYKQjSC7OELRVKVlmLYuJGyBx+kYOlStL/8gn7XLlR16FckhKiZ7OELRTJs2IDKYqHsjjsov/ZaTvfoQeQzz6ByOIIdTYgmSwq+UCTj+vXYExMp79ULAEdiImfWrgUp+EI0mF8Kvt1u5+9//zunT5+mvLycu+66i17n/3EBdu3axZo1a1Cr1QwYMID09HR/xBBNlPr0aUK2bKF0zBj3uH1lGk1wQgnRDPil4G/ZsoXw8HDGjRtHSUkJkyZN8hR8u93O8uXLmTdvHgaDgWnTptGrVy8i5QLV4jzjhx+icjgw33FHsKMI0az45aDt9ddfz/Dhwz23NZX2yrKzs4mPjycsLAytVkunTp3Yv3+/P2KIQHI6URUX+2RVxnXrKO/SBXunTj5ZnxDCzS8F32AwYDQaMZvNLFq0iPvuu89zn9lsxlSp46HRaKTsfGMs0XQZPvyQNldfjXbIEAwffdTg6ZOao0fR794dsGsWCNGS+O2gbX5+PgsXLmTw4MH069fPs9xoNGKxWDy3zWYzoaGhXtcTExPjr4h1ptVqFZHjYkrKpS4pAUB18CDRTz6Jq00bnI8+iuPRR+Hyy+u+ntdfB8Dw8MMYfPizKel3VZnkqjslZgLl5qqJXwp+UVERc+fO5dFHHyU1NbXKfW3btiU3N5fS0lIMBgP79+/n1ltv9bouJVxYQIkXOABl5Qo7e5ZWQPnevZR8/DGhy5cT8sILqOfPxzJoEGUjR2Lt1w/UtbypdLm47O23sfbuzRmTCXz4synpd1WZ5Ko7JWYCZebydgEUvxT8devWUVpaypo1a1izZg0AN998M1arlfT0dDIyMpg7dy5Op5MBAwYQHR3tjxgikCqmS+r1WNPTsaanozl+HNPbb2N6912M//kP9iuu4FxGBmX33osrKqraKrQ//oju118pevzxAIcXomXwS8F/5JFHeOSRR7ze36tXryrTNEUzUNHkrNIBekdSEiVTplAyYQLGTz7BtHw5EbNn0+qvf8X8hz9wbuRIyq+5xjP10rRuHS6tFvOwYcH4CYRo9qS1gvAJldPp/qKmefIhIZjvuIMz69dzauNGyu69F8Mnn3DZLbcQM3QopnfeQVVainH9eqwDBuCSd3xC+IUUfOEbdjsujab6iVIXP+zqqymeN4/fdu+m6IUXUJWXE/ncc8R164YmL48ymXsvhN9IawXhGw5Hvc6CdYWFUTZyJGUZGeh37sS0YgXaY8ewDh7sx5BCtGxS8IVPqBwO9x5+vb9RhS0tDVtamu9DCSGqkCEd4Rt2O2hl/0EIJZOCL3zD6ZTGZkIonBR84ROqioO2QgjFkoIvfMPhkCEdIRROCr7wDYej9rYJQoigk/9Q4RMqux2X7OELoWhS8IVvyEFbIRRPCr7wCTloK4TyScEXjWe3ozp3Tg7aCqFw8h8q6sdqRXfgALoffnB/7NuHbv9+VBYLtp49g51OCFELKfjCK1VZGdqffnIX9R9+QP/DD2gPHEB1vhWyMzyc8pQUzj30EOWpqVhvuAHpcymEcknBFwCoiovR/fjjhb32H35Ae+iQp+2xIzqa8tRULAMGUJ6SQnlqKo6kJJmKKUQTIgW/BVKfOeMp6hUFXnv0qOd+R3y8u7j/4Q/YUlMpT0nB2abNJVsfCyGUTQp+C6IqLOSyW29Fe/iwZ5k9OZnylBTK7ruP8ori3kQuyCyEqB8p+C2IuqAA7eHDlN11l7vAd+2KKyIi2LGEEAEiBb8FqbhwePk112Dr0yfIaYQQgSZH3FoQZ6tWgPsArRCi5ZGC35JotTjDwlAXFQU7iRAiCKTgtzDOyEgp+EK0UFLwWxhXRARqGdIRokXy60HbgwcPsnLlSmbOnFll+UcffcSmTZtodX5M+YknniAhIcGfUcR5zogIVLKHL0SL5LeC//7777N582YMBkO1+w4fPszYsWNp3769v55eeOGMjET766/BjiGECAK/DenExcUxceLEGu87cuQI69atY9q0aaxbt85fEUQNnJGRMqQjRAvlt4Lfu3dvNF76o/fp04c//vGPzJgxg59//plvv/3WXzHERVwREXLQVogWKuAnXrlcLoYNG4bJZAKgR48eHDlyhJ5eWuvGKOA0f61Wq4gcF2tILnXbtqisVmJCQ8FoVESmQJBc9aPEXErMBMrNVZOAF3yz2cyECRNYvHgxISEh7Nu3jwEDBnh9fH5+fgDT1SwmJkYROS7WkFwmrZZIoODXX90N0RSQKRAkV/0oMZcSM4Eyc3mbBBOwgr9161YsFgvp6encf//9zJo1C61WS2pqKj169AhUjBbPGRkJgLq42C8FXwihXH4t+LGxscydOxeAfv36eZb379+f/v37+/OphRfO883SZBxfiJZHTrxqYVyV9vCFEC2LFPwWpmJIR06+EqLlkYLfwsiQjhAtlxT8FsYVHo5LpZIhHSFaICn4LY1aLSdfCdFCyRWvWgqHA/233xKycSMqiwVVSUmwEwkhAkwKfjOmOnuWkC+/xPDZZ4Rs2oSmsBCXVoutd2/KHngg2PGEEAEmBb+Z0Rw9iuGzzzBs3Ij+v/9FZbfjjIzEMnAglkGDsN50E67zbamFEC2LFPymzm5Hv3s3IRs3Yti4Ed3BgwCUX3UVpU88gXXQIGw9e4KXRnZCiJZDCn4TVDFUo9myhfhPPkFdVOQZqikeMQJLejqOK64IdkwhhMJIwW8iNEePYji/F6//5htUdjuu1q0x33yze6jmxhtlqEYIUSsp+Eplt7tn1Zwfj/cM1XTsSOmTT2JNT6fVkCEUFRYGOagQoqmQgq8gnlk1Gzdi2LTJPVSj07mHah56yD1Uk5x84RtkXF4IUQ9S8IOspqEaR1QUlvR0LOnp7lk14eHBjimEaAak4Afa+aEaw8aNhGzciO78BcXLO3VyD9UMGoStRw/ZexdC+JwU/ABQnT1LyBdfuOfHXzxUk5FRfahGCCH8QAq+n2iOHLlwAtTFQzUVs2pkqEYIEUBS8H2ltqGaUaOwpqfLUI0QIqik4DeCt6Ea6/XXUyZDNUIIhZGCX0+aI0cuzKrZscM9VBMdLUM1QgjFk4J/KXY7qi1baLVmTY1DNZb0dMplqEYI0QRIwa+BqrjY01a4YqhGWzFUM3Kke6gmKSnYMYUQol6k4J/ndahm0CB0d91FfvfuuMLCgh1TCCEarOUWfLsd/a5dF2bVHDoEQHnnzu6hmkGDKO/eHTQaYmJicOXnBzmwEEI0jl8L/sGDB1m5ciUzZ86ssnzXrl2sWbMGtVrNgAEDSE9P99lzGteuJXz+fDQ5OTgSEijJysJ8551ApaGajRsxfPHFhVk1ffpQ9vDDWAYNwpGY6LMsQgihJH4r+O+//z6bN2/GYDBUWW6321m+fDnz5s3DYDAwbdo0evXqRWRkZKOf07h2LRGTJqE2mwHQZmcTMXGi+/J+v/1WbajGM6tGhmqEEC2A2l8rjouLY+LEidWWZ2dnEx8fT1hYGFqtlk6dOrF//36fPGf4/PmeYl9BbbViWrcOdUEBpaNGcfr99/ltzx6KlizBMmyYFHshRIvhtz383r17c+rUqWrLzWYzJpPJc9toNFJWVuZ1PTExMXV+Tk1OTo3LXSoVrr17CQFC6ry2C7Rabb1yBIoScykxE0iu+lJiLiVmAuXmqknAD9oajUYsFovnttlsJjQ01Ovj8+txsDQ2IQFtdna15Y6EhHqt52IxMTGN+n5/UWIuJWYCyVVfSsylxEygzFwJCQk1LvfbkI43bdu2JTc3l9LSUux2O/v376djx44+WXdJVhZOo7HKMqfRSElWlk/WL4QQTVnA9vC3bt2KxWIhPT2djIwM5s6di9PpZMCAAURHR/vkOSpm43ibpSOEEC2ZyuVyuYIdojY5XsblA0mJb9lAmbmUmAkkV30pMZcSM4EycylmSEcIIURwSMEXQogWQgq+EEK0EFLwhRCihZCCL4QQLYTiZ+kIIYTwDdnDF0KIFkIKvhBCtBBS8IUQooVouVe8uojT6eSNN97g2LFj6HQ6Ro0aRXx8vOf+rVu38sknn6BWq0lKSuLxxx9Hrfb/9vJSuSr84x//ICwsjAcffNDvmeqS69dff2XFihW4XC4iIyMZN24cer0+qJm2bNnCRx995LnwzuDBg/2a52LBuCBQY3IF6zVfW6YKgX69XypXMF7vDSF7+Oft3LmT8vJy5s6dywMPPMCKFSs899lsNlatWsWMGTP485//TFlZGbt37w56rgobN27k+PHjAclTl1wul4t//OMfjBkzhjlz5nDNNdcE5NTzS/2u3nrrLaZNm8acOXP48MMPKS0t9XumCu+//z7Lli2jvLy8yvKKCwJNnTqVWbNm8fnnn1NUVBT0XMF8zXvLVCEYr3fwnitYr/eGkIJ/3s8//8w111wDQMeOHTl0/hq34O53PWfOHEJC3N30nU4nOp0u6LkAfvnlFw4ePMigQYMCkqcuuXJzcwkPD+fjjz9mxowZlJaWeu3tEahMAMnJyZSVlWGz2QBQqVR+z1QhGBcEakyuYL7mvWWC4L3ewXuuYL3eG0IK/nkXX5hFrVbjcDg8X1dcgnHDhg1YLBZ+97vfBT1XYWEh7733Ho899lhAstQ119mzZzlw4ABDhgxh2rRp7Nu3jx9++CGomQASExOZPHkyEyZMoEePHrVeh8HXevfujUajqba8vhcEClSuYL7mvWUK5usdvOcK1uu9IWQM/zyj0Yi50uURXS5XlT+u0+nk7bffJjc3lwkTJgRs77C2XNu3b+fs2bPMmzePoqIirFYrbdu25aabbgpqrvDwcOLj47n88ssB6NatG4cPHyY1NTVomY4dO8bu3bv529/+hsFg4OWXX2b79u1cf/31fs10KfW9IFAgBes1700wX++1CdbrvSGk4J/XqVMnvv32W/r06cMvv/xCUlJSlftfe+01dDodzz33XMAOXF0q1+9//3t+//vfA/Dll1+SnZ0dsBd/bbni4uKwWCzk5eURHx/Pzz//zMCBA4OayWQyodfr0ev1qNVqIiIiOHfunN8zXUrlCwIZDAb279/PrbfeGuxYQPBe894E8/Vem2C93htCCv55aWlpfP/99/zpT3/C5XIxZswYz0Vb2rdvzxdffEHnzp2ZPXs24H7xpaWlBTVXoGdz1CfX6NGjeemllwD3eHqPHj2Cnik9PZ1p06ah1WqJi4sLarEIxAWBGpMrmK95b5mC+XqvSbBf7w0hrRWEEKKFCP77NCGEEAEhBV8IIVoIKfhCCNFCSMEXQogWQgq+EEK0EFLwhfCxc+fOMXHixGqtHYQINpmHL4QP7d69m+XLl3P69OlgRxGiGin4Qpz3448/8s477xATE0NOTg4Gg4HbbruNDRs2kJOTw3XXXce1117Lm2++yYsvvuj5nsq3N2zYwLhx41i8eHEwfxQhaiQFX4hKDh06xOOPP067du144YUXWL9+PTNnzqSsrIwnn3ySDh061Pr9U6dODVBSIepPCr4QlcTGxtKuXTvA3SPFZDKh1Wpp1aoVJpMpoD30hfA1OWgrRCUX93y/uB1uUlISlbuR2O32gOQSwhek4AtRD6GhoeTn51NcXIzL5WLbtm3BjiREncmQjhD1oFKpSE9PJysri6ioKHr06CHTL0WTId0yhRCihZAhHSGEaCGk4AshRAshBV8IIVoIKfhCCNFCSMEXQogWQgq+EEK0EFLwhRCihZCCL4QQLcT/BxWcsvIuG377AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "## sequence 1\n",
    "n1 = len(param1)-1\n",
    "plt.plot(param1[:,1],param1[:,2],color='orange',label=\"Path 1\")\n",
    "plt.plot(param1[n1,1],param1[n1,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",
    "n2 = len(param2)-1\n",
    "plt.plot(param2[:,1],param2[:,2],color='blue',label=\"Path 2\")\n",
    "plt.plot(param2[n2,1],param2[n2,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",
    "n3 = len(param3)-1\n",
    "plt.plot(param3[:,1],param3[:,2],color='red',label=\"Path 3\")\n",
    "plt.plot(param3[n3,1],param3[n3,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()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can see the starting points and their gradient helping to converge to one point"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Problem 3\n",
    "Confirm your answer to Question 2 using the built in Python or R optimizers. Also provide a confidence interval for each of the three parameters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "#First we redefine the log likelihood function that takes one parameter and the function calculates negative\n",
    "#log likelihood so that we minimize the function using optim\n",
    "\n",
    "def loglikehood_optim(param):\n",
    "    alpha = param[0]\n",
    "    mu1 = param[1]\n",
    "    mu2 = param[2]\n",
    "    z = 0\n",
    "    for i in range(10):\n",
    "        i_fact = scipy.special.factorial(i)\n",
    "        p = freq[i]*(np.log(alpha*np.exp(-mu1)*np.power(mu1,i)/i_fact + (1-alpha)*np.exp(-mu2)*np.power(mu2,i)/i_fact))\n",
    "        z=z+p\n",
    "    return -z"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "def gradient_ascent_fn(param):\n",
    "    lgrad_alpha= -lograd_alpha(param,freq)\n",
    "    lgrad_mu1 = -lograd_mu1(param,freq)\n",
    "    lgrad_mu2 = -lograd_mu2(param,freq)\n",
    "    return np.array((lgrad_alpha, lgrad_mu1, lgrad_mu2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Optimization terminated successfully.\n",
      "         Current function value: 1989.945860\n",
      "         Iterations: 67\n",
      "         Function evaluations: 121\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array([0.3598576 , 1.25603686, 2.66337702])"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#Now we input the initial value of parameters and the function that needs to be minimized in the 'optim' function\n",
    "param_initial = np.array((0.3,1.0,2.5))\n",
    "param = param_initial\n",
    "\n",
    "result = scipy.optimize.fmin(loglikehood_optim, param)\n",
    "result"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\jsruj\\Anaconda3\\envs\\Stat 654\\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": [
       "      fun: 1989.9458598829642\n",
       " hess_inv: array([[0.03778023, 0.06489559, 0.04602783],\n",
       "       [0.06489559, 0.12237734, 0.0755518 ],\n",
       "       [0.04602783, 0.0755518 , 0.06233602]])\n",
       "      jac: array([-9.78002144e-07,  4.75599606e-07,  1.69638014e-07])\n",
       "  message: 'Optimization terminated successfully.'\n",
       "     nfev: 17\n",
       "      nit: 13\n",
       "     njev: 17\n",
       "   status: 0\n",
       "  success: True\n",
       "        x: array([0.3598854 , 1.25609511, 2.66340436])"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#Using different in-built methods for optimization\n",
    "\n",
    "result = scipy.optimize.minimize(loglikehood_optim, param, method='BFGS',jac=gradient_ascent_fn)\n",
    "result\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* Now to get the confidence intervals, we need to calculate the standard error. \n",
    "* Standard error is obtained by taking the square root of the diagonal elements of inverse of hessian matrix\n",
    "* The standard error is then multiplied by 1.96 to obtain the 95% confidence interval"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Calculating confidence intervals"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>minimized</th>\n",
       "      <th>upper</th>\n",
       "      <th>lower</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>alpha</th>\n",
       "      <td>0.359885</td>\n",
       "      <td>0.740853</td>\n",
       "      <td>-0.021082</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>mu1</th>\n",
       "      <td>1.256095</td>\n",
       "      <td>1.941752</td>\n",
       "      <td>0.570439</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>mu2</th>\n",
       "      <td>2.663404</td>\n",
       "      <td>3.152761</td>\n",
       "      <td>2.174048</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "       minimized     upper     lower\n",
       "alpha   0.359885  0.740853 -0.021082\n",
       "mu1     1.256095  1.941752  0.570439\n",
       "mu2     2.663404  3.152761  2.174048"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sigma= np.sqrt(np.diagonal(result.hess_inv))\n",
    "upper = result.x+1.96*sigma\n",
    "lower = result.x-1.96*sigma\n",
    "import pandas as pd\n",
    "confidence_interval = pd.DataFrame({'minimized': result.x,\n",
    "                    'upper': upper,\n",
    "                    'lower': lower})\n",
    "confidence_interval.index=['alpha','mu1','mu2']\n",
    "confidence_interval"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>inbuilt_minimized</th>\n",
       "      <th>manual_minimized</th>\n",
       "      <th>upper</th>\n",
       "      <th>lower</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>$\u0007lpha$</th>\n",
       "      <td>0.359885</td>\n",
       "      <td>0.359885</td>\n",
       "      <td>0.740853</td>\n",
       "      <td>-0.021082</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>$\\mu_1$</th>\n",
       "      <td>1.256095</td>\n",
       "      <td>1.256094</td>\n",
       "      <td>1.941752</td>\n",
       "      <td>0.570439</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>$\\mu_2$</th>\n",
       "      <td>2.663404</td>\n",
       "      <td>2.663404</td>\n",
       "      <td>3.152761</td>\n",
       "      <td>2.174048</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "         inbuilt_minimized  manual_minimized     upper     lower\n",
       "$\u0007lpha$           0.359885          0.359885  0.740853 -0.021082\n",
       "$\\mu_1$           1.256095          1.256094  1.941752  0.570439\n",
       "$\\mu_2$           2.663404          2.663404  3.152761  2.174048"
      ]
     },
     "execution_count": 55,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "confidence_interval = pd.DataFrame({'inbuilt_minimized': result.x, 'manual_minimized' : param1[1000,:],\n",
    "                    'upper': upper,\n",
    "                    'lower': lower})\n",
    "confidence_interval.index=['$\\alpha$', '$\\mu_1$' , '$\\mu_2$']\n",
    "confidence_interval"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Problem 4\n",
    "Implement the EM algorithm approach detailed in Question 11 to fit the parameter  𝛼 ,  𝜇1  and  𝜇2 .\n",
    "\n",
    "Compare the number of steps and the computation time of the EM algorithm to the gradient ascent method you used in Problem 2.\n",
    "State which is faster and why you think this is the case.\n",
    "Lastly, name one advantage of each method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Values given in the textbook\n",
    "alpha = 0.3\n",
    "mu1 = 1\n",
    "mu2 = 2.5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [],
   "source": [
    "def Ez(alpha,mu1,mu2,k):\n",
    "    term1 = alpha*np.exp(-mu1)*np.power(mu1,k)\n",
    "    term2 = (1-alpha)*np.exp(-mu2)*np.power(mu2,k)\n",
    "    return term1/(term1+term2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([6.57619125e-01, 4.34481758e-01, 2.35073895e-01, 1.09469629e-01,\n",
       "       4.68661004e-02, 1.92888357e-02, 7.80587401e-03, 3.13704202e-03,\n",
       "       1.25718311e-03, 5.03252852e-04])"
      ]
     },
     "execution_count": 36,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Check function\n",
    "Ez_vec = np.array(list(map(lambda y: Ez(alpha,mu1,mu2,y), np.arange(10))))\n",
    "Ez_vec"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Run EM algorithm\n",
    "def em_poisson_alg(alpha,mu1,mu2,X,para,liketol):\n",
    "    old_likelihood = -np.inf\n",
    "    likelihood_diff = np.inf\n",
    "    iter_count = 0\n",
    "    para_diff = np.inf\n",
    "\n",
    "    while para_diff > para or likelihood_diff > liketol:\n",
    "        mu11 = mu1\n",
    "        mu22 = mu2\n",
    "        alpha_a = alpha\n",
    "        Ez_vec = np.array(list(map(lambda y: Ez(alpha,mu1,mu2,y), np.arange(10))))\n",
    "        alpha = sum(X*Ez_vec)/sum(X)\n",
    "        mu1 = sum(X*np.arange(10)*Ez_vec)/sum(X*Ez_vec)\n",
    "        mu2 = sum(X*np.arange(10)*(1-Ez_vec))/sum(X*(1-Ez_vec))\n",
    "        # Check convergence\n",
    "        para_diff = sum(np.absolute([alpha-alpha_a,mu1-mu11,mu2-mu22]))\n",
    "        new_likelihood = log_likelihood([alpha,mu1,mu2],X)\n",
    "        likelihood_diff = np.absolute(new_likelihood-old_likelihood)\n",
    "        old_likelihood = new_likelihood\n",
    "        iter_count = iter_count + 1\n",
    "    return [alpha,mu1,mu2,iter_count,para_diff,likelihood_diff]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [],
   "source": [
    "em_out = em_poisson_alg(alpha,mu1,mu2,freq,10e-6,10e-6)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[0.3593075138402641, 1.255089824982245, 2.662698781023993]"
      ]
     },
     "execution_count": 41,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "em_out[0:3]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The $\\alpha$, $\\mu_1$ and $\\mu_2$ values are 0.3593, 1.2550 and 2.6626 respectively which are very close to the actual value"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[0.3593075138402641,\n",
       " 1.255089824982245,\n",
       " 2.662698781023993,\n",
       " 1116,\n",
       " 9.968900464174535e-06,\n",
       " 3.8574398786295205e-08]"
      ]
     },
     "execution_count": 42,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "em_out"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [],
   "source": [
    "##Calculating time"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.02 s ± 110 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n"
     ]
    }
   ],
   "source": [
    "%timeit em_out = em_poisson_alg(alpha,mu1,mu2,freq,10e-6,np.inf) # No likelihood tolerance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\jsruj\\Anaconda3\\envs\\Stat 654\\lib\\site-packages\\ipykernel_launcher.py:8: RuntimeWarning: invalid value encountered in log\n",
      "  \n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "26.3 s ± 3.41 s per loop (mean ± std. dev. of 7 runs, 1 loop each)\n"
     ]
    }
   ],
   "source": [
    "param_initial_1 = np.array((0.3,1.0,2.5))\n",
    "%timeit param1 = gradient_ascent(param_initial_1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {},
   "outputs": [],
   "source": [
    "em_steps = np.zeros(100)\n",
    "for ii in np.arange(100):\n",
    "    alpha = np.random.uniform(0,1,1)[0]\n",
    "    mu1 = np.random.uniform(0,5,1)[0]\n",
    "    mu2 = np.random.uniform(0,5,1)[0]\n",
    "    em_steps[ii] = em_poisson_alg(alpha,mu1,mu2,freq,10e-6,np.inf)[3]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The EM algorithm seems to be a bit faster here. However this could be due to the implementation because it takes nearly 3 times as many iterations to converge (though we should not forget that the gradient ascent algorithm may do one or more halvings within each step)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Checking if the mixture model is a good fit:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.9997\n"
     ]
    }
   ],
   "source": [
    "deaths = np.arange(10)\n",
    "freq = np.array([162, 267, 271, 185, 111, 61, 27, 8, 3, 1])\n",
    "i, ni = deaths, freq\n",
    "i_fac = scipy.special.factorial(i)\n",
    "alpha0, mu10, mu20 = [0.35987258442400516, 1.2560722507101618, 2.6633886707631791]\n",
    "pois_mle = (alpha0*np.exp(-mu10)*mu10**i+(1-alpha0)*np.exp(-mu20)*mu20**i)/i_fac\n",
    "print(f'{np.sum(pois_mle):.4}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAECCAYAAADzStBRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3deVxU9f7H8dds7JuIgKDoRQQXFCuvuVWu5TWzGy1mWf70lhcxMzfU1HIjyqVcEtcWyTK9Zt5c8uZNrVzKJQU3xBVMUEFURIaBWX5/cB0lEFEZDjCf5+Ph4+GcM+fMmy/DfOZ7lu9XZbFYLAghhLBbaqUDCCGEUJYUAiGEsHNSCIQQws5JIRBCCDsnhUAIIeycFAIhhLBzUgiEEMLOaZUOcC/S09PveVsfHx+ysrIqME31JW1RnLTHTdIWxdWE9ggICLjtOukRCCGEnZNCIIQQdk4KgRBC2DmbnSMwm80sXbqU1NRUdDodUVFR+Pv7W9dv376djRs3olarCQoK4rXXXkOtVhMTE4OLiwsAvr6+REdH2yqiEEIIbFgI9uzZQ2FhIbGxsaSkpJCQkEBMTAwABQUFrFy5kpkzZ+Lo6Mjs2bP5/fffadmyJQCTJk2yVSwhhBB/YrNDQ8nJybRq1QqA0NBQTp48aV2n1WqZOnUqjo6OQFHvQafTkZqaisFgYNq0aUyePJmUlBRbxat08Ynx7EjfUWzZjvQdxCfGV3oW5wtr8N3VBt1qJ3x3tcH5wppKzyCEqDps1iPQ6/XWQzwAarUak8mERqNBrVbj5eUFwPfff09+fj4tW7bk7NmzPPXUU3Tt2pWMjAzi4uKYPXs2Go2m2L59fHzuOZdWq72v7e/VoyGP8vK3L/PlM1/SqUEntqVuI3pLNF8+82Wl5lGnrUCTMgaVKQ8AreEcXiljcHd3xxzUt9JyVEVKvTeqImmL4mp6e9isEDg7O6PX662PLRZLsQ90s9nM8uXLycjIYOTIkahUKurWrYu/vz8qlYqAgADc3Ny4fPlyiV/A/VzPq9T1wOGu4cR3iafvmr485PsQW85uoVO9Tvye+jumPBNNvJvgqHG0eQ7fxPHWInCDypQHiePJculu89evymrCteIVIT4xnkdDHiXcNdy6bEf6DhIzE4mOuLdzdgcOHGD48OFMnDiRLl26WJf/4x//oHHjxowdO5Z33nmHKVOmlLr9qVOnuHbtGhEREff0+uW1du1a1q9fz0svvWTNeeDAAaZOnUr9+vVRqVQYDAa6detGZGRkqfvIzs5m2bJlDB8+3KZZ75Yi9xGEhYWxf/9+AFJSUggKCiq2fvHixRQWFjJ69GjrIaKtW7eSkJAAFDWmXq+nVq1atopY6bppLnA8qJBNDps50xAa5u5k3I5x9FzbkyafNyG3IBeAg1kHScxMxGAyVHgGjaH0m/Fut1zYn4g6Ebz87cvWQ5k70ncQ9WMUEXXu70M4KCiILVu2WB+fOnWK/Px86+PbFQGAn3/+mdTU1Pt6/fL45ZdfePvtt4sVK4A2bdowe/ZsPvroI2bPns2qVavIzc0tdR/e3t5Vrgjcic16BG3atCEpKYkJEyZgsViIjo5m+/bt5OfnExwczNatW2nSpIn1l9+zZ0+6dOnC/PnzmThxIiqVisGDB5c4LFRdOV9Yg1vySHSWAlBBPa2JhXVMjHloElstdTmTcwY3BzcAPvz9Q35I/QGdWkdYrTAi6kTQ2q81L4S+cN85TI4BaA3nSl0u7Mdz658rsaxXcC/+r9n/8aDvgwS4B/DS9y/h5+LHhbwLhHqF8kfuHwBk52cz6L+Dim27utfqO75mo0aN+OOPP8jNzcXNzY3NmzfTtWtXLl68CEBkZCT/+te/GDZsGK+++iohISGMGDGC999/n02bNqHT6WjcuDGTJ08mISEBBwcHFi9eTFBQEP7+/ixatAidTkevXr3w9fXlk08+Qa1WExAQwMiRI9Fqb37cnT9/nunTp2MymQAYOnQoR48e5dixY8yYMYN33nmHunXrlvpz6PV6NBoNGo2G48ePM3fuXNRqNQ4ODowaNQqLxcKUKVOIj49n6dKl7N+/H7PZTNeuXXnuuedYu3Yt//nPf1Cr1bRo0YKoqKhS84SEhNCvXz/Cw8M5e/YstWrVYvLkyaSnp/PBBx9YM4wbN446deqU47d+ezYrBGq1mkGDir9ZAgMDrf9fuXJlqdsNGzbMVpEU5XB8clERuIXWUkDAuY/p1TGx2PKp7aYSGRLJwayDJGUlseH0BpKzk62FYPhPw9GpdbT0aUlLn5aEeYeV+7DSteCxeB6LQW2+edjOrHbmWvDY+/wJRU3i5eSFn4sf53LPEegWiKejZ4Xs95FHHuGXX36hR48eJCcn07dv32K9BI1Gw4QJExg3bhy1a9dm8ODB+Pv706NHD7y9vWnatOlt911QUMCCBQuwWCy8+uqrzJ07l1q1avHpp5+yadMmevXqZX3uggULiIyMpGPHjpw4cYIZM2awaNEifvzxR0aMGFGiCOzevZu33noLtVqNRqNh6NChODs7M3PmTEaPHk1ISAjbt28nPj6ewYMHW7f74YcfmDNnDrVr12bTpk0AbNq0iTfffJNmzZrx73//G5PJdNs8GRkZfPjhh/j6+vLGG29w7Ngxjh07RmhoKNHR0SQlJZGbm1t1C4EozsVY+rFnF2MWV/+0rJ57Peq51+Op4KeAovMrVwuuWv+fnZ/NnvN7+DL5SwB0ah0Dmw/knbbvAHDk0hFCvEJw0DiUeD29X9FxTfdT76MxpGNyDOBa8FjrcmEfyvoG76x1ZnzH8fRd05e3HniLhKMJDH9wOB0COgDg7eRdrh5Aabp27cpHH31E3bp1adGiRanP8ff3p0WLFhw+fJg2bdqUub9bp1yvX78+AFeuXOHSpUtMnjwZAIPBwF//+tdi26WmplrPN4SEhJCZmVnm67Rp04YxY8aUWH7p0iVCQkIAiIiIYMmSJcXWT5w4kSVLlpCdnc3DDz8MwJgxY1i5ciWLFi2iefPmWCyW2+bx9PTE19cXKLqvqqCggCeffJIVK1YQExODm5sbr732WpnZy0MKQSXJUXvhab5SYrnJMbCUZxenUqnwcvSy/n/ZE8uK3jzXUknKTOJg1kGa124OFHXbu6/pjoPagSbeTWjh04KIOhE8Gvgo9d3rE58YT0SdCDq02209ObojfQeJifH3fCJQ1Cw70ncQvSWahV0X0iGgA+0D2hP1Y5T18f0ICAggPz+fNWvW8Prrr5ORkVHiOUeOHOH06dNERESwatUq+vTpg0qlwmw2A+Dg4MClS5fw9/fnxIkTNGjQACg6CgFFH5516tRh2rRpuLm5sWPHDpydnYu9RoMGDUhKSqJDhw6cOHHins9F1q5dm5MnT9KoUSMSExOpV6+edV1BQQE//fQTEydOxGKxMGDAADp37sz69esZMWIEDg4OjB49mkOHDt1Vnu3bt9OiRQv69+/Pjz/+yIoVK0otUndDCkElsFgsTL3izFT3qzirbn6DuZ9DMiqVioYeDWno0ZDejXpblztpnFjYdSFJmUkkZSWx7tQ6vkz+kpmPzKRvk774ufjxyqZX6N+sP8M7DOdg+kHrH7kQAImZiXz5zJfWq4Y6BHRgYdeFJGYm3nchAOjcuTObN2+mfv36JQpBbm4u06dPZ+rUqfj5+TF48GBatWpFaGgoixYtokGDBrz44ouMHTsWf39/3N3dS+xfrVbzxhtvMG7cOCwWCy4uLowbN67YcwYPHszMmTNZtWoVRqPRerPr3Ro1ahRz5861XhU5evRo6zoHBwfc3d157bXXcHd3p3Xr1vj5+REcHExUVBReXl74+PjQrFkz/P39y50nLCyM9957j88//xyVSsWQIUPuKfutVJZb+1bVRHUbhjo1J5Un1jzBioineCL/p0o9JHOj5+Dp4Ektp1psO7uNQf8dxHXjdbQqLa46V5Z0X1Ihf+DVnVw+epO0RXE1oT3KunxUegSVoIFHA/a9vA+NSsNFrVOlvvaNnsMNnep34tj/HWPs9rEsT16Om86NdnXbVWomIUTVIqOP2pjBZMBiseCqc8WpkovA7ezM2MnGMxt5MuRJzl0/x6KDi5SOJIRQkBQCG5t3YB7d13Qn35h/5ydXghs3By3supBvnvuGae2nlToOkhDCfkghsKECUwHLjy4nwDWgyvQGEjMTrVd/qFQqBjQfwMKuC9lydgvV8HSREKICyDkCG9p4eiOZ+kwGNB+gdBSr0i4RddI4sfTgUprUasLzoc8rkEoIoSTpEdjQZ0c+o6FHQx6r95jSUcr0gO8DPOD7AJN+ncTFvItKxxFCVDIpBDZy6NIh9l7YS/9m/VGrqnYzq1VqZj46E71Rz/id45WOI6oAddoKfHe1oe62ehU2Z0VGRgbvvPMOb731Fm+88QYfffQReXl5JZ534sQJli1bdtv97N69m3Xr1t3167/44osUFBTcdv3ixYutw0CU5sKFC+zcufOuX7c6kENDNhLqFcrHnT+mc/3OSkcplxCvEEY8OIK4PXFsOL2BJ//ypNKRhEKcL6wpMWeF57GiG5zu9b4Xg8HA+PHjGTVqFM2aNQOKxtyZOnUqcXFxxZ4bEhJiHbahNHcadsJW9u/fT1paGu3bt1fk9W1JCoGNOGgceCbkGaVj3JWollH8J/U/nMstOTqpsB/up94vMWeF2qzH/dT791wIfv31VyIiIqxFAKBHjx589913pKenk5CQQE5ODjk5OfTp04etW7fyzjvvsGHDBtauXYu7uzs6nY7OnYu+WKWlpdG7d2+mTp1KnTp1SE9Pp2nTpgwfPpzMzEw++ugjCgoKyMnJ4dVXX6Vjx46l5vrpp59Yvnw5Xl5eFBYWEhQUhMlk4sMPP+TixYvk5OTw8MMPM3LkSL766isMBgPNmzfH1dXVOmR+fn4+48aNs451VB1JIbCBlSkryczLJDoiusofFrqVVq3l26e+RauWt4U9s8WcFenp6aXe2Vq3bl3rMNQPPPAAzz//PAcOHADg6tWrfP311yxZsgSdTseIESNKbP/HH38wY8YMHB0defnll8nOziYtLY0XXniBVq1acejQIT7//PPbFoLFixcTHx+Ph4eHdRiKzMxMmjVrxujRoykoKOD5558nJiaGl156ibS0NDp06MDatWt5++238fHxYfny5fz000/069fvnttHafIXX8HMFjNzfp+Dn4sfb7R6Q+k4d+1GEdhydgs6tY5HAh9ROJGobLaYs8LHx4fk5OQSy//44w/r6Jp//kZ97tw5GjRogJNT0aXXzZs3L7F9QECAdUpcb29vCgoK8Pb2Zvny5WzcuBEAo9FYaqbs7GxcXFzw9PQstn93d3eSk5PZv38/rq6uFBYWlvrzzJs3D2dnZ7KysggPDy/xnOqk+nxdrSa2nt1K6rVU/q/5/ykd5Z4ZzUam/DqFET+N4FrBNaXjiEp2LXgsFo1LsWX3O2dFx44d2bdvH0ePHrUu27BhA15eXtaewo3RQ28IDAwkLS0Ng8GA2WwutZCoVKoSyz777DMef/xx3n77bR544IHb3h/j6enJ9evXuXKlaFTgG/vftGkTbm5uTJgwgRdeeAGDoWh0gFtHQJ05cyZjxoxh7Nix1K5du9rfgyM9ggr2+ZHP8XPx428N/6Z0lHumVWuZ+ehM/v7d33lv93vEdYy780aixtD7RRaN6pk4vsIGSHR2diY2Npb4+HiuXr2KyWSiUaNGTJw48bbbeHp60rdvX9588008PDwwGAxoNBrrLF6389hjjzFv3jy+/PJL6tSpw9Wrf57xo4hGoyEmJoaYmBjc3d2tM5g9+OCDTJ06lYMHD+Lk5ERgYCAXL14kODiYL7/8ktDQULp37050dDRubm54e3tz6dKle26bqkBGH61AZ3LO0HFlR4Y/OJyRD420yWtUpDu1xbu73mXpoaWs7rXaLgamqwkjTFaUqtAWJpOJFStWWI+9Dxs2jIEDB9p8AvvSVIX2uF+KTF5vj/KN+XSp34WXm7ysdJQKMab1GBq4N2DUz6PQG/V33kCICqTRaNDr9QwaNIjo6GhCQkJo2bKl0rFqJOkR2LHytMXO9J2k5qTyYtiLpR6PrUnkvXGTtEVxNaE9ZD6CSnAg8wA+Tj7Uc6935ydXI+0D2tM+oOgGmhsnzIQQNYscGqoAFouFmF9iGLh5oNJRbGbtibU8s+4ZCky3v0VfCFE9SSGoAHsv7OXwpcO82vRVpaPYjKvOlT0X9vDxgY+VjiKEqGBSCCrAZ0c+w8PBg8gQ284/rKTuDbrzTKNnmHtgLsnZJa/nFkJUX1II7tOFvAtsOLWBF0JfwEXncucNqrEp7afg7uDOqJ9HYTKXfS23EKL6kEJwn/ac34NapaZ/s/5KR7E5bydvprWfxv7M/fx87mel4wghKohcNXSfegX3okNAB2o51VI6SqXoHdybv3j8hZZ15HpuIWoK6RHchxs3WdlLEYCisV1uFIETV05gtpgVTiSEuF9SCO7DK5te4a1tbykdQxFJmUl0Wd2F5UeXKx1FCHGfpBDco6PZR9mVsYuwWmFKR1FEC58WtKvbjtjdsTKRjRDVnBSCe/T54c9x0jjRJ6yP0lEUoVKpmPHIDEwWE2O3j632w/AKYc+kENyDq4arfHPiG55u9DTeTt5Kx1FMkEcQY/86li1nt/DNiW+UjiOEuEc2u2rIbDazdOlSUlNT0el0REVF4e/vb12/fft2Nm7ciFqtJigoiNdeew2gzG2qitXHV6M36hnQfIDSURQ3oNkANp3ZRJa+eg/IJYQ9s1kh2LNnD4WFhcTGxpKSkkJCQgIxMTEAFBQUsHLlSmbOnImjoyOzZ8/m999/x2Qy3XabquS5xs9Ry6kWLXxaKB1FcRq1hpU9V6JRa5SOIoS4RzY7NJScnEyrVq0ACA0N5eTJk9Z1Wq2WqVOn4ujoCBT1HnQ6XZnbVCWejp41ejiJu3WjCPyY9iNbzm5ROI0Q4m7ZrEeg1+utk0pD0XykJpMJjUaDWq3Gy8sLgO+//578/HxatmzJrl27brvNrXx8fO45l1arva/t3976Ng/5P8SzTZ+9531UFffbFrcymU3M+m4WF3IvcOD1A9Ryrn73VlRke1R30hbF1fT2sFkhcHZ2Rq+/OauVxWIp9oFuNptZvnw5GRkZjBw5EpVKdcdtbrifCSLuZ4KJs9fO8uGvHzK01VAeq/PYPWeoKip6so3p7afTc21P3tz4Jh899lGF7bey1ITJRyqKtEVxNaE9FJmqMiwsjP379wOQkpJCUFBQsfWLFy+msLCQ0aNHWw8R3WkbpX1x9AtUKhX9mvZTOkqVFO4TzuCIwaxKWcW2s9uUjiOEKCebTVV546qhtLQ0LBYL0dHRnD59mvz8fIKDgxk3bhxNmjSxznjVs2dPWrduXWKbwMDAEvtWYqrKfGM+rb9qTbu67VjSfck9v35VYotvOfnGfJ749gn0Rj3bnttWrUZkrQnf+iqKtEVxNaE9FJmqUq1WM2jQoGLLbv1QX7lyZanb/XmbquK7U99x2XDZLkYZvR9OWidmPTqLE1dO4Kx1VjqOEKIcZPTRcqrlWIvewb3pENBB6ShVXmu/1rT2aw2A2WJGrZL7FoWoyuQvtJy6N+jOgq4LZPL2u/Dvk/+m19pe1lFahRBVkxSCcth6ditXDVeVjlHteDt5k5iVyIf7PlQ6ihCiDFII7uCS/hIDfxjIrH2zlI5S7TwS+Ah9w/qy8OBCEjMTlY4jhLgNKQR38NWxrygwF/BK01eUjlItTXx4Ir7Ovoz8eSQFpgKl4wghSiGFoAxGs5GEIwl0DOhI41qNlY5TLXk6ehLXMY6j2UfZenar0nGEEKWQQlCGzambSb+eLqOM3qfHGzzOD5E/8ETDJ5SOIoQohRSCMuy7uI9At0C6BXVTOkq117x2cwCOZR/DZDYpnEYIcSspBGWY8PAENkduRquW2y0qwpFLR+i+pjtLDy1VOooQ4hZSCG7jxrXvno6eCiepOZp6N6Vz/c5M3zud01dPKx1HCPE/UghKca3gGq2/as3nRz5XOkqNolKpeL/j++jUOkb/MlrmORaiipBCUIpvjn/DFcMVInwilI5S49R1rcvEthPZlbGLL5O/VDqOEAIpBCVYLBY+P/I5ET4RPOD7gNJxaqSXwl4i2DOYg1kHiy3fkb6D+MR4hVIJYb+kEPzJjvQdHL9ynP7NZZRRW1GpVMR1iGPjmY3sSN8BFLV71I9RRNSRXpgQlU0uh/mTZUeXWUcaFbbTMbAjC7su5B8//IPH6j3GzoydLOy6UEZ3FUIBUgj+5N2H3+X4leMyln4laF+3PW46N9afXk9UyygpAkIoRA4N/Uk993p0rt9Z6Rh2YWfGTvKMeQB8dvgz62EiIUTlkkLwPwaTgegt0ey/uF/pKHZhR/oONu8awLlGOsyNIbm+gU07XpViIIQCpBD8z4bTG/j3yX+TU5CjdBS7YPnjS+J9CnE1ZqECGupgno8By9nlSkcTwu5IIfifzw9/TrBnMI8EPqJ0FLsQadqL1lJ8WGoXlYVI8z6FEglhv6QQAAezDrLv4j76N+sv8+tWEo0h/a6WCyFsRz71KOoNuGhdeL7x80pHsRsmx4BSl+eoPYn9LbaS0whh36QQAE28mzC45WAZYK4SXQsei1ld/BJds9qZVapWLEhaIFNbClGJ5D4C4PUWrysdwe7o/SIBcD/1PhpDOibHAK4Fj+WRWt2ok/oYb+94m3VPr5NDdUJUArv+KzOZTaw/tV7m0lWI3i+Si+12k9HpDy62243eLxIPBw8mPDyBA5kHWHFshdIRhbALdl0Itpzdwj9//Cc/pP6gdBRxi8iQSB72f5i43XFcL7yudBwhajy7PjS07Mgy/F38ZS7dKubGvAWZ+kxcda5KxxGixrPbHsGpq6fY+sdW+jXth06tUzqO+JPQWqHWsYeMZqPCaYSo2ey2ECw7sgydWsfLTV5WOooow9z9c3lm3TMy4b0QNmSXhcBisZCYmUjPv/TE18VX6TiiDEHuQfx+8Xe+OvaV0lGEqLHKfY4gMzOT3NzcYvPMBgcH2ySUralUKr596ltyC3OVjiLu4OlGT7M8eTnv73mfJ//yJN5O3kpHEqLGKVchWLlyJevWrcPT8+YNVyqVio8//thmwWzFYrGQV5iHi84Fdwd3peOIO1CpVMS2j+XxNY8TtzuOGY/OUDqSEDVOuQrBzz//zNy5c/H2Lv+3MbPZzNKlS0lNTUWn0xEVFYW/v3+x5xgMBqZNm0ZUVBSBgYEAxMTE4OLiAoCvry/R0dHlfs3y2PHHDv6+8u8s/9tyWvu1rtB9C9sI8w7jtRav8dnhzxj50Ej8Xf3vvJEQotzKVQh8fHzuqggA7Nmzh8LCQmJjY0lJSSEhIYGYmBjr+pMnT7JkyRIuXbpkXVZQUHRj16RJk+7qte7Ggn0LUKvUNK/d3GavISre8AeG0zesrxQBIWygXCeLw8PDWb58OcnJyZw6dcr6ryzJycm0atUKgNDQUE6ePFlsfWFhIaNGjbL2BABSU1OtvYTJkyeTkpJytz9PqeIT40lLfp/aOx9iZeFqTjYwkXliDvGJ8RWyf2F7bg5uhHiFAJBxPUPhNELULOXqEWzbtg2AXbt2WZfd6RyBXq+3HuIBUKvVmEwmNBoNAE2aNCmxjaOjI0899RRdu3YlIyODuLg4Zs+ebd3mBh8fn/LEtnrBy0jIyXk4qgEV1CYX5/R5uDQac9f7qkm0Wm21+/nf2/4ec/bM4dA/D1HHpU6F7rs6toetSFsUV9Pbo1yFYP78+Xe9Y2dnZ/R6vfWxxWIp8YH+Z3Xr1sXf3x+VSkVAQABubm5cvny5xC8gKyvrrrI0S1+O9k99Hxd10fKLDd68q33VJD4+Pnfdlkrr7NeZWEMso74fxazHZlXovqtje9iKtEVxNaE9AgJKH/odynloKD8/n6VLlzJkyBCioqKIj48nLy+vzG3CwsLYv79o/t+UlBSCgoLu+Dpbt24lISEBgOzsbPR6PbVq1SpPxDLJJCg1R+NajRnUYhBfp3zN3gt7lY4jRI1QrkKwbNkyCgsLGT16NDExMahUKj799NMyt2nTpg06nY4JEyawbNky+vfvz/bt2/nvf/972226dOnC9evXmThxIrNnz2bw4MF37EWUx+0mQbndclG1vfXgW/i7+vP2jrfljmMhKkC5Dg2dOHGCGTNuXr/9z3/+k5EjR5a5jVqtZtCgQcWW3Xpi+IZbrxDSarUMGzasPJHuyt5akbRMn4fLLWUvzwxJtSK5cz9FVDWuOlcmtZ3EuO3jOHn1JKG1QpWOJES1Vq4egclkwmw2Wx9bLBbU6uozOsV6gwdJAUMxOgZiQYXRMZCkgKGsN3goHU3co15/6cWOPjukCAhRAcrVIwgPD2f27Nl0794dlUrFDz/8QPPm1ec6/OiIopvSLjLWetInCKjYW9VEZVKpVHg6emIym9iZsZNHAh9ROpIQ1Va5vtb379+fevXqsWLFCpYvX05AQAD9+vWzdTYh7ujTw5/y4sYX2XN+j9JRhKi2VJZbR5GrJtLT7/1qn5pwGVhFqQltkVeYx2P/egwvRy++f+Z7tOp7n2upJrRHRZG2KK4mtEdZl4+W+VczceJEpk6dyquvvopKpSqxftmyZfefToj74KJzYVK7SQz67yASjiQwMHyg0pGEqHbKLAQjRowAYNaskjfuVMOOhKihejbsyWOBjzF973SeCn6qwu84FqKmK/McwY2buZYsWUKdOnWK/ZszZ06lBBTiTlQqFVPbT6Weez0u6i8qHUeIaqfMHsGsWbPIyMjgwoULjBo1yrrcZDKh1dr1vPeiimnk1YjNkZtLPYQphChbmZ/mr7zyCpmZmSxatIiBA28ee1Wr1dSrV8/m4YS4GyqVityCXL44+gWvt3j9vk4cC2FPyvxL8fX1xdfXl5CQEJo1a1ZZmYS4ZzvSdzBt9zR0Gh2vhb+mdBwhqoVy3Udw9uxZOTksqoXHGzxO53qdmbl3JhfyLigdR4hqoVx9Zy8vL0aMGEHjxo1xcnKyLr/1cJEQVYFKpWJK+yl0Xd2Vab9NY17neUpHEqLKKwOjlgUAAB4sSURBVFePIDQ0lPbt21OnTh3c3d2t/4SoioI9gxkcMZg1J9bwa8avSscRosorV4/g+eefJz8/n1OnTmE0GmncuDHOzs62zibEPRvaaih/XPsDH+eaO6uUEBXlroah9vT0xGw2c+nSJcaOHUtYWJit8wlxT5y1zsztPFfpGEJUC+UqBF988QVDhw4lPDwcgEOHDpGQkEBsbKxNwwlxvzKuZ/De7vcY32Y8/q7+SscRokoqVyHQ6/XWIgBFw1IbDAabhRKiohhMBjac3oDFYuHjLh8rHUeIKqlcJ4tVKhWZmZnWxxcvXqxWE9MI+9XQoyHREdF8e/JbdqbvVDqOEFVSuXoEzz77LOPHj6dFixYAJCUl8Y9//MOmwYSoKEMihrA6ZTXjd4znh2d/QKfWKR1JiCqlXIWgTZs21KtXj0OHDmE2m3nmmWdkiAlRbThrnZnSfgoDfhjA54c/5/UWrysdSYgqpdzHd86fP8+5c+e4cOECV69etWUmISrc4w0eZ2q7qTzb+FmlowhR5ZSrEKxatYovvvgCFxcXHB0dWbx4MRs3brR1NiEq1MDwgXg7ectwKUL8SbkKwS+//EJsbCx9+vThxRdfJDY2ls2bN9s6mxAV7lzuOXr9uxfbz21XOooQVUa5CoG7u3uxMYZcXV2LPRaiuvB28uZy/mUm7JxAgalA6ThCVAnlKgRNmjRh+vTp7Nu3jwMHDjBv3jxq167Nb7/9xm+//WbrjEJUmBsnjo9fOc4nhz5ROo4QVUK5rho6ffo0AOvXry+2fNOmTQA8/PDDFRxLCNvpFtSNxxs8zoe/f8jTjZ4mwC1A6UhCKKpcheDdd98FiqaotFgsMk2lqPYmt51M59WdWXhwIVPaTVE6jhCKKtcn+tWrV5k/fz6HDh3CZDLRrFkzhg4dire3t63zCWETQR5BrOi5gog6EUpHEUJx5TpH8Mknn9C4cWOWLFnC0qVLadq0KUuXLrV1NiFsqo1/Gxw1juQV5smJY2HXylUIMjIyeP7553F1dcXd3Z0XXniB8+fP2zqbEDY3Y+8M2n7dliUHl1iX7UjfQXxivIKphKhc5SoEJpOJgoKb35gMBgMqlcpmoYSoLO0D2pNTkMPMfTM5m3OWHek7iPoxSg4ZCbtSrnME7du3Z+rUqXTu3BmArVu3ypVCokboENCB9a1fIvTiMoL+E4KDSc2Gvw4hKKCD0tGEqDTlKgTPPfcctWvX5sCBA5jNZjp16kSXLl3K3MZsNrN06VJSU1PR6XRERUXh7198YhCDwcC0adOIiooiMDCwXNsIUZGcL6yh29VV3BiQtL7WTODFpVytFYreL1LZcEJUknIVgilTpvDOO+9YewTlsWfPHgoLC4mNjSUlJYWEhARiYmKs60+ePMmSJUu4dOlSubcRoqK5n3oftVlfbJnarMf91PtSCITdKNc5guvXr5Ofn39XO05OTqZVq1YAhIaGcvLkyWLrCwsLGTVqFIGBgeXeRoiKpjGk39VyIWqicvUInJycGDJkCEFBQcXGGBozZsxtt9Hr9bi4uFgfq9VqTCYTGo0GKBq24m63ucHHx6c8sUul1Wrva/uaRNoCcjSeeJqulFiebtJQ4Ggg0D2wlK1qPnlvFFfT2+OOhSAtLY3WrVsTERFxVzeQOTs7o9ff7HJbLJYSH+j3uk1WVla5c/yZj4/PfW1fk0hbgHPoVMzHYoodHjKpHJmQbWH/18/wzVPf4Kx1VjChMuS9UVxNaI+AgNsPpVLmoaGtW7fy7rvvsmvXLjZs2ICXlxedOnWy/itLWFgY+/fvByAlJYWgoKA7Br2XbYS4H3q/SK6GTcfoGIgFFUbHQHKazOTRvy4kMSuRmF9iZP4CUeOV2SP4/vvvmTVrFt7e3qSkpLBixQrrMfw7adOmDUlJSUyYMAGLxUJ0dDTbt28nPz+fbt26lXsbIWxN7xeJ3i+y2Le+J4DRD41mxr4ZNPNuxuCIwcqGFMKG7nho6MbhoNDQUHJycsq9Y7VazaBBg4otu/XE8A2TJk0qcxshlDLsgWEczT7K1ylfM6D5AJy0MgeHqJnKLAR/vnv4Tsf4hahJVCoVHz32EYXmQikCokYr9+T1ULIwCFHTuehc8HT0JN+YT9yeOK4ariodSYgKV2aPIDU1lf79+1sfGwwG+vfvj8ViQaVSsWzZMpsHFKIqSL6czMLEhRzOOsyyJ5ahUUvvWNQcZRaCefPmVVYOIaq0VnVaMa3DNMZuH0vcnjgmPDxB6UhCVJgyC0GdOnUqK4cQVd4rTV/hyKUjLEhaQFPvpjzb+FmlIwlRIe7qHIEQ9m5K+ym0q9uOab9NQ2/U33kDIaoBmXxYiLugU+tY3G0x2fnZdnnHsaiZpEcgxF3ydvImxCsEi8XC6uOryTfe3YCMQlQ1UgiEuEeJWYkM2zaMsdvHyjAUolqTQiDEPWpVpxUjHhzBv47/iyWHltx5AyGqKCkEQtyH4Q8Op2fDnkz9bSo//fGT0nGEuCdSCIS4D2qVmtmdZhNWK4whW4aQW5CrdCQh7ppcNSTEfXLVufLZ459x4soJ3BzclI4jxF2THoEQFaC+e3061y+a03v/xf2YzCaFEwlRflIIhKhAhy8dpvd3vZm+d7rSUYQoNykEQlSgZt7N6BvWl48TP2btibVKxxGiXKQQCFGBVCoV09pP42H/hxn580iSMpOUjiTEHUkhEKKCOWgcWNxtMbWdazNw80CuGK4oHUmIMslVQ0LYgI+zD58+/ik703fi6eCpdBwhyiSFQAgbCa8dTnjtcAAu5F3A19lXZvkTVZIcGhLCxk5cOUGnf3Xis8OfKR1FiFJJIRDCxoI9g2lXtx2Tfp3EL+d+UTqOECVIIRDCxtQqNXM7zSXEK4SoH6NIzUlVOpIQxUghEKISuDm48enjnwIw4IcBMiaRqFKkEAhRSRp6NGRB1wVE1IlAq5brNETVIe9GISrRo4GP8mjgowAUmApw0DgonEgI6REIoYhzuefo9k031p1ap3QUIaQQCKEEH2cfajnVYvhPwzl06RAAzhfW4LurDXW31cN3VxucL6xROKWwF1IIhFCAo8aRJd2W4OXoxT9++AfGswl4HotBaziHCgtawzk8j8VIMRCVQgqBEArxdfHl0+6fkqXPQnf8XdRmfbH1arMe91PvK5RO2BMpBEIoqGWdlsx6dBZ+qoJS12sM5yo5kbBHUgiEUNjfQ/6OXudT6ro8benLhahINrt81Gw2s3TpUlJTU9HpdERFReHv729dv3fvXr755hvUajWdO3emW7duAMTExODi4gKAr68v0dHRtoooRJVR0PhdHJNHorXc7BkYVQ4UNH5XwVTCXtisEOzZs4fCwkJiY2NJSUkhISGBmJgYAIxGI8uWLSMuLg4nJycmTpxI69atrQVg0qRJtoolRJWk94skz6jHfHQsgRozGWYt6ibvYfGLVDqasAM2OzSUnJxMq1atAAgNDeXkyZPWdefOncPf3x83Nze0Wi1hYWEcPXqU1NRUDAYD06ZNY/LkyaSkpNgqnhBVzk5VQ1pleNHuWhvqnTLywLYZ7LuwT+lYwg7YrEeg1+ut3/AB1Go1JpMJjUZTYp2zszN5eXk4Ojry1FNP0bVrVzIyMoiLi2P27NloNJpi+/bxuffjplqt9r62r0mkLYpTsj22pW4jeks0KyJX0KlBJ+btmUfMjzEM3TaUo4OPolFr7ryTCiTvjeJqenvYrBA4Ozuj19+8HM5isVg/0J2dncnPz7eu0+v1uLq6UrduXfz9/VGpVAQEBODm5sbly5dL/AKysrLuOZePj899bV+TSFsUp2R7/HziZ+K7xBPuGk5WVhZ9/9IXn+4+/JrxK5ezL5NXmEemPpMGHg0qJY+8N4qrCe0REBBw23U2OzQUFhbG/v37AUhJSSEoKMi6LjAwkIyMDHJzczEajRw9epTQ0FC2bt1KQkICANnZ2ej1emrVqmWriEJUGdER0XQI6FBsWfcG3ZnYdiIA7+99n+5rurPy2EosFosSEUUNprLY6F1146qhtLQ0LBYL0dHRnD59mvz8fLp162a9ashsNtO5c2d69OiB0Whk/vz5ZGVloVKpePnllwkLCyux7/T09HvOVRMqe0WRtiiuKrfHudxzDNs2jF0Zu+j5l5580PEDvJ28bfZ6VbktlFAT2qOsHoHNCoEtSSGoGNIWxVX19jCZTSw6uIjpe6dT26k2i7otorVfa5u8VlVvi8pWE9pDkUNDQoiKpVFriI6IZv3T66nrWhcf55p78lJULikEQlQz4T7hrHt6HQ09GmKxWJj22zSSs5OVjiWqMSkEQlRDKpUKgPTr6fzr+L/oubYnSw8txWwxK5xMVEdSCISoxgLdAvnx2R/pGNCRd3e9yyubXuFC3gWlY4lqRgqBENWcj7MPy55YRlyHOH7N+JW+G/tKz0DcFZmzWIgaQKVS8WqzV2kf0J4sfRZqlRqj2YjBZMBV56p0PFHFSY9AiBokxCuEtnXbAvDxgY95fM3j7L+4X+FUoqqTQiBEDdW2blsKzYU8/d3TfPT7RxjNRqUjiSpKCoEQNVTbum3ZHLmZ3sG9mblvJs+uf5bUnFSlY4kqSAqBEDWYp6MnH3f5mHmd53Hyykku5V9SOpKogqQQCGEHIkMi2d13Nw/6PgjAqpRVXDFcUTiVqCqkEAhhJ1x0RXOAnL12lphfYuj2TTd2pO9QOJWoCqQQCGFn6rvX57ve3+GidaHPhj5M+20aBpMBAOcLa/Dd1Qbdaid8d7XB+cIahdOKyiCFQAg71LJOSzY9s4l+TfuxIGkBfTb0wen8GjyPxaA1nEOFBa3hHJ7HYqQY2AG5oUwIO+Wic+H9ju/TpX4XcgtzUaWMR23WF3uO2qyHY+PBL1KhlKIySCEQws493uBxALz+eLPU9V7ma+hLXSNqCjk0JIQAwORY+sQlt1suag4pBEIIAK4Fj8Wsdi627LoZFhluzjdeDSc0FOUghUAIAYDeL5Ldvq9x1qjGAqQZ1axy7IZ/2EgA0nLS6LCyA+/tfo+DWQelKNQgUgiEEADsSN/Bk3u+5NdmX1P4nIHfmn3NqFO/W4e0zi3MJdgzmEVJi+jxbQ86rupI3J44svTVey5fIYVACPE/iZmJLOy6kA4BHQDoENCBhV0XkpiZCECz2s1Y/rfl7O+3n5mPzKSBewOWHlyKWlX0MfL7xd85fOmw9BSqIZWlGv7W0tPT73lbHx8fsrLkGwxIW/yZtMdN5W2L3IJc3BzcAOizoQ/b07cT7BnMU8FP0esvvWjq3dQ6rWZ1VhPeGwEBtz/pLz0CIcQ9u1EEAOK7xPNBxw8IcA1g3oF5dF/TnTe2vqFgOlFech+BEKJC1HauTb+m/ejXtB9Z+iw2nt6Ij7MPAFcNV3l+w/P0aNCDXsG9CK0VqnBacSspBEKICufj7MOrzV61Ps7UZ+Lh4MGHv3/IrN9nEVYrjKeCn+LlJi/j6+KrYFIBcmhICFEJQrxCWN1rNfte3se09tOo5ViLWftmkVOQA8Dxy8c5ceUEAPGJ8SVGRd2RvoP4xPhKz20vpEcghKg0fi5+DGg+gAHNB5CZl0kdlzoAzN4/m7Un19LUuyktfVpy/viH9PR3xcV4iTxtbTZf0NO93WcKp6+5pBAIIRRxowgATHh4Ag/4PsC6U+tQp69kjh+4GotGOHI1ZjG/to7zplNYLO1rxFVIVY1cPmrHpC2Kk/a4Scm2qL3zIRwLzpdYfqYQItI9CHIPooFHA3r9pRe9G/XGYrHwR+4fBLgGoFFrbJKpJrw3yrp8VHoEQogqxaHgQqnLG+iKptxMzUnlaPZRwmuHA3DZcJm2X7dFq9JSz70eDT0aEuQexDMhz9DGvw1GsxGDyYCrzrUyf4xqRQqBEKJKydPWxtVY8tt3ntaH2A6xJZbr1DpmPjKTM9fOkJqTSlpOGgcyD9DSpyVt/NtwNPsoPb7tgY+zD0HuQcUKRYhXCCazCZVKZb1D+ob4xHh6OebQ+vIaNIZ0fB0D2FsrkvUGD6Ijom328ytBCoEQokrZ4NCBSNP3aC0F1mVGlQMbHDrQsZTnuzu407dJ3xLLb4yR5O3kzdi/jiUtJ40zOWfYfX43a0+upY1/G0K8Qvhv2n+J3hJNkHsQQR5Fh50aujeki/o8LdM/Qfu/+qA1nKNl+jwIGGqLH/u2/lyQTDYoSDY7R2A2m1m6dCmpqanodDqioqLw9/e3rt+7dy/ffPMNarWazp07061btztuc4OcI6gY0hbFSXvcpHRbOF9Yg/up960ffNeCx6KvwFnSCkwFqFQqdGodhy8dZvXx1aTlpJF6LZUzOWfQG/Xom/nhVFjyMFWG2YG++e34uMvHeDt5s/XsVnam78RF54KrzhU3nRuuOld6/qUnOrWOjOsZXC+8jou2aL2LzgWdWlfurGnJ79MyfR4ut3RY8syQFDCUoCZjy70fRc4R7Nmzh8LCQmJjY0lJSSEhIYGYmBgAjEYjy5YtIy4uDicnJyZOnEjr1q05duzYbbcRQtgPvV9khX7w/5mDxsH6/+a1m9O8dnPrY4vFQpY+C8fdD5S6rZ+qgJzCHLTqoo/PA5kH+OTwJxhMhmLPOzXwFFD0jf7Tw58WW+esdeb4/x1HpVIxc99Mtv2xDVetK666on/eTt5MbjcZgPDMr4oVAQAXNbS+vIaLlL8QlMVmhSA5OZlWrVoBEBoaysmTJ63rzp07h7+/P25uReOUhIWFcfToUVJSUm67jRBCVAaVSkUdlzqYHAPQGs6VWG92CmR95/XWx8MfHM7wB4dTaC4krzCP64XXyTPm4aAuKjZ9wvrQ2q811wuvW/8Vmgutl8G669zx0Hlw3XidS/mXyCvMw0XnYt2/i/ESlHLFrMZw70dG/sxmhUCv1+PicvOHUavVmEwmNBpNiXXOzs7k5eWVuc2tfHx87jmXVqu9r+1rEmmL4qQ9bpK2gCMB/Qg5+UGJQzInAvrR5C7appNPJzrR6bbrx3cZz3jG33a9xaUe6P8oucKlfoX9jmxWCJydndHrb055bbFYrB/ozs7O5OfnW9fp9XpcXV3L3OZW93PsUuljn1WJtEVx0h43SVvAqitaegUMLXaSNqlWJOuvaImuxLZJ83yGltdLOUfg8TRBd5FDkXMEYWFh7Nu3j/bt25OSkkJQ0M15TwMDA8nIyCA3NxcnJyeOHj1K7969AW67jRBCVKYbV+RcZKy1MAYBlX3h6HqDB5RWkAweFZbF5lcNpaWlYbFYiI6O5vTp0+Tn59OtWzfrVUNms5nOnTvTo0ePUrcJDAwssW+5aqhiSFsUJ+1xk7RFcTWhPcrqEcgQE3ZM2qI4aY+bpC2KqwntITOUCSGEuC0pBEIIYeekEAghhJ2TQiCEEHauWp4sFkIIUXGkRyCEEHZOCoEQQtg5KQRCCGHn7GJimvLOc2AvjEYjCxYsIDMzk8LCQp599llat26tdCxFXb16lbFjxzJhwoRS72a3J99++y179+7FaDTyxBNP0KVLF6UjKcZoNDJ//nwyMzNRq9X885//rJHvD7voEdw6N8JLL71EQkKC0pEU9csvv+Du7s6UKVN4++23+eSTT5SOpCij0cjixYtxcHC485NruMOHD3Ps2DGmTp3K5MmTq/3dtPdr//79mEwmpk2bxnPPPceKFSuUjmQTdtEjKGtuBHvUrl072rZta31c2giv9uSLL76ge/furF27VukoiktMTCQoKIiZM2ei1+vp16+f0pEUVbduXcxmM2azmby8PLTamvmRaRc9gtvNc2CvnJycrEN+f/jhh7z44otKR1LMtm3b8PDwsH5RsHc5OTmcOnWKESNG8PrrrzN37lzs+QpzJycnMjMzGT58OIsWLeJvf/ub0pFswi4KQXnnObAnWVlZTJ48mUceeYSOHUubEtw+bN26laSkJCZNmsSZM2f4+OOPuXLlitKxFOPu7k5ERARarZaAgAAcHBzIyclROpZiNmzYQEREBHPmzGHGjBnMnz+fgoICpWNVuJrZz/mTsuZGsEdXrlwhNjaWgQMH0qJFC6XjKGry5MnW/0+aNInXX38dLy8vBRMpq0mTJmzcuJFevXpx+fJl8vPzcXd3VzqWYlxdXa2Hg9zc3DCZTJjNZoVTVTy7uLO4vPMc2IvPPvuMnTt3FmuDt99+2+5Plt4oBPb83gBYvnw5hw8fxmw207dvX7s+bJafn098fDxXrlzBaDTSs2fPGtmDtotCIIQQ4vbs4hyBEEKI25NCIIQQdk4KgRBC2DkpBEIIYeekEAghhJ2zi/sIRPVw8eJF3njjDaKioooNdPbdd99x9uxZhgwZUiGvM2TIEEaMGEGjRo0qZH9lycvLIy4ujuvXr9OnTx8efvhh67r58+eTlJSEh4cHUDTmUdOmTXnhhRfu616G33//nePHj9OnTx+2bdvGr7/+ytixY+/7ZxE1lxQCUaWoVCoSEhJo0qQJAQEBSse5b2fOnOHKlSvMmzev1PVPPvkkvXv3BorueP/222+JjY3lgw8+QK2+tw77iRMnyM3NvefMwv5IIRBVioODA7169WLOnDnExsaWGORr/vz51K9f3/rheevjIUOG0KFDBw4fPkxubi5PP/00ycnJnD59Go1GQ0xMDN7e3gD85z//ITU1lcLCQnr16mXtgezdu5c1a9ZgNBpxdHTklVdeITQ0lFWrVnH8+HGys7Np0KABb775ZrFcu3fvZvXq1VgsFpycnOjfvz8uLi4sWLCA7OxsRo8eTWxsbJk37alUKiIjI/npp59ISkqiVatWHDt2jC+//BKDwYBarea5557joYceIj8/n6VLl5KRkUFubi5OTk4MGzaM69evs3nzZsxmMy4uLtStW5fLly8TFxdHVlYWGo2GN998k3r16vHbb7+xZs0aVCoVarWafv360axZs4r8dYpqQgqBqHIiIyM5ePAgX331Fa+++updbXtjuPGdO3cyZ84cPvjgAxo2bMiMGTPYtm0bkZGRQFHB+eCDD8jOzmbMmDE0btwYrVbLihUrmDRpEu7u7pw9e5apU6cyd+5cADIzM5k1a1aJcarOnTvHkiVLmDZtGn5+fhw6dIjp06cze/ZsoqKi+PTTT5kxY0a5f4YGDRqQlpZGSEgI8fHxjB8/Hl9fX7Kzsxk/fjwNGjTgxIkTuLi4EBsbC8DixYvZtGkTAwcOpHv37ly7do2+ffuybds2Ll68yPDhw/H39+ezzz5j3bp1DB48mOXLlzN06FBCQ0NJTEzkyJEjUgjslBQCUeWo1WqGDh1KTEzMXQ9vcOMYvJ+fH15eXjRs2ND6+NbDJd26dQPA29ubli1bcvDgQTQaDVeuXGHKlCnW56lUKs6fPw9A48aNSx2s8NChQ7Ro0QI/Pz8AwsPD8fT05NSpU6hUqrvKf4ODgwMpKSlcuXKlRBFJTU2lbdu2+Pr68v3333P+/HmOHDlCaGhoqfsKCQmxTsTUsGFDfvvtNwDat2/PzJkzefDBB2nZsqW1lyXsjxQCUSX5+PgwaNAg5s+fz6OPPmpd/ucPVqPRWOyxTqez/r+sseNvPf5usVjQarWYTCbCw8MZPny4dV1WVhbe3t7s3r0bJyenUvdV2iBkZrMZk8l01+PXWywWTp8+TY8ePcjLyyMwMJD33nvPuj47OxsPDw9++OEH/vvf/9KjRw86duyIm5sbFy9eLHWftxavW9uvb9++dO7cmaSkJLZt28a6deuIi4u7q7yiZpDLR0WV1bZtW1q1asXGjRutyzw8PKwTC2VnZ3PkyJF72ve2bduAog/6gwcPEh4eTosWLUhKSuLcuXNA0dU3o0ePvuOwwy1atCAxMZELFy4ART2ES5cu0bhx47vKZDabWb16Ne7u7jRr1ozQ0FAyMjKsP+OZM2cYNmwY2dnZHDhwgE6dOtGlSxcCAgLYt2+ftSBpNJo7zrdhMpkYMmQIBQUFPP7447z22mukpaVRWFh4V5lFzSA9AlGlDRgwgOTkZOvjHj16MG/ePIYNG4avry/h4eH3tN/CwkLGjBmD0WhkwIAB1iuUBg0axOzZs4GiXkNMTMxtewI31KtXj9dee42ZM2diNptxcHBgzJgxxSZDup0NGzbwyy+/oFKpMJvNNGrUiHHjxgFFRW/kyJEsX76cwsJCzGYzb7zxBr6+vvTu3ZtFixaxZcsWLBYLoaGhpKWlAUWHpubOncunn35KcHBwqa+r0Wjo378/c+bMQavVolKpGDx4cLEelbAfMvqoEELYOTk0JIQQdk4KgRBC2DkpBEIIYeekEAghhJ2TQiCEEHZOCoEQQtg5KQRCCGHnpBAIIYSd+38yDI15PoVTmwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "\n",
    "plt.plot(i, pois_mle,'x--', color='green',label='Mixture of Poissons')\n",
    "plt.plot(np.arange(freq.size), freq/np.sum(freq), 'o', color='orange',label='Original data')\n",
    "plt.legend()\n",
    "plt.xlabel(\"Number of Deaths\")\n",
    "plt.ylabel(\"Proportion\");\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***Comparison of methods:<br>***\n",
    "\n",
    "***The gradient ascent optimization method has 5 steps: <br>***\n",
    "1. First calculate log likelihood\n",
    "2. Calculate gradient for all the 3 parameters. \n",
    "3. Then update the parameters and check condition if new likelihood is greater than the previous one. \n",
    "3a. We also will be checking for Nan values since gradient may be negative.<br>\n",
    "4. The EM algorithm involves calculating Ez, $\\alpha$ $\\mu1$ and $\\mu2$ and updating the value with each iteration until the increase in likelihood is less than a particular threshold.<br>\n",
    "\n",
    "Hence we can see that the EM algorithm requires less steps as compared to the gradient ascent method.<br>\n",
    "\n",
    "* Also, we see that the running time of poisson gradient is more than that of EM method. \n",
    "* (Here we observe a huge difference because we are running the posisson algorithm for 1000 times.Even without 1000 iterations that it runtime is more than EM Algorithm).\n",
    "* The gradient ascent method is slower since it involves computing the gradient at each iteration which is computationally expensive. <br>\n",
    "\n",
    "***Advantages of each method:***\n",
    "\n",
    "1. Loglikehood function has a close form solution and hence can converge faster than EM Method.\n",
    "    * This comes with caveat of choosing good starting values\n",
    "2. The EM method is more stable and tends to converge even if the starting values are off target."
   ]
  }
 ],
 "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.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}