MLE_tutorial.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy import optimize\n",
    "from scipy import stats\n",
    "import matplotlib.pyplot as plt\n",
    "import math\n",
    "from scipy.special import factorial"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# MLE for normally distributed data\n",
    "\n",
    "Let's walk through an example of using MLE to estimate the parameters for the distribution of some sample data. To start, let's [simulate some data from a normal distribution](https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.norm.html) with known parameters, which we will refer to as $\\mathbf{x}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "MU = 4\n",
    "SIGMA = 3\n",
    "N = 1000\n",
    "norm_sample = stats.norm.rvs(loc = MU, scale = SIGMA, size = N)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Visualize\n",
    "Let's take a look at what our sample data looks like."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAEWCAYAAABollyxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAGBNJREFUeJzt3X+8bXVd5/HXO0BRIIHuBQHJm0pOaA06dwh/TEORhGhhPTQxSzQd8qE2OWMm2g/RdMIardEpG0wUf0OWSYKjhJlpiV5JEQaUi17lcq/3XkB++Rv8zB/re2xz2OeeH/ucs8/99no+Hvux117ru9f6rB/7fdb+rr33SVUhSerX9027AEnSyjLoJalzBr0kdc6gl6TOGfSS1DmDXpI6Z9CvgCRvTvKKNvyfknxuGef9/iSnteGnJ/noMs77qUk+uFzzW8RyH5XkmiS3J3nCai9/KUb3w54mSSV50LTrWA1JzkzytmnXMW0G/Qqrqn+sqgfP126hB2RVPbaqzp20riQb2gt+75F5v72qTpx03kvwcuB/V9X+VfU34xokOTXJpUm+lmRnG35OkqxyrcDy7QdpNRj0e4gMet1f9weunGtikhcA/wv4I+C+wKHAs4FHAfdYjQKnZfQPsbRUvQbHqkrysCSXJbktyXnAviPTjk+ydeTxi5Jc39p+LskJSU4CXgI8uXVffKa1/XCSVyb5GPB14AFt3LPuuvi8LsktSa5OcsLIhC1Jfnrk8ei7ho+0+5vbMh8xuysoySOTfLLN+5NJHjky7cNJfj/Jx9q6fDDJut1so/+SZHOSm5JckOTwNv5a4AHA37Y67jnrefdhOON/TlW9u6puq8G/VNVTq+pbrd3jkvxLkluTXJfkzLn2wextk+TYJJvac3ckeU0bv2+StyW5McnNbRscOrL+z2rDD0zyodbuhiRvT3LgrGX9ZpLL27Y8L8m+jNH2wceS/HGSm4AzJ51/khcm2Z5kW5Jfnb19k7wlya4kX0ryOzMnFLNquTnJF9ox8fS2jXdmN91Xrd0X2vHxxSRPXcT2emFbn68leWOSQzN0l92W5O+SHNTazrwzPb2t3/YMJwZz1XRckn9q6/OZJMfP1bYrVeVtghvDGeWXgP8G7AM8EfgO8Io2/Xhgaxt+MHAdcHh7vAF4YBs+E3jbrHl/GPgy8BBg7zb/DwPPatOfDtwxsuwnA7cAB7fpW4CfHpnf95bRll3A3iPTnw58tA0fDHwV+JW27Ke0xz8wUtu1wA8D92qPz5pjG/0UcAPwcOCewOuAj4xMv0uds557UlvHvcdNH2l3PPCjDCcvPwbsAJ4wex+MWybwz8CvtOH9gePa8K8BfwvcG9gL+A/A94+s/8x+eBDwmLZu6xn+iP7JrGV9Aji8bdergGfPsR4z+/TX23a/1yTzb9tvB/BQYD/gHW2/P6hNfwvwXuCAdkx8HnjmrFqe0db/FQzH45+2Wk4EbgP2H7Me+wG3Ag9ujw8DHrKI7fVxhnduRwA7gcuAh7XnfAh46azj+J1tmT8K7BrZt2fyr8f8EcCNwMkMx8lj2uP1086Rlb55Rj+54xhC9k+q6jtV9W7gk3O0vZPhQD06yT5VtaWqrp1n/m+uqiur6o6q+s6Y6TtHln0e8DngcUtcl1GPA66pqre2Zb8TuBr42ZE2b6qqz1fVN4DzgWPmmNdTgXOq6rIazsBfDDwiyYYF1LEOuKGq7pgZMXJG9o0kPwFQVR+uqs9W1Xer6nKGF/5/XuC6fgd4UJJ1VXV7VX18ZPwPMITinVX1qaq6dfaTq2pzVV1cVd+qql3Aa8Ys+7VVta2qbmL44zHXtgLYVlWva9v9GxPO/xcZ9tMVVfU1huADIMleDCcHL67hndIW4NUMf9xnfLGq3lRVdwLnAUcCL2+1fBD4NkNwj/Nd4KFJ7lVV26vqykVsr9dV1Y6quh74R+DSGt7FfQt4D0Poj3pZVX2tqj4LvInhxGS2XwYuqqqL2nFyMbCJIfi7ZtBP7nDg+mqnDM2XxjWsqs3A8xlebDuTvGumC2M3rptn+rhlzzfPhTicu6/HlxjOimZ8ZWT46wxnw/POq6puZziTOmKO9qNuBNblrheNH1lVB7ZpM90MP57k71sXxC0MffhzdiXN8kyGdyZXt+6Zx7fxbwU+ALyrdQv8YZJ9Zj85ySFtX16f5FbgbWOWvdBtBbP2+YTzP3zW/Eb36Tr+9R3p6PTR/bJjZPgbAFU1e9zd1qX9UXkyw37YnuTCJP9uEeszexnzLXP2Oo57DdwfeFI7Sbg5yc3AoxnebXTNoJ/cduCI5C6f/vjBuRpX1Tuq6tEMB10Br5qZNNdT5ln+uGVva8NfY+h2mHHfRcx3W6tx1A8C18/zvHnnlWQ/hjPlhczrn4FvAafM0+4dwAXAkVV1H+DPgZntcpft0M5k1888rqprquopwCEM++PdSfZr75JeVlVHA48EHg88bcyy/4Bhe/5YVX0/w5njJJ8Gmr1vJpn/doaz8Bmjx+YNDO9a7j9r+lL28d1U1Qeq6jEMQXo18IY2abm3F9x9HbeNaXMd8NaqOnDktl9VnTXhstc8g35y/8zQj/lfk+yd5BeAY8c1TPLgJD+V4YLjNxnOTO5sk3cAG7L4T9Yc0pa9T5InAT8CXNSmfRo4tU3byHD9YMYuhrfWD5hjvhcBP5zkl9p6PRk4GnjfIuuDIYSfkeSYtu7/g+Gt+Jb5nlhVNwMvA/4syROT7J/k+5Icw9AnO+MA4Kaq+maSY4FfGpn2eWDfDBds9wF+h6ELDYAkv5xkfVV9F7i5jb4zyU8m+dH2h+FWhlC8k7s7ALid4cL2EcALF7BNFmOS+Z8PPD3J0UnuDbx0ZkLrjjkfeGWSA5LcH/jvDGfYE2kXT3+u/VH/Vqt/ZtutxPb63ST3TvIQhmsK541p8zbgZ5P8TJK9MlxsPz7J/ZZh+WuaQT+hqvo28AsMF66+yvB29a/naH5P4CyGM6mvMIT0S9q0v2z3Nya5bBElXAoc1eb5SuCJVXVjm/a7wANbXS9jCNyZur/e2n+svY09btZ63chwBvsChi6S3wIeX1U3LKK2mXld0mr5K4YzzAcCpy7i+X/IEEC/xXBNYgfwf4AXAf/Umj0HeHmS24DfYwiwmeff0qb/BcPZ6teA0U/hnARcmeR2ho9xnlpV32R4B/RuhpC/CvgHxofgyxguNN8CXMjc+3+pljz/qno/8CcMFzA3t/tRv86wPb4AfJThGDln8pL5PoZjZxtwE0Mf/HPatJXYXv/AsH6XAP+zXT+4i6q6juGd4UsYTnSuY/gj030O5q7du5K052gX9L8I7DN6wV531f1fMkn6t86gl6TO2XUjSZ3zjF6SOrcmfjBp3bp1tWHDhmmXIUl7lE996lM3VNX6+dqtiaDfsGEDmzZtmnYZkrRHSTL2W/iz2XUjSZ0z6CWpcwa9JHXOoJekzhn0ktQ5g16SOmfQS1LnDHpJ6pxBL0mdWxPfjJXWsg1nXDiV5W45azn+x7vkGb0kdc8zeu0RpnVWLfXAM3pJ6pxBL0mdM+glqXMGvSR1zqCXpM4Z9JLUOYNekjpn0EtS5wx6SeqcQS9JnTPoJalzBr0kdc6gl6TOGfSS1DmDXpI6Z9BLUucMeknqnEEvSZ0z6CWpc/MGfZIjk/x9kquSXJnkN9r4g5NcnOSadn9QG58kr02yOcnlSR6+0ishSZrbQs7o7wBeUFU/AhwHPDfJ0cAZwCVVdRRwSXsM8FjgqHY7HXj9slctSVqweYO+qrZX1WVt+DbgKuAI4BTg3NbsXOAJbfgU4C01+DhwYJLDlr1ySdKCLKqPPskG4GHApcChVbUdhj8GwCGt2RHAdSNP29rGzZ7X6Uk2Jdm0a9euxVcuSVqQBQd9kv2BvwKeX1W37q7pmHF1txFVZ1fVxqrauH79+oWWIUlapAUFfZJ9GEL+7VX11230jpkumXa/s43fChw58vT7AduWp1xJ0mIt5FM3Ad4IXFVVrxmZdAFwWhs+DXjvyPintU/fHAfcMtPFI0lafXsvoM2jgF8BPpvk023cS4CzgPOTPBP4MvCkNu0i4GRgM/B14BnLWrEkaVHmDfqq+ijj+90BThjTvoDnTliXJGmZ+M1YSeqcQS9JnTPoJalzBr0kdc6gl6TOGfSS1DmDXpI6Z9BLUucW8s1Y6Xs2nHHhtEuQtEie0UtS5wx6SeqcQS9JnTPoJalzBr0kdc6gl6TOGfSS1DmDXpI6Z9BLUucMeknqnEEvSZ0z6CWpcwa9JHXOoJekzhn0ktQ5g16SOmfQS1LnDHpJ6pxBL0mdM+glqXMGvSR1zqCXpM4Z9JLUOYNekjpn0EtS5wx6SeqcQS9JnZs36JOck2RnkitGxp2Z5Pokn263k0emvTjJ5iSfS/IzK1W4JGlhFnJG/2bgpDHj/7iqjmm3iwCSHA2cCjykPefPkuy1XMVKkhZv3qCvqo8ANy1wfqcA76qqb1XVF4HNwLET1CdJmtAkffTPS3J569o5qI07ArhupM3WNk6SNCVLDfrXAw8EjgG2A69u4zOmbY2bQZLTk2xKsmnXrl1LLEOSNJ8lBX1V7aiqO6vqu8Ab+Nfuma3AkSNN7wdsm2MeZ1fVxqrauH79+qWUIUlagCUFfZLDRh7+PDDziZwLgFOT3DPJDwFHAZ+YrERJ0iT2nq9BkncCxwPrkmwFXgocn+QYhm6ZLcCvAVTVlUnOB/4fcAfw3Kq6c2VKlyQtxLxBX1VPGTP6jbtp/0rglZMUJUlaPn4zVpI6Z9BLUucMeknqnEEvSZ0z6CWpcwa9JHXOoJekzhn0ktQ5g16SOmfQS1Ln5v0JBEnTseGMC6ey3C1nPW4qy9XK8Yxekjpn0EtS5wx6SeqcQS9JnTPoJalzBr0kdc6gl6TOGfSS1DmDXpI6Z9BLUucMeknqnEEvSZ0z6CWpcwa9JHXOoJekzhn0ktQ5g16SOmfQS1LnDHpJ6pxBL0mdM+glqXN7T7sALd6GMy6cdgmS9iCe0UtS5wx6SeqcQS9JnTPoJalz8wZ9knOS7Exyxci4g5NcnOSadn9QG58kr02yOcnlSR6+ksVLkua3kDP6NwMnzRp3BnBJVR0FXNIeAzwWOKrdTgdevzxlSpKWat6gr6qPADfNGn0KcG4bPhd4wsj4t9Tg48CBSQ5brmIlSYu31D76Q6tqO0C7P6SNPwK4bqTd1jbubpKcnmRTkk27du1aYhmSpPks98XYjBlX4xpW1dlVtbGqNq5fv36Zy5AkzVhq0O+Y6ZJp9zvb+K3AkSPt7gdsW3p5kqRJLTXoLwBOa8OnAe8dGf+09umb44BbZrp4JEnTMe9v3SR5J3A8sC7JVuClwFnA+UmeCXwZeFJrfhFwMrAZ+DrwjBWoWZK0CPMGfVU9ZY5JJ4xpW8BzJy1KkrR8/GasJHXOoJekzhn0ktQ5g16SOmfQS1LnDHpJ6pxBL0mdM+glqXMGvSR1zqCXpM4Z9JLUOYNekjpn0EtS5wx6SeqcQS9JnTPoJalzBr0kdc6gl6TOGfSS1DmDXpI6Z9BLUucMeknqnEEvSZ0z6CWpcwa9JHXOoJekzhn0ktQ5g16SOmfQS1LnDHpJ6pxBL0mdM+glqXMGvSR1zqCXpM7tPe0CJK0tG864cGrL3nLW46a27J55Ri9JnZvojD7JFuA24E7gjqramORg4DxgA7AF+MWq+upkZUqSlmo5zuh/sqqOqaqN7fEZwCVVdRRwSXssSZqSlei6OQU4tw2fCzxhBZYhSVqgSYO+gA8m+VSS09u4Q6tqO0C7P2TcE5OcnmRTkk27du2asAxJ0lwm/dTNo6pqW5JDgIuTXL3QJ1bV2cDZABs3bqwJ65AkzWGiM/qq2tbudwLvAY4FdiQ5DKDd75y0SEnS0i056JPsl+SAmWHgROAK4ALgtNbsNOC9kxYpSVq6SbpuDgXek2RmPu+oqv+b5JPA+UmeCXwZeNLkZUqSlmrJQV9VXwD+/ZjxNwInTFKUJGn5+M1YSeqcQS9JnTPoJalzBr0kdc6gl6TO+Xv0E5jm73ZL0kJ5Ri9JnTPoJalzBr0kdc6gl6TOGfSS1DmDXpI6Z9BLUucMeknqnEEvSZ0z6CWpcwa9JHXOoJekzhn0ktQ5g16SOmfQS1LnDHpJ6pxBL0mdM+glqXMGvSR1zqCXpM4Z9JLUOYNekjpn0EtS5wx6Serc3tMuYFIbzrhw2iVI0pq2xwe9pH5M68Rty1mPm8pyV4tdN5LUOYNekjpn0EtS5wx6SeqcF2Ml/Zs3zU/vrcaF4BU7o09yUpLPJdmc5IyVWo4kafdWJOiT7AX8KfBY4GjgKUmOXollSZJ2b6XO6I8FNlfVF6rq28C7gFNWaFmSpN1YqT76I4DrRh5vBX58tEGS04HT28Pbk3xugfNeB9wwcYUrZ63XB9a4HNZ6fbD2a1zr9cEq1JhXTfT0+y+k0UoFfcaMq7s8qDobOHvRM042VdXGpRa20tZ6fWCNy2Gt1wdrv8a1Xh/sGTUuxEp13WwFjhx5fD9g2wotS5K0GysV9J8EjkryQ0nuAZwKXLBCy5Ik7caKdN1U1R1Jngd8ANgLOKeqrlym2S+6u2eVrfX6wBqXw1qvD9Z+jWu9PtgzapxXqmr+VpKkPZY/gSBJnTPoJalzaz7ok5yZ5Pokn263k+doN5WfXEjyR0muTnJ5kvckOXCOdluSfLatw6ZVqm232yTJPZOc16ZfmmTDatTVln1kkr9PclWSK5P8xpg2xye5ZWTf/95q1TdSw273Wwavbdvw8iQPX+X6HjyyfT6d5NYkz5/VZtW3Y5JzkuxMcsXIuIOTXJzkmnZ/0BzPPa21uSbJaatY35p9LU+sqtb0DTgT+M152uwFXAs8ALgH8Bng6FWq70Rg7zb8KuBVc7TbAqxbxe027zYBngP8eRs+FThvFes7DHh4Gz4A+PyY+o4H3jetY28h+w04GXg/w3dHjgMunWKtewFfAe4/7e0I/ATwcOCKkXF/CJzRhs8Y91oBDga+0O4PasMHrVJ9a/K1vBy3NX9Gv0BT+8mFqvpgVd3RHn6c4TsDa8FCtskpwLlt+N3ACUnGfdlt2VXV9qq6rA3fBlzF8I3qPc0pwFtq8HHgwCSHTamWE4Brq+pLU1r+91TVR4CbZo0ePd7OBZ4w5qk/A1xcVTdV1VeBi4GTVqO+NfxantieEvTPa2+nzpnj7d64n1yYRmj8KsPZ3TgFfDDJp9rPP6y0hWyT77VpB/gtwA+sQm130bqMHgZcOmbyI5J8Jsn7kzxkVQsbzLff1sqxB8O7snfOMW3a2xHg0KraDsMfeuCQMW3WyvZcS6/lia2J36NP8nfAfcdM+m3g9cDvM2zc3wdezbAT7jKLMc9dts+N7q6+qnpva/PbwB3A2+eYzaOqaluSQ4CLk1zdzipWykK2yYput4VIsj/wV8Dzq+rWWZMvY+iGuL1dm/kb4KjVrI/599vUtyFA+2LizwEvHjN5LWzHhZr69lyDr+WJrYmgr6qfXki7JG8A3jdm0or+5MJ89bULRo8HTqjWiTdmHtva/c4k72HoWlnJg2Mh22SmzdYkewP34e5vt1dMkn0YQv7tVfXXs6ePBn9VXZTkz5Ksq6pV+yGsBey3tfJzH48FLquqHbMnrIXt2OxIclhVbW/dWzvHtNnKcE1hxv2AD69CbcCafS1PbM133czq7/x54Ioxzab2kwtJTgJeBPxcVX19jjb7JTlgZpjhos+49VhOC9kmFwAzn2p4IvChuQ7u5dauBbwRuKqqXjNHm/vOXDNIcizD8XrjatTXlrmQ/XYB8LT26ZvjgFtmuidW2VOYo9tm2ttxxOjxdhrw3jFtPgCcmOSg1k17Yhu34tbwa3ly074aPN8NeCvwWeByhgPlsDb+cOCikXYnM3xy41qGLpXVqm8zQ5/ip9vtz2fXx/DJl8+025WrVd+4bQK8nOFABtgX+Mu2Dp8AHrCK2+3RDG/JLx/ZdicDzwae3do8r22vzzBcHHvkKh97Y/fbrBrD8E92rm3H6cbVrLHVcG+G4L7PyLipbkeGPzrbge8wnKU/k+H6zyXANe3+4NZ2I/AXI8/91XZMbgaesYr1rdnX8qQ3fwJBkjq35rtuJEmTMeglqXMGvSR1zqCXpM4Z9JLUOYNekjpn0EtS5wx6aYwk/7H9kN6+7duQVyZ56LTrkpbCL0xJc0jyCoZvD98L2FpVfzDlkqQlMeilObTfCPok8E2Gnw24c8olSUti1400t4OB/Rn+A9a+U65FWjLP6KU5JLmA4T9z/RDDj+k9b8olSUuyJn6PXlprkjwNuKOq3pFkL+CfkvxUVX1o2rVJi+UZvSR1zj56SeqcQS9JnTPoJalzBr0kdc6gl6TOGfSS1DmDXpI69/8BdD/TX0JkKZoAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.hist(norm_sample)\n",
    "plt.xlabel('x')\n",
    "plt.title('distribution of Gaussian random sample')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Define likelihood function\n",
    "\n",
    "As you recall, the [likelihood function for a continuous distribution](https://en.wikipedia.org/wiki/Likelihood_function#Continuous_probability_distribution) such as this one:\n",
    "\n",
    "$$\\mathcal{L}(\\theta\\mid \\mathbf{x})=f_{\\theta}(\\mathbf{x})=P_{\\theta}(X=\\mathbf{x})$$\n",
    "or,\n",
    "$$\\mathcal{L}(\\theta\\mid \\mathbf{x})=\\prod_{i=1}^{N} f_{\\theta}(x_i) = \\prod_{i=1}^{N} P_{\\theta}(X=x_i).$$\n",
    "\n",
    "$f_{\\theta}(x_i)=P_{\\theta}(X=x)$ is essentially another way of writing the PDF for the continuous distribution given the parameterization $\\theta$. \n",
    "\n",
    "Here, because we are fitting the data to a [normal distribution](https://en.wikipedia.org/wiki/Normal_distribution), $\\theta=[\\mu,\\sigma]$. The PDF for the normal distribution is defined as:\n",
    "\n",
    "$$f_{\\mu,\\sigma}(x)=f(x\\mid\\mu,\\sigma)=\\frac{1}{\\sqrt{2\\pi\\sigma^2}}\\exp{\\left(-\\frac{(x-\\mu)^2}{2\\sigma^2}\\right)}.$$\n",
    "\n",
    "In practice, it is more common to optimize over the log-likelihood, because the product of many probability values can lead to floating point errors in computation:\n",
    "\n",
    "$$\\log (\\mathcal{L}(\\theta\\mid \\mathbf{x}))=\\sum_{i=1}^{N} \\log (f_{\\theta}(x_i)) = \\sum_{i=1}^{N} \\log (P_{\\theta})(X=x_i).$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Optimize \n",
    "\n",
    "Many built-in optimizers don't have a maximization function, only a minimization function. Therefore, instead of maximizing the log-likelihood, we will minimize the negative log-likelihood."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# define negative log-likelihood function\n",
    "def ll_norm(params, x):\n",
    "    '''params is a list of parameters to estimate: [mu, sigma]\n",
    "    x is list of normally distributed values described by estimated params'''\n",
    "    \n",
    "    mu = params[0]\n",
    "    sigma = params[1]\n",
    "    loglik = np.log((1/((2*math.pi*sigma**2)**0.5))*np.exp(-((x-mu)**2)/(2*sigma**2)))\n",
    "    return -loglik.sum()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Quick intro to essential arguments for [`scipy.optimize.minimize`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.minimize.html):\n",
    "- `fun`: function to minimize\n",
    "- `x0`: initial guesses for parameters to estimate, in the form of a list, tuple, or numpy array\n",
    "- `args`: any other variables to pass to the function, given in a list, tuple, or numpy array\n",
    "- `bounds`: bounds for parameters to be estimated, given as a list of lists or tuple of tuples, corresponds to params in `x0`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "      fun: 2509.8651698184385\n",
       " hess_inv: <2x2 LbfgsInvHessProduct with dtype=float64>\n",
       "      jac: array([ 0.00040927, -0.00027285])\n",
       "  message: b'CONVERGENCE: REL_REDUCTION_OF_F_<=_FACTR*EPSMCH'\n",
       "     nfev: 45\n",
       "      nit: 14\n",
       "   status: 0\n",
       "  success: True\n",
       "        x: array([4.02069038, 2.97703045])"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# minimize negative log-likelihood\n",
    "norm_res = optimize.minimize(fun = ll_norm,\n",
    "              x0 = [0,1],\n",
    "              args = (norm_sample),\n",
    "            bounds = ((None,None),(1e-6,None)))\n",
    "\n",
    "norm_res"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Interpreting the results \n",
    "- `fun`: minimum value of function at estimated parameters\n",
    "- `nfev`: number of function evaluations\n",
    "- `nit`: number of iterations\n",
    "- `success`: bool - did the optimizer run into an issue?\n",
    "- `x`: array of estimated parameters that minimize the function, corresponding to `x0`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "mu = 4, sigma = 3\n",
      "mu_est = 4.0207, sigma_est = 2.9770\n"
     ]
    }
   ],
   "source": [
    "print('mu = %d, sigma = %d' % (MU, SIGMA))\n",
    "print('mu_est = %.4f, sigma_est = %.4f' % (norm_res.x[0], norm_res.x[1]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "How do the estimated parameters $\\hat\\mu$ and $\\hat\\sigma$ compare to the actual parameter values used to generate the sample data? How will changing the size of the sample affect the precision of the estimation? "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# MLE for Poisson distributed data\n",
    "\n",
    "Now let's see how we can use MLE to estimate the parameters of the distribution from which a set of [Poisson distributed data](https://en.wikipedia.org/wiki/Poisson_distribution) came from. Poisson-distributed data is discrete and nonnegative. The Poisson distribution takes only one parameter, $\\lambda>0$, which represents the rate at which events occur. \n",
    "\n",
    "It is a discrete probability distribution that expresses the probability of a given number of events, $k$, occurring in a fixed interval of time or space, if the events occur with a known constant rate $\\lambda$ and independently of the time since the last event. \n",
    "\n",
    "Here, we will randomly generate a set of Poisson distributed data for a specified value of $\\lambda$.\n",
    "\n",
    "Note: The [scipy notation for the Poisson distribution](https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.poisson.html) uses $\\mu$ in place of $\\lambda$!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "LAMBDA = 0.5\n",
    "N = 1000\n",
    "poisson_sample = stats.poisson.rvs(mu = LAMBDA, size = N)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Visualize"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAEWCAYAAABollyxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAGU9JREFUeJzt3XmUXGWdxvHvIwnLCBKWBkMWGjQioAIhYhRnROI4LEpyFEaQgQSj0QF3z2hkHBVFjHPmiCKKJwoSkHVQTIQoMIEMokMkIKsBaTCQNpGEbGwGDf7mj/s2FJ3qrlvdVV1dr8/nnDp1l7fu/b11q5+69VZ1lSICMzPL10taXYCZmTWXg97MLHMOejOzzDnozcwy56A3M8ucg97MLHMO+iEi6UJJZ6bpv5f0QAO3/TNJ09P0DEm3NHDbJ0q6vlHbq2O/h0p6UNJTkqY1cLunS/p+o7Y3nEk6TFJ3q+sYKpIWS3p/q+sYjka0uoC/RRHxC2CfWu0kfRF4ZUT8S43tHdmIuiR1Ar8HRkbE5rTtS4BLGrH9On0JODcivlltpaTlwO7Ac8DTwELgIxHxVH8bjYizGlyn2bDnM/o2pkKux3BP4L4abd4ZEdsDE4HXA59relVNJGmrVtdgeco1JFpO0kGS7pD0pKQrgG0r1r3oJbWkz0j6Q2r7gKQpko4ATgfek4Yv7kptF0v6iqRfAs8Ae1d5ySpJ35K0UdL9kqZUrFgu6W0V81+U9MM0e3O63pD2+cbeQ0GS3iTptrTt2yS9qWLdYklflvTL1JfrJe3az330AUldktZJWiBpj7T8IWBv4Kepjm36u68j4g/Az4DXpNvvkba3Lm3/A9X6K2lbST+UtFbShtSf3dO6GZIeTv34vaQT0/KXSPqcpEckrZZ0kaQd07pOSSFpuqRHJT0u6d/76f+Fks6TtFDS08BbJR0t6TeSnpC0Ir2qo8z2JW2Xtrle0m8pnvwq97dvOkYbJN0n6ZhetXxHxTDgU+kYvlzSN9L27pd0UB/9kKSz0/2xUdLdknqORZn+nJLWrZf0IUmvT9vYIOncivYzUl1VH9tV6nqfpGVpu9dJ2rOvttmLCF8afAG2Bh4BPgGMBI4F/gKcmdYfBnSn6X2AFcAeab4TeEWa/iLww17bXgw8CuxPMfQ2Mi17f1o/A9hcse/3ABuBndP65cDbKrb3/D7SvgMYUbF+BnBLmt4ZWA+clPZ9QprfpaK2h4BXAdul+Tl93EeHA49TnI1vA3wLuLli/YvqrHL759cD4yjO/r+c5v8X+A7Fk+uBwBpgSpX+fhD4KfB3wFbAwcDLgJcCTwD7pHajgf3T9PuALoonou2BHwMX97r/vpf6fwDwLLBvH324MB2bQylOuraleGy8Ns2/DngMmFZm+8Ac4BfpOI0D7uWFx9nIVPfpFI/Pw4EnK/p4YToeB6c6bqQYxjs53TdnAjf10Y9/Am4HRgEC9gVGVzzWa/Xnu2mfbwc2AT8BdgPGAKuBt5R8bC/mhb+Daam/+1I8Vj8H/KrV2dCqi8/om2MyxQPxGxHxl4i4Critj7bPUQTdfpJGRsTyiHioxvYvjIj7ImJzRPylyvrVFfu+AngAOHqAfal0NPBgRFyc9n0ZcD/wzoo2P4iI30XEn4ArKYK2mhOBCyLijoh4Fvgs8EYV7xOU9RNJG4BbKML9LEnjgDcDn4mITRFxJ/B9iien3v4C7ELxPshzEXF7RDyR1v0VeI2k7SJiVUT0DCOdCHw9Ih6O4v2AzwLHS6p8v+uMiPhTRNwF3EURyH2ZHxG/jIi/pnoXR8Q9af5u4DLgLb1u09f2/xn4SkSsi4gVwDkVt5lM8cQ0JyL+HBE3AtdQPFn3uDrdB5uAq4FNEXFRRDwHXAFUPaNP9+MOwKsBRcSyiFgFULI/X059v57i/ZbLImJ1FK/UftFrv2Uf2x8Evppq2QycBRz4t3pW76Bvjj2AP0Q6tUgeqdYwIrqAj1Ocaa6WdHnPEEY/VtRYX23ftbZZxh5s2Y9HKM68evyxYvoZinCpua0Ummt7bauWaRExKiL2jIhT05PLHsC6iHiynxp7XAxcB1wuaaWk/0xPtk9TnC1+CFgl6VpJr65Wd5oeQfHGcI+y9wH0OpaS3iDpJklrJG1MNfQe/upr+3v02l5lnXsAKyLir73WV94vj1VM/6nKfNV+pCeNc4FvA49JmivpZXX0p579ln1s7wl8Mw3/bADWUbzaqOfxlQ0HfXOsAsZIUsWy8X01johLI+LNFA/OAL7Ws6qvm9TYf7V9r0zTT1MMVfR4eR3bXZlqrDQe+EON29XclqSXUpxdD2Rbvbe7s6QdKpZVrTGdFZ4REfsBbwLeQTFUQURcFxH/SDFscz/FcMkWdadtb+bF4VSP3vf5pcACYFxE7EgxrKEtblXdKoohm8raeqwExunFb94P9NhtISLOiYiDKYYUXwX8W1o1mP5U099ju9IK4IPpRKDnsl1E/GoQ+25bDvrm+D+KP/6PShoh6V3AIdUaStpH0uEq3nDcRHEG81xa/RjQqfo/WbNb2vdIScdRjFMuTOvupBhqGClpEsX7Bz3WUAxZ7N3HdhcCr5L03tSv9wD7UQwB1OtS4BRJB6a+nwUsiYjlA9jW89KQxa+Ar6p4s/V1wEyqfERU0lslvVbFp12eoBiCeE7S7pKOSU8+zwJP8cIxuQz4hKS9JG2f6r4iDQ80wg4Ur0g2SToEeG8dt70S+KyknSSNBT5SsW4JxZP8p9OxP4xiyO3ywRac3jx9g6SRaR+beOH+Gkx/qunvsV3puxT3xf6pxh1T+79JDvomiIg/A++iePNoPcUwwI/7aL4NxZtoj1O8JN+N4g0zgP9O12sl3VFHCUuACWmbXwGOjYi1ad1/AK9IdZ1BEbg9dT+T2v8yveSd3KtfaynOej9FMczyaeAdEfF4HbX1bGtRquVHFGeirwCOr3c7fTiB4o2+lRRjzV+IiBuqtHs5cBVFyC+jGOf/IcXfxafS7ddRjCmfmm5zAcWQz80Ub1Zu4sWBOlinAl+S9CTweYrwLusMiqGM3wPXpzqB5x+TxwBHUjwuvgOcHBH3N6Dml1G84lmf9r8W+K+0bjD9qaa/x/bzIuJqilfGl0t6guKN6Yb8v0k70ouHu8zMhidJMyg+VfPmVtfSbnxGb2aWOQe9mVnmPHRjZpY5n9GbmWVuWHx75a677hqdnZ2tLsPMrK3cfvvtj0dER612wyLoOzs7Wbp0aavLMDNrK5Kq/sd9bx66MTPLnIPezCxzDnozs8w56M3MMuegNzPLnIPezCxzpYJe0ihJV6XfaFym4rdEd5Z0g6QH0/VOqa0knaPitzrvljSxuV0wM7P+lD2j/ybw84h4NcXPli0DZgOLImICsCjNQ/FVoBPSZRZwXkMrNjOzutQM+vSTYP8AnA/F91pHxAZgKjAvNZtH8WO8pOUXReFWYJSk0Q2v3MzMSinzn7F7U/zy0A8kHUDxa+8fA3av+AHgVZJ2S+3H8OLfrexOy1ZVblTSLIozfsaP7/NX9mrqnH3tgG87WMvnNOL3ts3MmqvM0M0IYCJwXkQcRPFTYbP7aV/t9yC3+IrMiJgbEZMiYlJHR82vajAzswEqE/TdQHdELEnzV1EE/2M9QzLpenVF+8ofKB5L9R/vNTOzIVAz6CPij8AKSfukRVOA31L8svv0tGw6MD9NLwBOTp++mQxs7BniMTOzoVf22ys/AlwiaWvgYeAUiieJKyXNBB4Fen5hfSFwFNAFPJPamplZi5QK+oi4E5hUZdWUKm0DOG2QdZmZWYP4P2PNzDLnoDczy5yD3swscw56M7PMOejNzDLnoDczy5yD3swscw56M7PMOejNzDLnoDczy5yD3swscw56M7PMOejNzDLnoDczy5yD3swscw56M7PMOejNzDLnoDczy5yD3swscw56M7PMOejNzDLnoDczy5yD3swscw56M7PMOejNzDJXKuglLZd0j6Q7JS1Ny3aWdIOkB9P1Tmm5JJ0jqUvS3ZImNrMDZmbWv3rO6N8aEQdGxKQ0PxtYFBETgEVpHuBIYEK6zALOa1SxZmZWv8EM3UwF5qXpecC0iuUXReFWYJSk0YPYj5mZDULZoA/gekm3S5qVlu0eEasA0vVuafkYYEXFbbvTMjMza4ERJdsdGhErJe0G3CDp/n7aqsqy2KJR8YQxC2D8+PElyzAzs3qVOqOPiJXpejVwNXAI8FjPkEy6Xp2adwPjKm4+FlhZZZtzI2JSREzq6OgYeA/MzKxfNYNe0ksl7dAzDbwduBdYAExPzaYD89P0AuDk9OmbycDGniEeMzMbemWGbnYHrpbU0/7SiPi5pNuAKyXNBB4FjkvtFwJHAV3AM8ApDa/azMxKqxn0EfEwcECV5WuBKVWWB3BaQ6ozM7NB83/GmpllzkFvZpY5B72ZWeYc9GZmmXPQm5llzkFvZpY5B72ZWeYc9GZmmXPQm5llzkFvZpY5B72ZWeYc9GZmmXPQm5llzkFvZpY5B72ZWeYc9GZmmXPQm5llzkFvZpY5B72ZWeYc9GZmmXPQm5llzkFvZpY5B72ZWeYc9GZmmXPQm5llzkFvZpa50kEvaStJv5F0TZrfS9ISSQ9KukLS1mn5Nmm+K63vbE7pZmZWRj1n9B8DllXMfw04OyImAOuBmWn5TGB9RLwSODu1MzOzFikV9JLGAkcD30/zAg4HrkpN5gHT0vTUNE9aPyW1NzOzFih7Rv8N4NPAX9P8LsCGiNic5ruBMWl6DLACIK3fmNq/iKRZkpZKWrpmzZoBlm9mZrXUDHpJ7wBWR8TtlYurNI0S615YEDE3IiZFxKSOjo5SxZqZWf1GlGhzKHCMpKOAbYGXUZzhj5I0Ip21jwVWpvbdwDigW9IIYEdgXcMrNzOzUmqe0UfEZyNibER0AscDN0bEicBNwLGp2XRgfppekOZJ62+MiC3O6M3MbGgM5nP0nwE+KamLYgz+/LT8fGCXtPyTwOzBlWhmZoNRZujmeRGxGFicph8GDqnSZhNwXANqMzOzBvB/xpqZZc5Bb2aWOQe9mVnmHPRmZplz0JuZZc5Bb2aWOQe9mVnmHPRmZplz0JuZZc5Bb2aWOQe9mVnmHPRmZplz0JuZZc5Bb2aWOQe9mVnmHPRmZplz0JuZZc5Bb2aWOQe9mVnmHPRmZplz0JuZZc5Bb2aWOQe9mVnmHPRmZplz0JuZZa5m0EvaVtKvJd0l6T5JZ6Tle0laIulBSVdI2jot3ybNd6X1nc3tgpmZ9afMGf2zwOERcQBwIHCEpMnA14CzI2ICsB6YmdrPBNZHxCuBs1M7MzNrkZpBH4Wn0uzIdAngcOCqtHweMC1NT03zpPVTJKlhFZuZWV1GlGkkaSvgduCVwLeBh4ANEbE5NekGxqTpMcAKgIjYLGkjsAvweK9tzgJmAYwfP35wvbAh0zn72pbsd/mco1uyX7MclHozNiKei4gDgbHAIcC+1Zql62pn77HFgoi5ETEpIiZ1dHSUrdfMzOpU16duImIDsBiYDIyS1POKYCywMk13A+MA0vodgXWNKNbMzOpX5lM3HZJGpentgLcBy4CbgGNTs+nA/DS9IM2T1t8YEVuc0ZuZ2dAoM0Y/GpiXxulfAlwZEddI+i1wuaQzgd8A56f25wMXS+qiOJM/vgl1m5lZSTWDPiLuBg6qsvxhivH63ss3Acc1pDozMxs0/2esmVnmHPRmZplz0JuZZc5Bb2aWOQe9mVnmHPRmZplz0JuZZc5Bb2aWOQe9mVnmHPRmZplz0JuZZc5Bb2aWOQe9mVnmHPRmZplz0JuZZc5Bb2aWOQe9mVnmHPRmZplz0JuZZc5Bb2aWOQe9mVnmHPRmZplz0JuZZc5Bb2aWOQe9mVnmHPRmZpmrGfSSxkm6SdIySfdJ+lhavrOkGyQ9mK53Sssl6RxJXZLuljSx2Z0wM7O+lTmj3wx8KiL2BSYDp0naD5gNLIqICcCiNA9wJDAhXWYB5zW8ajMzK61m0EfEqoi4I00/CSwDxgBTgXmp2TxgWpqeClwUhVuBUZJGN7xyMzMrpa4xekmdwEHAEmD3iFgFxZMBsFtqNgZYUXGz7rSs97ZmSVoqaemaNWvqr9zMzEopHfSStgd+BHw8Ip7or2mVZbHFgoi5ETEpIiZ1dHSULcPMzOpUKugljaQI+Usi4sdp8WM9QzLpenVa3g2Mq7j5WGBlY8o1M7N6lfnUjYDzgWUR8fWKVQuA6Wl6OjC/YvnJ6dM3k4GNPUM8ZmY29EaUaHMocBJwj6Q707LTgTnAlZJmAo8Cx6V1C4GjgC7gGeCUhlZsZmZ1qRn0EXEL1cfdAaZUaR/AaYOsy8zMGsT/GWtmljkHvZlZ5hz0ZmaZc9CbmWXOQW9mljkHvZlZ5hz0ZmaZc9CbmWXOQW9mljkHvZlZ5hz0ZmaZc9CbmWXOQW9mljkHvZlZ5hz0ZmaZc9CbmWXOQW9mljkHvZlZ5hz0ZmaZc9CbmWXOQW9mljkHvZlZ5hz0ZmaZc9CbmWXOQW9mlrmaQS/pAkmrJd1bsWxnSTdIejBd75SWS9I5krok3S1pYjOLNzOz2sqc0V8IHNFr2WxgUURMABaleYAjgQnpMgs4rzFlmpnZQNUM+oi4GVjXa/FUYF6angdMq1h+URRuBUZJGt2oYs3MrH4DHaPfPSJWAaTr3dLyMcCKinbdadkWJM2StFTS0jVr1gywDDMzq6XRb8aqyrKo1jAi5kbEpIiY1NHR0eAyzMysx0CD/rGeIZl0vTot7wbGVbQbC6wceHlmZjZYAw36BcD0ND0dmF+x/OT06ZvJwMaeIR4zM2uNEbUaSLoMOAzYVVI38AVgDnClpJnAo8BxqflC4CigC3gGOKUJNZuZWR1qBn1EnNDHqilV2gZw2mCLMjOzxqkZ9GbWGp2zr23JfpfPObol+7Xm8VcgmJllzkFvZpY5B72ZWeYc9GZmmXPQm5llzkFvZpY5B72ZWeYc9GZmmXPQm5llzkFvZpY5B72ZWeYc9GZmmXPQm5llzkFvZpY5B72ZWeYc9GZmmXPQm5llzkFvZpY5B72ZWeYc9GZmmXPQm5llzkFvZpY5B72ZWeZGtLoAM7NW65x9bcv2vXzO0U3fh8/ozcwy15Sgl3SEpAckdUma3Yx9mJlZOQ0PeklbAd8GjgT2A06QtF+j92NmZuU044z+EKArIh6OiD8DlwNTm7AfMzMrQRHR2A1KxwJHRMT70/xJwBsi4sO92s0CZqXZfYAHBrjLXYHHB3jb4cZ9GX5y6Qe4L8PVYPqyZ0R01GrUjE/dqMqyLZ5NImIuMHfQO5OWRsSkwW5nOHBfhp9c+gHuy3A1FH1pxtBNNzCuYn4ssLIJ+zEzsxKaEfS3ARMk7SVpa+B4YEET9mNmZiU0fOgmIjZL+jBwHbAVcEFE3Nfo/VQY9PDPMOK+DD+59APcl+Gq6X1p+JuxZmY2vPg/Y83MMuegNzPLXNsEfa2vVZC0jaQr0volkjqHvspySvRlhqQ1ku5Ml/e3os5aJF0gabWke/tYL0nnpH7eLWniUNdYVom+HCZpY8Ux+fxQ11iGpHGSbpK0TNJ9kj5WpU1bHJeSfWmX47KtpF9Luiv15YwqbZqXYREx7C8Ub+o+BOwNbA3cBezXq82pwHfT9PHAFa2uexB9mQGc2+paS/TlH4CJwL19rD8K+BnF/1ZMBpa0uuZB9OUw4JpW11miH6OBiWl6B+B3VR5fbXFcSvalXY6LgO3T9EhgCTC5V5umZVi7nNGX+VqFqcC8NH0VMEVStX/earVsviIiIm4G1vXTZCpwURRuBUZJGj001dWnRF/aQkSsiog70vSTwDJgTK9mbXFcSvalLaT7+qk0OzJden8SpmkZ1i5BPwZYUTHfzZYH/Pk2EbEZ2AjsMiTV1adMXwDenV5WXyVpXJX17aBsX9vFG9NL759J2r/VxdSSXvofRHH2WKntjks/fYE2OS6StpJ0J7AauCEi+jwujc6wdgn6Ml+rUOqrF4aBMnX+FOiMiNcB/8MLz/Ltpl2OSRl3UHyvyAHAt4CftLiefknaHvgR8PGIeKL36io3GbbHpUZf2ua4RMRzEXEgxbcFHCLpNb2aNO24tEvQl/lahefbSBoB7MjwfClesy8RsTYink2z3wMOHqLaGi2br8OIiCd6XnpHxEJgpKRdW1xWVZJGUgTjJRHx4ypN2ua41OpLOx2XHhGxAVgMHNFrVdMyrF2CvszXKiwApqfpY4EbI72rMczU7Euv8dJjKMYm29EC4OT0KY/JwMaIWNXqogZC0st7xkslHULxt7O2tVVtKdV4PrAsIr7eR7O2OC5l+tJGx6VD0qg0vR3wNuD+Xs2almFt8Zux0cfXKkj6ErA0IhZQPCAultRF8Sx4fOsq7lvJvnxU0jHAZoq+zGhZwf2QdBnFpx52ldQNfIHiTSYi4rvAQopPeHQBzwCntKbS2kr05VjgXyVtBv4EHD9MTyQOBU4C7knjwQCnA+Oh7Y5Lmb60y3EZDcxT8cNMLwGujIhrhirD/BUIZmaZa5ehGzMzGyAHvZlZ5hz0ZmaZc9CbmWXOQW9mljkHvVkVkjr7+iZLs3bjoDczy5yD3qwGSXtL+o2k17e6FrOBcNCb9UPSPhTftXJKRNzW6nrMBqItvgLBrEU6gPnAuyPivlYXYzZQPqM369tGiu8HP7TVhZgNhs/ozfr2Z2AacJ2kpyLi0lYXZDYQDnqzfkTE05LeAdwg6emImN/qmszq5W+vNDPLnMfozcwy56A3M8ucg97MLHMOejOzzDnozcwy56A3M8ucg97MLHP/DxjrOtgaHwnZAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.hist(poisson_sample)\n",
    "plt.xlabel('k')\n",
    "plt.title('distribution of Poisson random sample')\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Define likelihood function\n",
    "\n",
    "For a [discrete distribution](https://en.wikipedia.org/wiki/Likelihood_function#Discrete_probability_distribution) such as the Poisson, the likelihood function doesn't change much, except it takes the products of the probability mass functions (PMF) rather than the probability distribution function (PDF):\n",
    "\n",
    "$$\\mathcal{L}(\\theta\\mid \\mathbf{x})=p_{\\theta}(\\mathbf{x})=P_{\\theta}(X=\\mathbf{x})$$\n",
    "or,\n",
    "$$\\mathcal{L}(\\theta\\mid \\mathbf{x})=\\prod_{i=1}^{N} p_{\\theta}(x_i) = \\prod_{i=1}^{N} P_{\\theta}(X=x_i).$$\n",
    "\n",
    "Again, we will need to define the [PMF for the Poisson distribution](https://en.wikipedia.org/wiki/Poisson_distribution#Probability_of_events_for_a_Poisson_distribution) to define our likelihood function. Because the Poisson is parameterized by $\\lambda$, $\\theta=[\\lambda]$ for our likelihood function.\n",
    "$$p_{\\lambda}(k)=p(k \\mid \\lambda)=\\exp^{-\\lambda}\\frac{\\lambda^k}{k!}$$\n",
    "\n",
    "Again, we will minimize over the negative log-likelihood function. The log-likelihood function for the Poisson is expressed as:\n",
    "$$\\log (\\mathcal{L}(\\theta\\mid \\mathbf{k}))=\\sum_{i=1}^{N} \\log (p_{\\theta}(k_i)) = \\sum_{i=1}^{N} \\log (P_{\\theta}(X=k_i))$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "def ll_pois(params, k):\n",
    "    '''params is list of parameters to estimate: [lambda]\n",
    "    k is list of Poisson distributed values described by estimated parameter'''\n",
    "    \n",
    "    lmbd = params[0]\n",
    "    loglik = np.log(np.exp(-lmbd)*(lmbd**k)/factorial(k))\n",
    "    return -loglik.sum()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note: When you want to set bounds but you only have one parameter to estimate, you need to format it as demonstrated below, or you will run into an error."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "      fun: 937.3870827020155\n",
       " hess_inv: <1x1 LbfgsInvHessProduct with dtype=float64>\n",
       "      jac: array([-0.0001819])\n",
       "  message: b'CONVERGENCE: REL_REDUCTION_OF_F_<=_FACTR*EPSMCH'\n",
       "     nfev: 56\n",
       "      nit: 17\n",
       "   status: 0\n",
       "  success: True\n",
       "        x: array([0.51099991])"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pois_res = optimize.minimize(fun = ll_pois,\n",
    "              x0 = [1e-6],\n",
    "              args = (poisson_sample),\n",
    "            bounds = ((1e-6,None),))\n",
    "\n",
    "pois_res"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "lambda = 0.5\n",
      "lambda_est = 0.5110\n"
     ]
    }
   ],
   "source": [
    "print('lambda = %.1f' % LAMBDA)\n",
    "print('lambda_est = %.4f' % pois_res.x[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "scipy has built-in [pdf](https://docs.scipy.org/doc/scipy-1.1.0/reference/generated/scipy.stats.rv_continuous.pdf.html) and [logpdf](https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.rv_continuous.logpdf.html) functions for a number of distributions, so in practice you could use those functions instead of implementing them yourself."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}