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": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "# define parameters\n",
    "MU = 4\n",
    "SIGMA = 3\n",
    "N = 1000\n",
    "\n",
    "# add some noise\n",
    "noise = np.random.randn(N)\n",
    "\n",
    "norm_sample = stats.norm.rvs(loc = MU, scale = SIGMA, size = N) + noise"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Visualize\n",
    "Let's take a look at what our sample data looks like."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAEWCAYAAACdaNcBAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAWxUlEQVR4nO3dfbRddX3n8fdHoKBAqzSBQkCCiE5BO7GToVY6HXymYgfsEg21Fjs42CW2dcaxBVc7ojXT1KkPHac6gxWhIg8ZHwZaaAVRSn0mUKpERKNECUmTACIPKkr8zh9733pycy733px7c7i/vF9r3XX22fu39/7uvc/9nN/5nXPPTVUhSWrLo8ZdgCRp7hnuktQgw12SGmS4S1KDDHdJapDhLkkNMtznQZLzk7yln/53SW6dw23/bZLT+ulXJPnUHG77ZUmumqvtzWK/xyX5WpL7k5y8q/e/Mwavw0KTpJI8cdx17ApJzkly4bjrGAfDfZ5V1T9U1ZOnazfTB2FV/UpVXTBqXUmW9r/kew5s+4NV9bxRt70T3gz8r6rar6r+37AGSVYk+XySB5Js6adfnSS7uFZg7q6DNF8M9wUinVav1+HA2qkWJnkd8OfA/wB+BjgI+G3gOOAndkWB4zL45CvNRqthsUsleVqSG5Pcl+RSYJ+BZccn2TBw/w+S3NG3vTXJs5OcALwBeGk/NPFPfdtrk6xM8mngu8AT+nmv3H73eVeS7yT5SpJnDyxYn+Q5A/cHXx1c19/e0+/zFycP8yR5RpLr+21fn+QZA8uuTfLHST7dH8tVSRY9zDn6T0nWJbk7yeVJDunnfx14AvDXfR17T1rvp+h69q+uqg9V1X3V+ceqellVPdi3OzHJPya5N8ntSc6Z6hpMPjdJjk2ypl93c5K39/P3SXJhkruS3NOfg4MGjv+V/fSRST7Rt7szyQeTPHbSvv5rki/25/LSJPswRH8NPp3kHUnuBs4ZdftJXp9kU5KNSf7j5POb5K+SbE3yzSR/ONGJmFTLPUm+0T8mXtGf4y15mKGpvt03+sfHbUleNovz9fr+eB5I8r4kB6UbCrsvyceTPK5vO/EK9Iz++Dal6wxMVdPTk3ymP55/SnL8VG0XvKryZ4Qfup7jN4H/DOwFvBj4IfCWfvnxwIZ++snA7cAh/f2lwJH99DnAhZO2fS3wLeAYYM9++9cCr+yXvwJ4aGDfLwW+AxzQL18PPGdge/+yj37fBew5sPwVwKf66QOAbwMv7/d9an//pwdq+zrwJODR/f1VU5yjZwF3Aj8P7A28C7huYPl2dU5a94T+GPcctnyg3fHAU+k6LD8HbAZOnnwNhu0T+Czw8n56P+Dp/fSrgL8GHgPsAfwb4CcHjn/iOjwReG5/bIvpnjjfOWlfXwAO6c/rLcBvT3EcE9f0d/rz/uhRtt+fv83AU4B9gYv66/7EfvlfAZcB+/ePia8Cp0+q5bf6438L3ePxL/pangfcB+w35Dj2Be4FntzfPxg4Zhbn63N0r9CWAFuAG4Gn9et8AnjjpMfxxf0+nwpsHbi25/Djx/wS4C7gBXSPk+f29xePO0fm48ee++ieThes76yqH1bVh4Drp2i7je7BeXSSvapqfVV9fZrtn19Va6vqoar64ZDlWwb2fSlwK3DiTh7LoBOBr1XVB/p9Xwx8BfjVgTbvr6qvVtX3gNXAsim29TLgvKq6sbqe9tnALyZZOoM6FgF3VtVDEzMGel7fS/LLAFV1bVV9qap+VFVfpPtl//czPNYfAk9Msqiq7q+qzw3M/2m6INxWVTdU1b2TV66qdVV1dVU9WFVbgbcP2ff/rKqNVXU33RPGVOcKYGNVvas/798bcfsvobtON1fVA3RhB0CSPeg6BGdX94poPfA2uif0CbdV1furahtwKXAY8Oa+lquAH9CF9TA/Ap6S5NFVtamq1s7ifL2rqjZX1R3APwCfr+7V2oPAR+mCftCbquqBqvoS8H66zshkvwFcWVVX9o+Tq4E1dGHfHMN9dIcAd1TfNeh9c1jDqloHvJbuF2xLkksmhicexu3TLB+27+m2OROHsONxfJOu9zPhnwemv0vX6512W1V1P12PackU7QfdBSzK9m/8PqOqHtsvmxhC+IUkn+yHF75DNyY/5TDRJKfTvQL5Sj/08sJ+/geAjwGX9C/535pkr8krJzmwv5Z3JLkXuHDIvmd6rmDSNR9x+4dM2t7gNV3Ej195Di4fvC6bB6a/B1BVk+ftcCz9E8lL6a7DpiRXJPlXszieyfuYbp+Tj3HY78DhwCl9x+CeJPcAv0T3qqI5hvvoNgFLku0+tfH4qRpX1UVV9Ut0D7QC/nRi0VSrTLP/Yfve2E8/QDekMOFnZrHdjX2Ngx4P3DHNetNuK8m+dD3imWzrs8CDwEnTtLsIuBw4rKp+CvjfwMR52e489D3WxRP3q+prVXUqcCDd9fhQkn37V0NvqqqjgWcALwR+c8i+/4TufP5cVf0kXQ9xlE/xTL42o2x/E11ve8LgY/NOulcnh09avjPXeAdV9bGqei5deH4FeG+/aK7PF+x4jBuHtLkd+EBVPXbgZ9+qWjXivh+RDPfRfZZuXPJ3k+yZ5NeAY4c1TPLkJM9K96bh9+l6INv6xZuBpZn9J2IO7Pe9V5JTgJ8FruyX3QSs6Jctp3s/YMJWupfNT5hiu1cCT0ry6/1xvRQ4GvibWdYHXfD+VpJl/bH/d7qX2eunW7Gq7gHeBLw7yYuT7JfkUUmW0Y2xTtgfuLuqvp/kWODXB5Z9Fdgn3ZuuewF/SDc8BkCS30iyuKp+BNzTz96W5JlJnto/GdxLF4Tb2NH+wP10b04vAV4/g3MyG6NsfzXwiiRHJ3kM8MaJBf1Qy2pgZZL9kxwO/Be6nvRI+jdA/0P/RP5gX//EuZuP8/VHSR6T5Bi69wguHdLmQuBXkzw/yR7p3jA/Psmhc7D/RxzDfURV9QPg1+jefPo23UvRj0zRfG9gFV2P6Z/pgvkN/bL/29/eleTGWZTweeCofpsrgRdX1V39sj8CjuzrehNdyE7U/d2+/af7l6hPn3Rcd9H1VF9HN/zx+8ALq+rOWdQ2sa1r+lo+TNeTPBJYMYv130oXOr9P9x7DZuD/AH8AfKZv9mrgzUnuA/4bXWhNrP+dfvlf0vVKHwAGPz1zArA2yf10H7lcUVXfp3ul8yG6YL8F+HuGB9+b6N4s/g5wBVNf/52109uvqr8F3kn3JuS6/nbQ79Cdj28An6J7jJw3esk8iu6xsxG4m25M/dX9svk4X39Pd3zXAH/Wvx+wnaq6ne4V4BvoOje30z2xNJmD2X64VpIWjv5N+duAvQbfdFejz1iStLsz3CWpQQ7LSFKD7LlLUoMeEV9KtGjRolq6dOm4y5CkBeWGG264s6oWD1v2iAj3pUuXsmbNmnGXIUkLSpKhfw0PDstIUpMMd0lqkOEuSQ0y3CWpQYa7JDXIcJekBhnuktQgw12SGmS4S1KDHhF/oSpNZ+lZV4xt3+tXzcX/G5d2LXvuktQgw12SGmS4S1KDDHdJapDhLkkNMtwlqUGGuyQ1yHCXpAYZ7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDTLcJalB04Z7ksOSfDLJLUnWJvm9fv45Se5IclP/84KBdc5Osi7JrUmeP58HIEna0Uy+z/0h4HVVdWOS/YEbklzdL3tHVf3ZYOMkRwMrgGOAQ4CPJ3lSVW2by8IlSVObtudeVZuq6sZ++j7gFmDJw6xyEnBJVT1YVbcB64Bj56JYSdLMzGrMPclS4GnA5/tZr0nyxSTnJXlcP28JcPvAaht4+CcDSdIcm3G4J9kP+DDw2qq6F3gPcCSwDNgEvG2i6ZDVa8j2zkiyJsmarVu3zrpwSdLUZhTuSfaiC/YPVtVHAKpqc1Vtq6ofAe/lx0MvG4DDBlY/FNg4eZtVdW5VLa+q5YsXLx7lGCRJk8zk0zIB3gfcUlVvH5h/8ECzFwE399OXAyuS7J3kCOAo4AtzV7IkaToz+bTMccDLgS8luamf9wbg1CTL6IZc1gOvAqiqtUlWA1+m+6TNmX5SRpJ2rWnDvao+xfBx9CsfZp2VwMoR6pIkjcC/UJWkBs1kWEb6F0vPumLcJUiaAXvuktQgw12SGmS4S1KDDHdJapDhLkkNMtwlqUGGuyQ1yHCXpAYZ7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDTLcJalBhrskNchwl6QG+c86pGmM6x+UrF914lj2qzbYc5ekBhnuktQgw12SGmS4S1KDDHdJapDhLkkNMtwlqUGGuyQ1yHCXpAYZ7pLUoGnDPclhST6Z5JYka5P8Xj//gCRXJ/laf/u4gXXOTrIuya1Jnj+fByBJ2tFMeu4PAa+rqp8Fng6cmeRo4Czgmqo6Crimv0+/bAVwDHAC8O4ke8xH8ZKk4aYN96raVFU39tP3AbcAS4CTgAv6ZhcAJ/fTJwGXVNWDVXUbsA44dq4LlyRNbVZj7kmWAk8DPg8cVFWboHsCAA7smy0Bbh9YbUM/b/K2zkiyJsmarVu3zr5ySdKUZhzuSfYDPgy8tqrufbimQ+bVDjOqzq2q5VW1fPHixTMtQ5I0AzMK9yR70QX7B6vqI/3szUkO7pcfDGzp528ADhtY/VBg49yUK0maiZl8WibA+4BbqurtA4suB07rp08DLhuYvyLJ3kmOAI4CvjB3JUuSpjOT/8R0HPBy4EtJburnvQFYBaxOcjrwLeAUgKpam2Q18GW6T9qcWVXb5rxySdKUpg33qvoUw8fRAZ49xTorgZUj1CVJGoF/oSpJDTLcJalBhrskNchwl6QGGe6S1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQYa7JDXIcJekBhnuktQgw12SGmS4S1KDDHdJapDhLkkNMtwlqUGGuyQ1yHCXpAYZ7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDTLcJalBhrskNchwl6QGGe6S1KBpwz3JeUm2JLl5YN45Se5IclP/84KBZWcnWZfk1iTPn6/CJUlTm0nP/XzghCHz31FVy/qfKwGSHA2sAI7p13l3kj3mqlhJ0sxMG+5VdR1w9wy3dxJwSVU9WFW3AeuAY0eoT5K0E0YZc39Nki/2wzaP6+ctAW4faLOhn7eDJGckWZNkzdatW0coQ5I02c6G+3uAI4FlwCbgbf38DGlbwzZQVedW1fKqWr548eKdLEOSNMxOhXtVba6qbVX1I+C9/HjoZQNw2EDTQ4GNo5UoSZqtnQr3JAcP3H0RMPFJmsuBFUn2TnIEcBTwhdFKlCTN1p7TNUhyMXA8sCjJBuCNwPFJltENuawHXgVQVWuTrAa+DDwEnFlV2+andEnSVKYN96o6dcjs9z1M+5XAylGKkiSNxr9QlaQGGe6S1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQYa7JDXIcJekBhnuktQgw12SGmS4S1KDDHdJapDhLkkNMtwlqUGGuyQ1yHCXpAYZ7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDTLcJalBhrskNchwl6QGGe6S1CDDXZIaZLhLUoP2nK5BkvOAFwJbquop/bwDgEuBpcB64CVV9e1+2dnA6cA24Her6mPzUrnUuKVnXTG2fa9fdeLY9q25MZOe+/nACZPmnQVcU1VHAdf090lyNLACOKZf591J9pizaiVJMzJtuFfVdcDdk2afBFzQT18AnDww/5KqerCqbgPWAcfOUa2SpBmadlhmCgdV1SaAqtqU5MB+/hLgcwPtNvTzdpDkDOAMgMc//vE7WcbuaZwv1yUtDHP9hmqGzKthDavq3KpaXlXLFy9ePMdlSNLubWfDfXOSgwH62y39/A3AYQPtDgU27nx5kqSdsbPhfjlwWj99GnDZwPwVSfZOcgRwFPCF0UqUJM3WTD4KeTFwPLAoyQbgjcAqYHWS04FvAacAVNXaJKuBLwMPAWdW1bZ5ql2SNIVpw72qTp1i0bOnaL8SWDlKUZKk0fgXqpLUIMNdkhpkuEtSgwx3SWqQ4S5JDTLcJalBhrskNchwl6QGGe6S1CDDXZIaZLhLUoMMd0lqkOEuSQ0y3CWpQYa7JDXIcJekBhnuktQgw12SGmS4S1KDDHdJapDhLkkNMtwlqUGGuyQ1yHCXpAYZ7pLUIMNdkhpkuEtSgwx3SWqQ4S5JDdpzlJWTrAfuA7YBD1XV8iQHAJcCS4H1wEuq6tujlSlJmo256Lk/s6qWVdXy/v5ZwDVVdRRwTX9fkrQLzcewzEnABf30BcDJ87APSdLDGDXcC7gqyQ1JzujnHVRVmwD62wOHrZjkjCRrkqzZunXriGVIkgaNNOYOHFdVG5McCFyd5CszXbGqzgXOBVi+fHmNWIckacBIPfeq2tjfbgE+ChwLbE5yMEB/u2XUIiVJs7PT4Z5k3yT7T0wDzwNuBi4HTuubnQZcNmqRkqTZGWVY5iDgo0kmtnNRVf1dkuuB1UlOB74FnDJ6mZKk2djpcK+qbwD/esj8u4Bnj1KUJGk0/oWqJDXIcJekBhnuktQgw12SGmS4S1KDDHdJapDhLkkNMtwlqUGGuyQ1yHCXpAYZ7pLUoFG/z323tvSsK8ZdgiQNZc9dkhpkuEtSgxyWkbSDcQ05rl914lj22yJ77pLUIMNdkhpkuEtSgwx3SWqQ4S5JDTLcJalBhrskNchwl6QGGe6S1CDDXZIaZLhLUoMMd0lqkOEuSQ3yWyElPWL4bZRzZ97CPckJwJ8DewB/WVWr5mtf/kckSdrevAzLJNkD+AvgV4CjgVOTHD0f+5Ik7Wi+eu7HAuuq6hsASS4BTgK+PE/7k6SdNs5X//M1JDRf4b4EuH3g/gbgFwYbJDkDOKO/e3+SW+eplplaBNw55hrGzXPgOZjgedhF5yB/OtLqh0+1YL7CPUPm1XZ3qs4Fzp2n/c9akjVVtXzcdYyT58BzMMHzsPDPwXx9FHIDcNjA/UOBjfO0L0nSJPMV7tcDRyU5IslPACuAy+dpX5KkSeZlWKaqHkryGuBjdB+FPK+q1s7HvubQI2aIaIw8B56DCZ6HBX4OUlXTt5IkLSh+/YAkNchwl6QGGe69JOckuSPJTf3PC8Zd066U5IQktyZZl+SscdczDknWJ/lSf/3XjLueXSXJeUm2JLl5YN4BSa5O8rX+9nHjrHG+TXEOFnQmGO7be0dVLet/rhx3MbuKXxexnWf213/Bfr55J5wPnDBp3lnANVV1FHBNf79l57PjOYAFnAmGu2Dg6yKq6gfAxNdFaDdQVdcBd0+afRJwQT99AXDyLi1qF5viHCxohvv2XpPki/1LtKZfhk4y7OsiloyplnEq4KokN/Rfj7E7O6iqNgH0tweOuZ5xWbCZsFuFe5KPJ7l5yM9JwHuAI4FlwCbgbWMtdtea9usidhPHVdXP0w1PnZnkl8ddkMZqQWfCbvXPOqrqOTNpl+S9wN/MczmPJH5dBFBVG/vbLUk+Sjdcdd14qxqbzUkOrqpNSQ4Gtoy7oF2tqjZPTC/ETNiteu4Pp38AT3gRcPNUbRu0239dRJJ9k+w/MQ08j93rMTDZ5cBp/fRpwGVjrGUsFnom7FY992m8NckyuuGI9cCrxlvOrrNAvy5irh0EfDQJdL8XF1XV3423pF0jycXA8cCiJBuANwKrgNVJTge+BZwyvgrn3xTn4PiFnAl+/YAkNchhGUlqkOEuSQ0y3CWpQYa7JDXIcJekBhnuktQgw12SGmS4S0Mk+bf9F0bt0//16tokTxl3XdJM+UdM0hSSvAXYB3g0sKGq/mTMJUkzZrhLU+i/Z+d64PvAM6pq25hLkmbMYRlpagcA+wH70/XgpQXDnrs0hSSX0/1XqiOAg6vqNWMuSZoxvxVSGiLJbwIPVdVF/f+Y/UySZ1XVJ8ZdmzQT9twlqUGOuUtSgwx3SWqQ4S5JDTLcJalBhrskNchwl6QGGe6S1KD/DwSVPUvWxXVBAAAAAElFTkSuQmCC\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) is:\n",
    "\n",
    "$$\\mathcal{L}(\\theta\\mid \\mathbf{x})=f_{\\theta}(\\mathbf{x})=P_{\\theta}(X=\\mathbf{x})$$\n",
    "which can also be expressed as,\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": 13,
   "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": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "      fun: 2570.9590091352106\n",
       " hess_inv: <2x2 LbfgsInvHessProduct with dtype=float64>\n",
       "      jac: array([0.00000000e+00, 4.54747351e-05])\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.13222182, 3.16458052])"
      ]
     },
     "execution_count": 14,
     "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": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "mu = 4, sigma = 3\n",
      "mu_est = 4.1322, sigma_est = 3.1646\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? \n",
    "\n",
    "##### MLE with smaller sample\n",
    "Let's try taking a subset of 100 points from our sample data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "      fun: 255.70201754960678\n",
       " hess_inv: <2x2 LbfgsInvHessProduct with dtype=float64>\n",
       "      jac: array([1.42108547e-05, 5.11590770e-05])\n",
       "  message: b'CONVERGENCE: REL_REDUCTION_OF_F_<=_FACTR*EPSMCH'\n",
       "     nfev: 48\n",
       "      nit: 15\n",
       "   status: 0\n",
       "  success: True\n",
       "        x: array([4.75188616, 3.12077813])"
      ]
     },
     "execution_count": 16,
     "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[:100]),\n",
    "            bounds = ((None,None),(1e-6,None)))\n",
    "\n",
    "norm_res"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What about a subset of 10 points?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "      fun: 22.34460957137232\n",
       " hess_inv: <2x2 LbfgsInvHessProduct with dtype=float64>\n",
       "      jac: array([-7.10542736e-07,  3.55271368e-07])\n",
       "  message: b'CONVERGENCE: NORM_OF_PROJECTED_GRADIENT_<=_PGTOL'\n",
       "     nfev: 54\n",
       "      nit: 16\n",
       "   status: 0\n",
       "  success: True\n",
       "        x: array([5.32609591, 2.26035613])"
      ]
     },
     "execution_count": 17,
     "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[:10]),\n",
    "            bounds = ((None,None),(1e-6,None)))\n",
    "\n",
    "norm_res"
   ]
  },
  {
   "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": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "LAMBDA = 0.5\n",
    "N = 1000\n",
    "\n",
    "poisson_sample = stats.poisson.rvs(mu = LAMBDA, size = N)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Visualize"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAEWCAYAAACdaNcBAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAYlklEQVR4nO3dfbRddX3n8ffHJALDg0C5xJAEAk5Egq2gEak4rRIrsVDI6sgYqjZYbOySOthxjU0cZyrVWDprltXW0i6KlCgPMaNSUrFCVoRB1AIBQQ0hJUIkMSkJzw9KNOEzf+xf8HBzbu65uefkcH/5vNbKOvvht3/7u889+Zx9fudhyzYREVGXl/S7gIiI6L6Ee0REhRLuEREVSrhHRFQo4R4RUaGEe0REhRLue4ikyyV9skz/J0lrutj3v0iaV6bPlXRLF/t+l6QbutXfCPZ7iqT7JD0taU4X+/2opEu71d+LmaQ3S9rQ7zr2FEk3SXpfv+t4sRjf7wL2Rra/BRw7XDtJHwf+o+13D9Pf27tRl6RpwAPABNvbSt9XAld2o/8R+nPgc7Y/226lpHXARGA78AzwdeCDtp/eVae2P9XlOiNelHLmPoapUevf8Chg1TBtfsf2AcBrgdcDH+t5VT0kaVy/a4h61BoMfSfpREl3SnpK0peAfVvWveDlsqQ/lfST0naNpFmSZgMfBd5ZhibuLm1vkrRI0reBnwLHtHk5Kkl/I+kJSfdKmtWyYp2kt7bMf1zSFWX25nL7eNnnrw8e5pH0Rkm3l75vl/TGlnU3SfqEpG+XY7lB0mG7uI/+UNJaSY9KWibpiLL8R8AxwD+XOvbZ1X1t+yfAvwCvLtsfUfp7tPT/h+2OV9K+kq6Q9Iikx8vxTCzrzpV0fzmOByS9qyx/iaSPSfqxpM2SviDpZWXdNEmWNE/Sg5IelvQ/dnH8l0v6O0lfl/QM8BZJp0v6nqQnJa0vr97opH9J+5U+H5N0D80TXuv+jit/o8clrZJ05qBaLlYzxPd0+Ru+XNJnSn/3SjpxiOOQpL8q98cTkr4vacffopPjeW9Z95ikP5L0+tLH45I+19L+3FJX28d2m7r+QNLq0u/1ko4aqm2VbOdfl/8BLwV+DPwJMAF4B/AL4JNl/ZuBDWX6WGA9cESZnwa8okx/HLhiUN83AQ8Cx9MMq00oy95X1p8LbGvZ9zuBJ4BDy/p1wFtb+nt+H2XfBsa3rD8XuKVMHwo8Bryn7PucMv8rLbX9CHglsF+Zv2iI++hU4GGas+59gL8Bbm5Z/4I622z//HpgKs1Z/ifK/P8DLqZ5Qj0B2ALManO87wf+GfgPwDjgdcBBwP7Ak8Cxpd0k4Pgy/QfAWponnwOArwJfHHT//UM5/tcAW4HjhjiGy8vf5hSaE619aR4bv1rmfw14CJjTSf/ARcC3yt9pKvBDfvk4m1Dq/ijN4/NU4KmWY7y8/D1eV+r4Js0Q3e+X++aTwI1DHMdpwB3AwYCA44BJLY/14Y7n78s+3wY8C/wTcDgwGdgM/GaHj+2b+OX/gznleI+jeax+DPhOv7NhT/7LmXtvnEzz4PuM7V/Y/jJw+xBtt9OE2wxJE2yvs/2jYfq/3PYq29ts/6LN+s0t+/4SsAY4fTePpdXpwH22v1j2fTVwL/A7LW3+0fa/2f4ZsJQmXNt5F3CZ7TttbwUWAr+uZty/U/8k6XHgFppA/5SkqcCbgD+1/aztu4BLaZ6QBvsF8Cs072tst32H7SfLuueAV0vaz/Ym2zuGiN4FfNr2/W7G9xcCcyW1vn91oe2f2b4buJsmhIdyre1v236u1HuT7R+U+e8DVwO/OWibofr/L8Ai24/aXg/8dcs2J9M8GV1k++e2vwl8jeYJeodryn3wLHAN8KztL9jeDnwJaHvmXu7HA4FXAbK92vYmgA6P5xPl2G+gef/katub3bwi+9ag/Xb62H4/8Bellm3Ap4AT9qaz94R7bxwB/MTlFKL4cbuGttcCH6I5o9wsacmO4YldWD/M+nb7Hq7PThzBzsfxY5ozrB3+vWX6pzSBMmxfJSgfGdTXcObYPtj2UbY/UJ5QjgAetf3ULmrc4YvA9cASSRsl/e/yBPsMzVnhHwGbJF0n6VXt6i7T42ne3N2h0/sABv0tJb1B0o2Stkh6otQweGhrqP6PGNRfa51HAOttPzdofev98lDL9M/azLc9jvJE8Tngb4GHJF0i6aARHM9I9tvpY/so4LNlaOdx4FGaVxUjeXyNaQn33tgETJaklmVHDtXY9lW230TzgDTwlztWDbXJMPtvt++NZfoZmmGIHV4+gn43lhpbHQn8ZJjthu1L0v40Z9G709fgfg+VdGDLsrY1lrO/C23PAN4InEEzDIHt623/Fs2QzL00QyE71V363sYLA2kkBt/nVwHLgKm2X0YzZKGdtmpvE81wTGttO2wEpuqFb8Dv7t9uJ7b/2vbraIYLXwn897JqNMfTzq4e263WA+8vT/47/u1n+zuj2PeYknDvje/S/If/r5LGS/pd4KR2DSUdK+lUNW8aPktzprK9rH4ImKaRfyLm8LLvCZLOphl3/HpZdxfNMMIESTNp3g/YYQvNcMQxQ/T7deCVkn6vHNc7gRk0L+9H6irgvZJOKMf+KeBW2+t2o6/nleGI7wB/oeYN018DzqPNxzklvUXSr6r5lMqTNMML2yVNlHRmecLZCjzNL/8mVwN/IuloSQeUur9UXvp3w4E0rzyelXQS8Hsj2HYpsFDSIZKmAB9sWXcrzRP7R8rf/s00w2lLRltweQP0DZImlH08yy/vr9EcTzu7emy3+nua++L4UuPLSvu9RsK9B2z/HPhdmjeAHqN5if/VIZrvQ/NG2MM0L7cPp3nTC+D/lttHJN05ghJuBaaXPhcB77D9SFn3P4FXlLoupAnZHXX/tLT/dnk5e/Kg43qE5uz2wzRDKB8BzrD98Ahq29HXilLLV2jOOF8BzB1pP0M4h+bNuo00Y8d/Znt5m3YvB75ME+yracbtr6D5f/Hhsv2jNGPEHyjbXEYznHMzzRuOz/LCEB2tDwB/Lukp4H/RBHanLqQZpngAuKHUCTz/mDwTeDvN4+Ji4Pdt39uFmg+ieWXzWNn/I8D/KetGczzt7Oqx/Tzb19C8Al4i6UmaN5e78n2QsUIvHL6KiHhxknQuzadh3tTvWsaCnLlHRFQo4R4RUaEMy0REVChn7hERFXpR/CrkYYcd5mnTpvW7jIiIMeWOO+542PZAu3UvinCfNm0aK1eu7HcZERFjiqS233yHDMtERFQp4R4RUaGEe0REhRLuEREVSrhHRFQo4R4RUaGOwl3SwZK+XK5ZuFrNtTUPlbRc0n3l9pCW9gvVXLtyjaTTeld+RES00+mZ+2eBb9h+Fc0lvVYDC4AVtqcDK8o8kmbQ/HTr8cBs4GLlqu4REXvUsOFeLpf1G8DnofldaNuPA2cBi0uzxTQXpKUsX2J7q+0HaC5S2/ZCFRER0RudfEP1GJor9PyjpNfQXOX8AmBiy0VwN0k6vLSfDPxry/YbaHPdQknzgfkARx455BXoOjJtwXWj2n53rbuoG9ecjojovk6GZcYDrwX+zvaJNJfRWrCL9u2uj7jTT0/avsT2TNszBwba/jRCRETspk7CfQOwwfatZf7LNGH/kKRJAOV2c0v71ov0TqH9BWwjIqJHhg132/8OrJd0bFk0C7iH5orm88qyecC1ZXoZzQWY95F0NM31Dm/ratUREbFLnf4q5AeBKyW9FLgfeC/NE8NSSecBDwJnA9heJWkpzRPANuB829vbdxsREb3QUbjbvguY2WbVrCHaL6K5MnlERPRBvqEaEVGhhHtERIUS7hERFUq4R0RUKOEeEVGhhHtERIUS7hERFUq4R0RUKOEeEVGhhHtERIUS7hERFUq4R0RUKOEeEVGhhHtERIUS7hERFUq4R0RUKOEeEVGhhHtERIUS7hERFUq4R0RUKOEeEVGhhHtERIUS7hERFUq4R0RUKOEeEVGhjsJd0jpJP5B0l6SVZdmhkpZLuq/cHtLSfqGktZLWSDqtV8VHRER7Izlzf4vtE2zPLPMLgBW2pwMryjySZgBzgeOB2cDFksZ1seaIiBjGaIZlzgIWl+nFwJyW5Utsb7X9ALAWOGkU+4mIiBHqNNwN3CDpDknzy7KJtjcBlNvDy/LJwPqWbTeUZRERsYeM77DdKbY3SjocWC7p3l20VZtl3qlR8yQxH+DII4/ssIyIiOhER2futjeW283ANTTDLA9JmgRQbjeX5huAqS2bTwE2tunzEtszbc8cGBjY/SOIiIidDBvukvaXdOCOaeBtwA+BZcC80mwecG2ZXgbMlbSPpKOB6cBt3S48IiKG1smwzETgGkk72l9l+xuSbgeWSjoPeBA4G8D2KklLgXuAbcD5trf3pPqIiGhr2HC3fT/wmjbLHwFmDbHNImDRqKuLiIjdkm+oRkRUKOEeEVGhhHtERIUS7hERFUq4R0RUKOEeEVGhhHtERIUS7hERFUq4R0RUKOEeEVGhhHtERIUS7hERFUq4R0RUKOEeEVGhhHtERIUS7hERFUq4R0RUKOEeEVGhhHtERIUS7hERFUq4R0RUKOEeEVGhhHtERIUS7hERFUq4R0RUKOEeEVGhjsNd0jhJ35P0tTJ/qKTlku4rt4e0tF0oaa2kNZJO60XhERExtJGcuV8ArG6ZXwCssD0dWFHmkTQDmAscD8wGLpY0rjvlRkREJzoKd0lTgNOBS1sWnwUsLtOLgTkty5fY3mr7AWAtcFJ3yo2IiE50eub+GeAjwHMtyyba3gRQbg8vyycD61vabSjLXkDSfEkrJa3csmXLiAuPiIihDRvuks4ANtu+o8M+1WaZd1pgX2J7pu2ZAwMDHXYdERGdGN9Bm1OAMyX9NrAvcJCkK4CHJE2yvUnSJGBzab8BmNqy/RRgYzeLjoiIXRv2zN32QttTbE+jeaP0m7bfDSwD5pVm84Bry/QyYK6kfSQdDUwHbut65RERMaROztyHchGwVNJ5wIPA2QC2V0laCtwDbAPOt7191JVGRETHRhTutm8CbirTjwCzhmi3CFg0ytoiImI35RuqEREVSrhHRFQo4R4RUaGEe0REhRLuEREVSrhHRFQo4R4RUaGEe0REhRLuEREVSrhHRFQo4R4RUaGEe0REhRLuEREVSrhHRFQo4R4RUaGEe0REhRLuEREVSrhHRFQo4R4RUaGEe0REhRLuEREVSrhHRFRofL8LiJGbtuC6vu173UWn923fEdG5nLlHRFQo4R4RUaFhw13SvpJuk3S3pFWSLizLD5W0XNJ95faQlm0WSloraY2k03p5ABERsbNOzty3Aqfafg1wAjBb0snAAmCF7enAijKPpBnAXOB4YDZwsaRxvSg+IiLaGzbc3Xi6zE4o/wycBSwuyxcDc8r0WcAS21ttPwCsBU7qatUREbFLHY25Sxon6S5gM7Dc9q3ARNubAMrt4aX5ZGB9y+YbyrLBfc6XtFLSyi1btozmGCIiYpCOwt32dtsnAFOAkyS9ehfN1a6LNn1eYnum7ZkDAwOdVRsRER0Z0adlbD8O3EQzlv6QpEkA5XZzabYBmNqy2RRg46grjYiIjnXyaZkBSQeX6f2AtwL3AsuAeaXZPODaMr0MmCtpH0lHA9OB27pdeEREDK2Tb6hOAhaXT7y8BFhq+2uSvgsslXQe8CBwNoDtVZKWAvcA24DzbW/vTfkREdHOsOFu+/vAiW2WPwLMGmKbRcCiUVcXERG7Jd9QjYioUMI9IqJCCfeIiAol3CMiKpRwj4ioUMI9IqJCCfeIiAol3CMiKpRwj4ioUMI9IqJCCfeIiAol3CMiKpRwj4ioUMI9IqJCCfeIiAol3CMiKpRwj4ioUMI9IqJCCfeIiAol3CMiKpRwj4ioUMI9IqJCCfeIiAol3CMiKpRwj4ioUMI9IqJCw4a7pKmSbpS0WtIqSReU5YdKWi7pvnJ7SMs2CyWtlbRG0mm9PICIiNhZJ2fu24AP2z4OOBk4X9IMYAGwwvZ0YEWZp6ybCxwPzAYuljSuF8VHRER7w4a77U227yzTTwGrgcnAWcDi0mwxMKdMnwUssb3V9gPAWuCkbhceERFDG9GYu6RpwInArcBE25ugeQIADi/NJgPrWzbbUJYN7mu+pJWSVm7ZsmXklUdExJA6DndJBwBfAT5k+8ldNW2zzDstsC+xPdP2zIGBgU7LiIiIDnQU7pIm0AT7lba/WhY/JGlSWT8J2FyWbwCmtmw+BdjYnXIjIqITnXxaRsDngdW2P92yahkwr0zPA65tWT5X0j6SjgamA7d1r+SIiBjO+A7anAK8B/iBpLvKso8CFwFLJZ0HPAicDWB7laSlwD00n7Q53/b2rlceERFDGjbcbd9C+3F0gFlDbLMIWDSKuiIiYhTyDdWIiAol3CMiKpRwj4ioUMI9IqJCCfeIiAol3CMiKpRwj4ioUMI9IqJCCfeIiAol3CMiKpRwj4ioUMI9IqJCCfeIiAol3CMiKpRwj4ioUMI9IqJCCfeIiAol3CMiKpRwj4ioUMI9IqJCCfeIiAol3CMiKpRwj4ioUMI9IqJCCfeIiAoNG+6SLpO0WdIPW5YdKmm5pPvK7SEt6xZKWitpjaTTelV4REQMrZMz98uB2YOWLQBW2J4OrCjzSJoBzAWOL9tcLGlc16qNiIiODBvutm8GHh20+CxgcZleDMxpWb7E9lbbDwBrgZO6VGtERHRod8fcJ9reBFBuDy/LJwPrW9ptKMt2Imm+pJWSVm7ZsmU3y4iIiHa6/Yaq2ixzu4a2L7E90/bMgYGBLpcREbF3291wf0jSJIByu7ks3wBMbWk3Bdi4++VFRMTu2N1wXwbMK9PzgGtbls+VtI+ko4HpwG2jKzEiIkZq/HANJF0NvBk4TNIG4M+Ai4Clks4DHgTOBrC9StJS4B5gG3C+7e09qj0iIoYwbLjbPmeIVbOGaL8IWDSaoiIiYnTyDdWIiAol3CMiKpRwj4ioUMI9IqJCCfeIiAoN+2mZiBeDaQuu69u+1110et/2HbG7cuYeEVGhhHtERIUS7hERFUq4R0RUKOEeEVGhhHtERIUS7hERFUq4R0RUKOEeEVGhhHtERIUS7hERFUq4R0RUKOEeEVGhhHtERIUS7hERFUq4R0RUKOEeEVGhhHtERIUS7hERFUq4R0RUqGcXyJY0G/gsMA641PZFvdpXRI1yUfAYjZ6cuUsaB/wt8HZgBnCOpBm92FdEROysV2fuJwFrbd8PIGkJcBZwT4/2FxEV6OerlX7p1ask2e5+p9I7gNm231fm3wO8wfYft7SZD8wvs8cCa0axy8OAh0ex/Viztx0v5Jj3FjnmkTnK9kC7Fb06c1ebZS94FrF9CXBJV3YmrbQ9sxt9jQV72/FCjnlvkWPunl59WmYDMLVlfgqwsUf7ioiIQXoV7rcD0yUdLemlwFxgWY/2FRERg/RkWMb2Nkl/DFxP81HIy2yv6sW+iq4M74whe9vxQo55b5Fj7pKevKEaERH9lW+oRkRUKOEeEVGhMR3ukmZLWiNpraQF/a6n1yRdJmmzpB/2u5Y9RdJUSTdKWi1plaQL+l1Tr0naV9Jtku4ux3xhv2vaEySNk/Q9SV/rdy17iqR1kn4g6S5JK7va91gdcy8/cfBvwG/RfPTyduAc29V+C1bSbwBPA1+w/ep+17MnSJoETLJ9p6QDgTuAOZX/nQXsb/tpSROAW4ALbP9rn0vrKUn/DZgJHGT7jH7XsydIWgfMtN31L26N5TP353/iwPbPgR0/cVAt2zcDj/a7jj3J9ibbd5bpp4DVwOT+VtVbbjxdZieUf2PzLKxDkqYApwOX9ruWWozlcJ8MrG+Z30Dl/+n3dpKmAScCt/a3kt4rQxR3AZuB5bZrP+bPAB8Bnut3IXuYgRsk3VF+kqVrxnK4D/sTB1EPSQcAXwE+ZPvJftfTa7a32z6B5tvdJ0mqdhhO0hnAZtt39LuWPjjF9mtpfkH3/DL02hVjOdzzEwd7iTLu/BXgSttf7Xc9e5Ltx4GbgNl9LqWXTgHOLOPPS4BTJV3R35L2DNsby+1m4Bqa4eauGMvhnp842AuUNxc/D6y2/el+17MnSBqQdHCZ3g94K3Bvf6vqHdsLbU+xPY3m//E3bb+7z2X1nKT9y4cEkLQ/8Daga5+EG7PhbnsbsOMnDlYDS3v8Ewd9J+lq4LvAsZI2SDqv3zXtAacA76E5m7ur/PvtfhfVY5OAGyV9n+YkZrntvebjgXuRicAtku4GbgOus/2NbnU+Zj8KGRERQxuzZ+4RETG0hHtERIUS7hERFUq4R0RUKOEeEVGhhHvEECRN25t+gTPqknCPiKhQwj2iA5KOKb81/vp+1xLRiYR7xDAkHUvz2zbvtX17v+uJ6MT4fhcQ8SI3AFwL/Ofaf94i6pIz94hde4LmugGn9LuQiJHImXvErv0cmANcL+lp21f1u6CITiTcI4Zh+5lyQYnlkp6xfW2/a4oYTn4VMiKiQhlzj4ioUMI9IqJCCfeIiAol3CMiKpRwj4ioUMI9IqJCCfeIiAr9fx3wwNb7UHzYAAAAAElFTkSuQmCC\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 statistical distributions, so in practice you could use those functions instead of implementing them yourself."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}