Maximum_likelihood.ipynb

{
 "metadata": {
  "name": "",
  "signature": "sha256:341b8fa35083e0a672401367dad8ff956fd442e65a7d8020c6fe3a4642aad208"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Maximum Likelihood Estimation"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Maximum likelihood estimation is one of the key techniques employed in statistical signal processing for a wide variety of applications from signal detection to parameter estimation. In the following, we consider a simple experiment and work through the details of maximum likelihood estimation to ensure that we understand the concept in one of its simplest applications."
     ]
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "Setting up the Coin Flipping Experiment"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Suppose we have coin and want to estimate the probability of heads ($p$) for it. The coin is Bernoulli distributed:\n",
      "\n",
      "$$ \\phi(x)= p^x (1-p)^{(1-x)} $$\n",
      "\n",
      "where $x$ is the outcome, *1* for heads and *0* for tails. The $n$ independent flips, we have the likelihood:\n",
      "\n",
      "$$ \\mathcal{L}(p|\\mathbf{x})= \\prod_{i=1}^n p^{ x_i }(1-p)^{1-x_i} $$\n",
      "\n",
      "This is basically notation. We have just substituted everything into $ \\phi(x)$ under the independent-trials assumption. \n",
      "\n",
      "The idea of *maximum likelihood* is to maximize this as the function of $p$ after plugging in all of the $x_i$ data. This means that our estimator, $\\hat{p}$ , is a function of the observed $x_i$ data, and as such, is a random variable with its own distribution."
     ]
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "Simulating the Experiment"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We need the following code to simulate coin flipping."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%matplotlib inline\n",
      "from __future__ import division\n",
      "from scipy.stats import bernoulli \n",
      "import numpy as np\n",
      "\n",
      "p_true=1/2 # this is the value we will try to estimate from the observed data\n",
      "fp=bernoulli(p_true)\n",
      "\n",
      "def sample(n=10):\n",
      "    'simulate coin flipping'\n",
      "    return fp.rvs(n)# flip it n times\n",
      "\n",
      "xs = sample(100) # generate some samples"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 1
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Now, we can write out the likelihood function using `sympy`"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import sympy\n",
      "from sympy.abc import x, z\n",
      "p=sympy.symbols('p',positive=True)\n",
      "\n",
      "L=p**x*(1-p)**(1-x)\n",
      "J=np.prod([L.subs(x,i) for i in xs]) # objective function to maximize"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 2
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Below, we find the maximum using basic calculus. Note that taking the `log` of $J$ makes the maximization problem tractable but doesn't change the extrema."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "logJ=sympy.expand_log(sympy.log(J))\n",
      "sol=sympy.solve(sympy.diff(logJ,p),p)[0]\n",
      "\n",
      "x=linspace(0,1,100)\n",
      "plot(x,map(sympy.lambdify(p,logJ,'numpy'),x),sol,logJ.subs(p,sol),'o',\n",
      "                                          p_true,logJ.subs(p,p_true),'s',)\n",
      "xlabel('$p$',fontsize=18)\n",
      "ylabel('Likelihood',fontsize=18)\n",
      "title('Estimate not equal to true value',fontsize=18)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stderr",
       "text": [
        "C:\\Python27\\lib\\site-packages\\numpy\\__init__.py:1: RuntimeWarning: divide by zero encountered in log\n",
        "  \"\"\"\n"
       ]
      },
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 3,
       "text": [
        "<matplotlib.text.Text at 0xcb707f0>"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAZgAAAElCAYAAADZb/T+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XtcTPn/B/DXTLqJJOU2leiyGXKLcltyGZdaVuzPyjUW\nLWvZ9V3W5bsrrKzL7ndZslhC1rL4Vrtrta65rlqxRHxFoUbZQomi2+f3x2cbjaaaqZk5U72fj8d5\n0DlnznnPmTrvOZ+riDHGQAghhGiZWOgACCGE1E6UYAghhOgEJRhCCCE6QQmGEEKITlCCIYQQohOU\nYAghhOgEJZgaLCAgAGIxfYSkYo6OjujXr5/QYRAA3t7eaN26tdBh6A3dnbQoOjoaYrG43MXY2Fjj\nY+7YsQPr1q1TuU0kEkEkElU37GoLCgpCZGSk0GForKJrW5uo+3uSnZ2NoKAgnDp1Sidx/PXXXwgK\nCsK9e/d0cvyawhD+ZvWlntAB1EZjx46Fj49PmfVVedrYsWMH7t27hzlz5pTZtnXrVmzevLlKMWrT\nsmXLEBAQgLffflvoUDRS0bWtTdTtS/3kyRMsW7YMYrEYffv21Xocf/31F5YtW4b+/fujVatWWj8+\nMTyUYHSgS5cuGDt2rNaOV943nnr16OOrrrr0bVJduh7cQ93jFxUVIT8/H+bm5jqNh+gOFZEJaNeu\nXfD09ETjxo3RoEEDODk5Yfz48cjMzATAy85Pnz6Nu3fvKhW1nT59GoDqOpiSdY8fP8aUKVNga2sL\nS0tLvP3220hLSwMAbN68GW3btoW5uTnatm2Ln3/+uUxsISEhGDRoEOzs7GBqaoqWLVtiwoQJSsUb\nJXEB/GmgdIylHTt2DIMGDULjxo1hbm6Ojh07qv3kVXKOpUuX4tdff0W3bt1gbm6Oli1bYv78+Sgq\nKirzmtOnT0Mmk8HKygr169eHh4cHtm/frrRPZde2Ipq8n61bt8LNzQ1mZmZwcXHBunXrEBoaWuZc\nFdWnicViTJ48WWmdOp+PJqKjo9GmTRsAwNKlSxXXo3R9QWFhIVatWgWpVApzc3PY2Nhg5MiRuHbt\nWqXHDwoKwpQpUwAA/fr1Uxy/ZF3J78/x48exfPlyODk5wdzcHPv37y/3GpR+3eufW3Z2Nj799FM4\nOzvDzMwMTZs2xdixY5GcnFxprJ9++inEYjHi4+PLbMvOzoa5uTn8/PwU6/bt24fhw4ejVatWMDMz\ng62tLfz8/FS+XpXy6shKitx37typtP7ly5cIDg5Gu3btYG5ujsaNG2P48OH466+/1DqfPtFXYB14\n/vy5IkmUZmpqioYNGwIAwsLCEBAQgD59+mD58uUwNzfH/fv3cfjwYWRkZMDGxgbr1q3DwoULkZmZ\niW+++UZxnLZt2yr+X9438CFDhsDe3h7Lly9HYmIi1q9fj+HDh8PPzw87duzA1KlTYWpqivXr1+Od\nd97BrVu34OjoqHj9V199hR49ekAmk8Ha2hrx8fH4/vvvceLECcTHx8Pa2hpNmzZFWFgYJkyYgD59\n+mD69Oll4tiyZQvef/999OzZE//+979hYWGBI0eOYMaMGbhz5w5Wr16t1jX97bffEBISghkzZmDq\n1KmIiIjA2rVr0bhxYyxcuFCx3y+//AI/Pz+0bNkSn3zyCRo2bIgff/wRU6dORVJSEr744gsAKPfa\nurm5VRiHJu/nm2++wdy5c9GpUyesXLkSz58/x9q1a2Fra6vy2BU9Tb2+TZ3PRxNSqRT/+c9/8PHH\nH2PkyJEYOXIkAKBBgwaKfcaNG4f9+/dj0KBB+OCDD5CWloaNGzeiR48eOHPmDDp16lTu8UeNGoX0\n9HRs2bIFixcvVvwOOzk5Ke33ySefoLCwEIGBgbC0tMQbb7xR7jUoT3Z2Nnr27ImUlBS89957aNeu\nHR48eICQkBB4eXnh4sWLcHBwKPf1AQEBWLNmDXbt2oU1a9Yobfvpp5/w8uVLBAQEKNZt3LgRNjY2\nCAwMRPPmzXH79m1s2bIFvXr1wqVLl+Ds7FxhvJXVkZXeVlBQgCFDhuCPP/7AxIkTMXv2bGRlZWHr\n1q3o1asXTp8+DQ8Pj0qukB4xojUnT55kIpGo3GXYsGGKff38/FijRo1YUVFRhcfs27cva926tcpt\nkyZNYiKRSOW6WbNmKa2fO3cuE4lEzMHBgeXk5CjWX716lYlEIrZw4UKl/XNzc8uc7/jx40wkErHV\nq1crrReJRGzy5Mll9n/w4AEzNTVl48aNK7Ntzpw5zMjIiCUlJal8byWSk5OZSCRiDRo0YPfu3VPa\n1r59e9aiRQvFz4WFhczBwYE1btyYpaWlKdbn5+ezXr16MSMjI5aYmKhYX9G1VUWT9/PkyRNWv359\n1q5dO5aXl6fYLzU1lTVo0ICJxWJ26tQpxXpVn2UJVddXk8+nVatWrF+/fpW+v5JrvXTp0jLbjhw5\nwkQiERszZozS+itXrrB69eqxN998s9Ljh4aGMpFIpPS+X9/m5uamdL1KlPc7puqYs2fPZvXr12dX\nr15V2vfevXvM0tKSBQQEVBprt27dWMuWLcv8ffbu3ZvZ2tqygoICxTpVn8WNGzeYqakpmzlzptJ6\nVb9z5X0+JfeTnTt3KtZ9/fXXTCQSsSNHjijt+/TpU+bg4MC8vb0rfW/6REVkOhAYGIhjx46VWVas\nWKHYx8rKCs+fP8evv/6qkzLvjz76SOnn3r17AwAmTZqk9K3U3d0dlpaWuH37ttL+JeXexcXFyM7O\nRmZmJjp06IBGjRohNjZWrRgOHDiA/Px8TJkyBZmZmUrLW2+9heLiYhw7dkytY40YMaLMt05vb2+k\np6cjNzcXABAXF4eUlBRMmTIFzZs3V+xnbGyM+fPno7i4uFqt3TR5P0eOHEFeXh4++OADmJmZKY4h\nkUgwbty4an/m2vh8NBEeHg4AWLx4sdL6Dh06YNiwYTh79iwePXpU7fPMmDFD6XppijGGH374AX36\n9EHLli2VPqP69evDy8sLR44cqfQ4kyZNQlpaGo4ePapYl5ycjPPnz8Pf31+p/rPks2CM4enTp8jM\nzISNjQ1cXV21/lns3r0bbdu2RZcuXZTe28uXLzFw4ECcPXsWL1++1Oo5q4OKyHTAxcUF/fv3r3Cf\nRYsW4fTp0xgxYgSaNGmCvn37YujQoXj33XeVEkBVlZSnl2jcuDEAqGyDb2VlVebmcOLECSxbtgyx\nsbF48eKF0rYnT56oFcONGzcAAAMHDlS5XSQS4e+//1brWK+/HwBo0qQJAODRo0eoX7++ony9Xbt2\nZfaVSqUAoFYZfHk0eT9JSUkAVBe5lS7irCptfD6aSE5OhpGRkcrYpVIpIiIikJycrPhMqsrV1bVa\nr8/IyMDjx4/x+++/l1sUaWRkVOlx/P398a9//Qu7du3C4MGDAfA6U8YYJk6cqLTv5cuX8dlnn+HU\nqVN4/vy50jZVv7fVcePGDbx48aLCYtbMzExIJBKtnreqKMEIxNnZGQkJCTh+/DiOHz+OU6dOYdq0\naViyZAlOnz5d7V/M8sp0y/vjKv2N+s8//8SgQYPg6uqKVatWoXXr1opvaWPGjEFxcbFaMZQcMyws\nDC1atFC5j7qdziq6KejiCbCi81T0fqr6uZX3eRUWFpZZp63PxxDVr19fo/1fvz4ln5FMJsOnn35a\n5Tisra3h4+ODiIgIPH/+HBYWFggLC4NUKlWq47h//z769OkDKysrfP7553jjjTdgYWEBgJcivJ5w\nVNHks2eMoUOHDvj666/LPZ6NjU2l59QXSjACMjExwdChQzF06FAAwOHDh+Hr64uvv/4aGzZsACBM\nM9o9e/aguLgYhw8fVuqv8Pz5czx+/Fjt45R8G23SpEmlT3TaUHJzV9WqKSEhQWkfQPNrq8n7Kam8\nvnHjRpkWQiWxlFZSKZ+VlQUrKyvF+pInodK09fm8rqLr0aZNGxQVFSEhIQHu7u5K2xISEiASiSr9\nslCd32Vra2uV7+3162NrawsrKytkZ2dX+3du0qRJiIiIwE8//QRXV1ckJSVh1apVSvuEh4crirpf\n7zuUmZmpVhNra2trlcWLqj57V1dX/P333+jXr1+NaGJPdTACUdXKrHPnzgCUizgaNGhQ4U1D1S9Z\ndX/xSp4WXv8mHBwcrPJpoUGDBir/QEaPHg1TU1MsWbKkTDEOwFv75OfnVyvW0jw8PODg4IDQ0FA8\nfPhQsb6goABr1qyBWCxW6gxa2bV9nSbvRyaTwdzcHBs3bkReXp5in9TUVOzZs6fMZ1TSWqp0mT/A\nW4u9TtPPR10lRbOqPsuSZrkrV65UWn/t2jX8/PPP6N27d6XFYxUdvzKurq44f/680rV88uQJQkND\nla6lWCzGuHHjEBsbi4MHD6o8VkZGhlrn9PX1hY2NDXbt2oVdu3ZBLBZj/PjxSvuU91ls3bpV6Xew\nIm+88QZu3ryJBw8eKNa9fPkSGzduLLPvxIkTkZ6eXu4TjLrn1Bd6gtGBuLg47N69W+U2Pz8/WFhY\nKPpR9O7dG/b29sjKylK06Z8wYYJi/x49euDQoUOYNWsWevToASMjIwwYMEBRBqvqhlLdIqORI0fi\nm2++gY+PD6ZPnw5jY2McPXoU8fHxsLGxKXP87t2749ixY1i9ejXs7e0hEokwZswYSCQSbNq0CVOn\nTkXbtm0xYcIEODg4ICMjA/Hx8YiMjMSNGzcqbDKqCbFYjA0bNsDPzw/dunXD9OnT0aBBA+zbtw8x\nMTFYvHixUrPYyq7t6zR5P1ZWVli+fDk++eQT9OzZExMmTEBubi42b94MV1dXXL58WenY/v7+WLRo\nEaZPn46bN2+icePGiIqKUnkz1vTzUVeTJk3g7OyMvXv3wsnJCU2bNoWFhQWGDRuGgQMHYvTo0di7\ndy+ePHkCX19fpKenY+PGjahfvz7Wr19f6fE9PT0hFouxYsUKPH78GBYWFmjTpg08PT0rfe2sWbMw\nfvx49O/fH+PHj0dWVha+//57ODo6lrmprlixAufOncPo0aMxevRoeHl5wcTEBPfu3cNvv/2Grl27\nIjQ0tNJz1qtXD/7+/tiwYQPi4uIgk8nKFI36+PhgwYIFmDBhAmbNmgUrKyucO3cOhw8fhpOTU7nF\nXK+/t71792LgwIEIDAxEfn4+du/erbK4cM6cOTh69CjmzZuHEydOoF+/frC0tMT9+/dx/PhxmJub\n48SJE5W+N73Rd7O12iw6OpqJRCImFotVNlMWi8Xszp07jDHGtm7dymQyGWvevDkzMTFhLVq0YL6+\nviw6OlrpmLm5uey9995jzZo1Y0ZGRkrNWwMCAphYLFbaX9U6xniTR7FYrNTksYSjo2OZZpIRERHM\nw8ODWVhYMBsbG+bv78/u37+vct/ExEQ2aNAgZmlpqXifpZ07d475+fmxpk2bMhMTE9ayZUvWv39/\n9vXXX7MXL15UeE0rajobFBTExGJxmebLp06dYjKZjFlaWjIzMzPWpUsXtn379jKvr+jaVkST97N5\n82b2xhtvMFNTU+bi4sLWrVtXbnPdmJgY1qtXL2ZmZsZsbGxYYGAgy8rKUtlEV5PPR9W68sTGxrJe\nvXoxCwsLJhKJlJrUFhYWslWrVrG2bdsyU1NT1qRJE+bn58euXbum1rEZY2znzp1MKpUyExMTpfcV\nGhpa6fVfs2YNa9WqFTM1NWVSqZSFhoayHTt2qHxdbm4uW758OXN3d2fm5uasYcOGTCqVsunTp7PY\n2Fi1442Li1P8Tu/Zs0flPqdPn2a9e/dmDRs2ZFZWVuytt95i169fZ97e3mWaJKtaV3Jd3njjDWZi\nYsLatGnD1qxZw06cOFGmmTJj/HNYv34969atG7OwsGAWFhbM1dWVjR8/nh09elTt96YPIsb0VEOq\nRd9++y1CQkJgZGQEX19fRbnoypUrsX37dhgZGWH9+vUYNGiQwJESUtaOHTswZcoUREdHo0+fPkKH\nQ4jO1LgispMnT+Lnn3/G1atXYWxsrChPTUhIwL59+5CQkAC5XI6BAwfi1q1bNJw9IYQIpMbdfTdt\n2oSFCxcqhr4vKS+PjIyEv78/jI2N4ejoCGdnZ510OCOEEKKeGpdgEhMTcfr0aXTv3h3e3t64ePEi\nAODBgwews7NT7GdnZwe5XC5UmIRUqCY0MSWkugyyiEwmkyE9Pb3M+hUrVqCwsBBPnjzBhQsX8Oef\nf2L06NEq24sD9EdMDFNAQIDSYImE1FYGmWBe7wtQ2qZNmxQjvXbr1g1isVgxNEJKSopiv9TUVJXD\nJTg7O+POnTvaD5oQQmoxJyenMmMWVqbGFZGNGDFC0c771q1byM/Ph42NDYYPH469e/ciPz8fycnJ\nSExMVNm+/s6dO2CM0cIYlixZIngMhrLQtaBrQdei4qUqX8wN8gmmIlOmTMGUKVPg7u4OExMT7Nq1\nCwAfcG/06NGQSqWoV68eQkJCqIiMEEIEVOMSjLGxMcLCwlRuW7RoERYtWqTniAghhKhS44rIiPZ4\ne3sLHYLBoGvxCl2LV+haVE+N7MlfHSKRCHXsLRNCSLVV5d5JTzCEEEJ0ghIMIYQQnahxlfyE1CZB\nAQHA3btK6x49foS43Ecw6e0KU5EpZo+dDV+ZryDxEVIdlGAIEdLduwg6darMau9WwKnWaQCAOxt5\n/wNKMqSmoQRDiI4VFwPPnr1anj8HcnOBFy8AdSZ3vNP5Dhb/51u8zPGFiQlgZsYXc3PAwgJo0IAv\nDRsC/0ywSIhBoARDiIYYAzIzAbmcL+npwMOH/N+MDJ40Hj0CHj8GsrKAp095MmjYkCcCCwugfn2e\nJKzUHI9VnvECYWFAfj5PTC9f8iT1/Lly8mrQALCyAho3BmxsgCZNAFtboFkzoHlzvrRsCdjZ8fU0\nmwXRJUowhKjw+DFw+zZw5w7/9+5dvty7B6Sm8gRhZwdIJPym3awZ0Lo14On56sZubc1v9JaWQL1y\n/tKCvAGULSErw8PdDOHbK96nqAjIyQGePOFLZiZPdBkZPAFeuMCTYElifPqUx+/o+Gpxdn61WFtr\ncMEIUYESDKmzGAPS0oBr1/iSkADcvAn873/8ScHZGXBy4ouXF/Duu0CrVoC9PU8w+uJ0yQkfzvqw\n0v2MjPjTi5UVT3aVefECSEnhSfPuXSApCYiM5Ek1MREwMQHc3PgilQLu7nxp1qz674nUDdTRktQJ\nxcU8cVy8CPz116vFyAho354vUumrG2qzZoA+hrIrrxXZpdzHMO7tAjOxGT70/1DvFfyM8aeemzf5\ncv06EB/PFyMjoFMnoHNnoEsXoGtXoE0b/VwvIpyq3DspwZBa6e+/gT/+AM6fB2JigEuXeJ1D1678\nxti5M9CxIy/eIupjDHjwgCfnS5f4cvEirw/y9ORPer16Ad278zonUntQglEDJZjahzH+EHD6NHDq\nFP83MxPo0YMvXl48sTRpInSktdeDB8Cff/Kkfu4ccPky4OIC9O0LeHsDffpQnU5NRwlGDZRgaoeM\nDOD4ceDYMf5vXh6/mfXtC7z5JtCuHbWQEtLLl/zp5tQp4ORJnnicnYFBg/jSqxdgaip0lEQTlGDU\nQAmmZioq4kVdhw8DUVG8ErpvX2DgQL64uVEdgCHLzwdiY4EjR4CjR3mdjrc34OvLFzs7oSMklaEE\nowZKMDXH8+f8hvTzz8ChQ7y+xMcHGDqUF32ZmAgdIamqx4+B338Hfv2V/9uqFTBiBODnx58+6cuC\n4akTCSY2NhazZs1CQUGBYubKbt26AQBWrlyJ7du3w8jICOvXr8egQYPKvJ4SjGHLyeHJ5MAB/k23\nWzfg7beBYcN4Pw1S+xQWAmfPAuHhQEQE74A6ejRvFt6+vdDRkRJ1IsF4e3tj4cKFGDx4MA4fPozV\nq1fj5MmTSEhIwNixY/Hnn39CLpdj4MCBuHXrFsSvFcRTgjE8L17woq89e/gTS+/ewDvvAMOHU8V8\nXcMYb5W2bx/w00+8Jdr48Xyxtxc6urqtTswH06JFC2RnZwMAsrKyIJFIAACRkZHw9/eHsbExHB0d\n4ezsjNjYWCFDJRVgjFf8BgbyoUu+/RYYPBhITuZPMJMnU3Kpi0Qi/tS6di1vGbh5M+8I2qkTMGAA\nEBbGm0STmqHG9eT/8ssv0bt3b3zyyScoLi7GH3/8AQB48OABunfvrtjPzs4OcrmaAz0RvXn4EAgN\n5YtIBAQEAFevUiUvKUss5k+zvXsD33zD62u2bwc++ggYMwaYNo0nHmK4DDLByGQypKenl1m/YsUK\nrF+/HuvXr4efnx/279+PKVOm4OjRoyqPIyqnpjAoKEjxf29vb5p3W8cYA06c4N9Gjx4FRo4Eduzg\nnfGoMpeow8yMF5u+8w4f3iY0lBeh2tsDH3zA11OjD+2Kjo5GdHR0tY5R4+pgLC0t8fTpUwAAYwxW\nVlbIzs7Gl19+CQBYsGABAGDIkCFYunQpvLy8lF5PdTD68+wZL9L49ls+vMiMGcC4cUCjRkJHRmqD\nwkLgl1+AjRt5s+cZM4CZM/lgo0T76kQdjLOzM079M0HTiRMn4OrqCgAYPnw49u7di/z8fCQnJyMx\nMRGenp5ChlpnpaYC8+bxpqfHjgEhIbwYbOZMSi5Ee+rV482aSzrbpqTw0QPef5/3kyLCM8gisops\n2bIFH3zwAV6+fAlzc3Ns2bIFACCVSjF69GhIpVJF8+XyisiIbly/DqxZw/utTJwIxMVR02KiH1Ip\nsHUr8MUX/ImmZ0/eAXfRIj4CNBFGjSsiqy4qItO+y5f5H/a5c8Ds2byoonFjoaMidVlODvDdd8DX\nX/Ox6IKCqEFAddWJIjJiOC5f5hWtb73Fx/9KSuLfGCm5EKE1bMiLaZOSgH79+OgP777Lpx4g+kMJ\nhmjs5k3e09rHB5DJ+ARVH32k30m4CFGHuTkwZw6flbRzZ/5F6L33+OjPRPcowRC1paUBU6fyP1IP\nD/5H++GHvAkpIYbMwgJYsIBX/tva8nqZJUt4S0eiO5RgSKWePweWLuXjQllb8z/STz/lf7SE1CRW\nVsCXX/KpBG7fBt54A/jhB95Xi2gfVfKTcjEG7N3Ly7LffBNYuZJahZHa5Y8/gFmz+JelDRuADh2E\njshw1YnBLquLEox6rl3jf3jZ2fwPr1cvoSMiRDeKingT588/ByZN4k/rVJ9YFrUiI9WWm8uLv/r3\nB/7v//jItpRcSG1mZMQ7Z167BsjlvCi4nNGniIYowRCFY8d4EcH9+0B8PB/jychI6KgI0Y+mTfmU\nERs38oE0p0wB/hmVilQRJRiC7Gz+x/Tee8D69cCPPwLNmgkdFSHCGDqUf8EyNuZfuE6eFDqimosS\nTB134gT/IzIx4UUEPj5CR0SI8Bo25KN/h4Twyc4+/hh4+VLoqGoequSvo/LyeL+AgweB778HhgwR\nOiJCDNOjR8D06XwCtL17+YCadRFV8hO13LjBx2dKT+ejHFNyIaR8TZoABw7wIuSePXk9DVEPPcHU\nIYzxGQEXLOB9Wt57jyb8IkQTf/3FxzTr3x9Yt65uTXJG/WDUUFcTTG4ub4p5+TKwbx8f3pwQormn\nT4EJE4DMTP5k06KF0BHpBxWREZWSkvijfXExEBNDyYWQ6rC0BMLDgcGDgW7dgAsXhI7IcFGCqeV+\n/x3o0YMXh4WFUQ9lQrRBLOY9/zdt4lNW/PST0BEZJoNMMPv370e7du1gZGSES5cuKW1buXIlXFxc\n4ObmhiNHjijWx8XFwd3dHS4uLpgzZ46+QzY4jPE+LQEBwP79fNRjqm8hRLuGDeO9/j/5hA+iWQdL\n3ytkkAnG3d0d4eHh6NOnj9L6hIQE7Nu3DwkJCYiKisLMmTMVZYIzZszAtm3bkJiYiMTERERFRQkR\nukEoLOTjiG3eDJw/D7x2GQkhWtSxIx80c98+3py5sFDoiAyHQSYYNzc3uLq6llkfGRkJf39/GBsb\nw9HREc7OzoiJiUFaWhpycnLg6ekJAJg4cSIiIiL0HbZBePqUzzB5+zZPLq1bCx0RIbWfRAKcPg2k\npPAx/F68EDoiw2CQCaY8Dx48gJ2dneJnOzs7yOXyMuslEgnkcrkQIQrq77/59LCtWgGHDgGNGgkd\nESF1R8OGwM8/A/Xq8S95NJkZUE+oE8tkMqSnp5dZHxwcjGHDhun03EFBQYr/e3t7w9vbW6fn04e7\nd4FBg4AxY/hw41TfQoj+mZjw3v6BgXw68UOH+CR9NVF0dDSio6OrdQzBEszRKoyHLZFIkJKSovg5\nNTUVdnZ2kEgkSE1NVVovkUjKPU7pBFMbXLvGB+ibP59X5hNChGNkxOeX+de/+Je+Y8f4TJo1zetf\nvpcuXarxMQy+iKx0x57hw4dj7969yM/PR3JyMhITE+Hp6YnmzZvD0tISMTExYIwhLCwMI0aMEDBq\n/blyhX9TWrWKkgshhkIkAr76Cujdm/eXyc4WOiJhGGSCCQ8Ph729PS5cuABfX18MHToUACCVSjF6\n9GhIpVIMHToUISEhEP1TFhQSEoKpU6fCxcUFzs7OGFIHBtj66y/+y7t+PTB2rNDREEJKE4mA//yH\nd8YcOhTIyRE6Iv2joWJqqEuX+ND6GzcCo0YJHQ0hpDzFxcDMmUBCAu/4bG4udERVQ2ORqaE2JJhr\n14CBA/lcFSNHCh0NIaQyxcV8Xpm8PD5+WU2cKZbGIqsDkpL48Pr/+Q8lF0JqCrEY2LGDN13+4IO6\n0+OfEkwNkpbGK/QXLwb8/YWOhhCiCRMTPsFfbCzwxRdCR6MfgjVTJpp58oRX6E+ZAsyYIXQ0hJCq\nsLQEfvuNj27eujUvNqvNqA6mBsjP58Vi7u7AN99QJ0pCarpr1/ikZb/+CvwzwpXBo0p+NdS0BMMY\nMG0aHwYmPLxmVg4SQsqKjOT1MbGxQMuWQkdTuarcO6mIzMCtXQvExQFnzlByIaQ2eftt/iTj5wec\nOgWYmQkdkfbRE4wBCw/nvfMvXABKjeVJCKklGOPjBzZoAGzbJnQ0FaNmyrXIjRt8bomICEouhNRW\nIhFPLGfPArt3Cx2N9tETjAHKyeEVf/Pm8VZjhJDa7coV3nn63DlAxVRYBoEq+dVg6Amm5JHZ0pKP\nyEoIqRvNdBceAAAgAElEQVQ2beKz0F64YJj1MZRg1GDoCeabb4CwMP5NxhB/yQghusEYMHo00LQp\nH2PQ0FCCUYMhJ5jYWGDYMP4NhqY6JqTuycoCOnYEtmzhHasNCVXy12DPnvFevRs3UnIhpK6ysgK+\n/573fasNc8jQE4yBCAwEXr7kA+IRQuq26dN5kZkh1cNqtYhs586dism8NDFx4kSNX6NPhphgfv4Z\n+OgjPoGYpaXQ0RBChPb0KR8aypCKyqp072TlEIlEVVq04aeffmJSqZSJxWJ28eJFxfojR44wDw8P\n5u7uzjw8PNiJEycU2y5evMjat2/PnJ2d2ezZs8s9dgVvWRBpaYw1b87YmTNCR0IIMSS//86YvT1j\nWVlCR8JV5d5Z7lAxJ06cUPq5oKAAn376KR4/foz3338fbdu2BQAkJCTgu+++g42NDVatWqVpUlTJ\n3d0d4eHhCAwMVHqKsrW1xa+//ormzZvj+vXrGDx4MFJTUwEAM2bMwLZt2+Dp6QkfHx9ERUXViGmT\nZ80CJk/mc3cTQkiJQYP4smQJb11aI6mbiT777DPWtm1blp2dXWZbdnY2c3NzY59//rnGGa4i3t7e\nLC4uTuW24uJiZm1tzfLz89mDBw+Ym5ubYtuPP/7IAgMDVb5Og7escxERjLm6MpaXJ3QkhBBDlJHB\nmK0tY1evCh1J1e6darci27FjByZPngxLFZUElpaWmDJlCnbosYb64MGD8PDwgLGxMeRyOexKjaci\nkUggl8v1FktVPH3Kn162bKH+LoQQ1WxsgKAgPiahgVUdq0Xt0ZQzMjJQVFRU7vaioiI8fPhQ7RPL\nZDKkp6eXWR8cHIxhw4ZV+Nrr169jwYIFOHr0qNrnKy0oKEjxf29vb3h7e1fpONWxeDGvvOvbV++n\nJoTUIIGBvDXZvn18lA99iY6ORnR0dPUOou6jTqdOnZiTkxN79OhRmW2ZmZmsTZs2rFOnTho/QlVE\nVRFZSkoKc3V1ZefPn1ese72IbM+ePQZdRPbHH7xiX8WlJISQMs6eZUwiYSwnR7gYqnLvVPsJJigo\nCH5+fnBzc8PkyZPh5uYGALhx4wZCQ0Px+PFjHDhwoHrZTnUCVPw/KysLvr6+WLVqFXr06KFY36JF\nC1haWiImJgaenp4ICwvD7NmztR6LNhQV8W8kX38NWFsLHQ0hpCbo1YvPgPnFF8CXXwodjQY0yUb/\n/e9/mb29fZmmyfb29uzAgQMaZ7eKzmNnZ8fMzMxYs2bN2JAhQxhjjC1fvpxZWFiwTp06KZaMjAzG\n2Ktmyk5OTuzDDz8s99gavmWt27yZsb59GSsuFjQMQkgNI5cz1rgxY6mpwpy/KvdOjXvyFxUVIS4u\nDklJSQCANm3aoGvXrhCLa8aoM0J2tHz6FHjjDeDQIaBLF0FCIITUYPPn8/vId9/p/9w02KUahEww\nCxcC6elAaKggpyeE1HCPHvH5YmJiAGdn/Z5bLwnm9u3biIyMRHJyMgD+BPP222/DyclJoxMLRagE\nc/cu4OEBxMcDLVvq/fSEkFpi2TLg1i39z4Cp8wTz73//G19++SWKi4uV1ovFYixcuBDLly/X6ORC\nECrBvPsu0K4d8Pnnej81IaQWycnhTy/HjvHxyvRFp8P1b9++HcHBwejevTsiIiJw69Yt3Lp1CxER\nEejRowdWrFiBUCr7USkmBjh/HvjkE6EjIYTUdA0bAgsWAJ99JnQklVP7Caak1/yZM2dgbGystK2g\noAB9+vRBfn4+4uLidBKotgjxBDN0KDBiBG+eTAgh1fXiBX+K+eUXoHNn/ZxTp08wN27cgL+/f5nk\nAgDGxsZ49913kZCQoNHJ64ILF4CEBD6gJSGEaIOZGTBnDrB2rdCRVEztjpYmJibIyckpd/uzZ89g\nYmKilaBqk6VLgUWLALo0hBBtmj4daNMGuHcPaNVK6GhUU/sJplu3btiyZYvK8cMePnyILVu2wMvL\nS6vB1XQxMfT0QgjRjUaN+L1l3TqhIymf2nUwp0+fRv/+/RUjJ7dr1w4AcO3aNYSGhiInJwfHjx9H\nnz59dBpwdemzDsbHBxg+HHj/fb2cjhBSx6SkAB07AklJgJWVbs+l82bKv/zyC2bNmoWUlBSl9Q4O\nDtiwYQPeeustjU4uBH0lmJgY4P/+D0hMBExNdX46QkgdNXEiIJXylmW6pJeOliVDxZTuaOnh4UFD\nxbzm7bf5cPwzZ+r8VISQOuzqVWDIECA5WbdfZmmoGDXoI8HcuQN0784r3+rX1+mpCCEEgwcD/v5A\nQIDuzkFDxahBHwnmo494M8IaNaw2IaTG+vVXIDiYd+jWFRoqRg26TjBPnwKOjsCVK4C9vc5OQwgh\nCoWFvKnykSN8SCpdoKFiDMCOHYBMRsmFEKI/9erx4rFt24SO5DXqThzTpUsX5uXlxfLz88tsy8/P\nZ927d2ddunRR93AV+umnn5hUKmVisbjMlMmMMXbv3j1mYWHB1q5dq1hXMuGYs7Mzmz17drnH1uAt\na6yoiDEnJ8bOndPZKQghRKXbtxmztWXsxQvdHL8q906DHCrG3d0d4eHh5fapmTt3Lnx9fZXWzZgx\nA9u2bUNiYiISExMRFRWllVg0cegQ0LgxUGo2Z0II0QsnJ6B9e+Dnn4WO5BW1E4w+h4pxc3ODq6ur\nym0RERFo06YNpFKpYl1aWhpycnLg6ekJAJg4cSIiIiK0Eosm1q3j4wOJRHo/NSGEYOpU4PvvhY7i\nlRo1VMyzZ8+wevVqBAUFKa2Xy+Wws7NT/CyRSCCXy3Uay+tu3+aTiY0erdfTEkKIwsiRQFwc7yJh\nCNQe7PKzzz5D//79IZVKyx0qZrcGU6zJZDKVySo4OBjDhg1T+ZqgoCB8/PHHqF+/vmDTHpdn1y5g\n7Fga1JIQIhwzM94fJjQUeO17uCDUTjB9+vRBeHg4Zs2aha+//lppm4ODA3bt2qXROGRHjx5VP8p/\nxMbG4uDBg5g/fz6ysrIgFothbm6OkSNHIjU1VbFfamoqJBJJuccp/QTk7e0Nb29vjWMprbgY2LkT\niIys1mEIIaTa3nuPzz+1ZEn1iuujo6MRHR1drVgMeqiYfv36Ye3atfDw8CizbenSpWjYsCHmzp0L\nAPDy8sL69evh6ekJX19fzJ49G0OGDCnzOl30gzl5kneuvHJFq4clhBCNMQa4uQE//AB07aq941bl\n3qn2E0wJIyMjeHp6KirUdSE8PByzZ89GZmYmfH190blzZxw+fLjC14SEhCAgIAB5eXnw8fFRmVx0\nZedOYNIkvZ2OEELKJRIBo0YBBw5oN8FUKRZNn2AAIDc3F48ePVKZzRwcHLQSmK5o+wnm2TPAzg74\n3/+AZs20dlhCCKmyuDhgzBjg1i3ttWrV6RNMUVERVq9ejW+//VZl5XxJAEVFRRoFUNMdPAi8+SYl\nF0KI4ejShQ8fEx8PdOggXBxqJ5iFCxdi7dq1aNeuHUaNGoUmTZqU2UdUBzuA7NxJQ/ITQgyLSMSb\nLB88KGyCUbuIrGXLlujYsWOldSGGTptFZHfv8jJOuZwmFSOEGJY//uAdL69f187xdDrY5ZMnTzBi\nxAiNg6rN9u0D3nmHkgshxPB4eQFZWcDNm8LFoHaCad++PdLS0nQZS40TGQn4+QkdBSGElCUWvyom\nEywGdXdcsmQJNm3ahPv37+synhojPR24cQPo10/oSAghRLVRo4RNMOVW8i9dulSp0p4xBkdHR7Rr\n1w4jRoxAmzZtYGRkVOZ1n3/+uW4iNTC//MLnwaahYQghhurNN4HUVCApCWjTRv/nL7eSv6o981+f\n7dLQaKuS39cXmDCBtzUnhBBDNWkS0L07MGNG9Y6j1SmT7969W6UgHB0dq/Q6fdFGgsnJASQSICUF\naNRIS4ERQogO7NrF54g5cKB6x9FqgqmttJFgDhwAtm4Ffv9dS0ERQoiOyOW8L0xGBq/4ryqdNlMm\nr0RG8tFKCSHE0EkkgK0t8Ndf+j93uU8wJZX8ixcvhpGRUZlK//IYeiV/dZ9gCgqA5s2Bq1f5B0cI\nIYZu1iygVStg3ryqH0OrRWQllfwvXryAiYmJ2pX+tb2S/8QJYMECIDZWi0ERQogOhYcD331XvWJ9\nrQ52mZSUBAAw+acdbsnPdV1EBBWPEUJqFm9vYOJE4OVL/Y48QpX8GnJ1BfbvBzp21GJQhBCiY926\nAWvXAn37Vu31VMmvY3I58Pgx4O4udCSEEKKZAQOA48f1e85yi8h27txZpeH3J06cWK2AAGD//v0I\nCgrCzZs38eeff6JLly6KbVevXkVgYCBycnIgFotx8eJFmJiYIC4uDgEBAXjx4gV8fHywbt26asfx\nuuhonv11MDs0IYTo1IABQFAQsGyZ/s6p1Z782ppw7ObNmxCLxQgMDMRXX32lSDCFhYXw8PDA7t27\n4e7ujidPnqBRo0YQi8Xw9PTEhg0b4OnpCR8fH8yePVvltMnVKSKbNo23J//ww2q9PUII0bu8PN5c\n+cEDwNJS89drtZL/xIkTmkegJW5ubirXHzlyBB06dID7P2VUjRs3BgCkpaUhJycHnp6eAPhTVERE\nhMoEUx3R0cDs2Vo9JCGE6IW5OeDpCZw+Dbz1ln7OWW6C8fb21k8EGkhMTIRIJMKQIUOQkZGBMWPG\nYN68eZDL5bCzs1PsJ5FIIJfLtXru1FQ+t0K7dlo9LCGE6M2AAbyrheAJRtdkMhnS09PLrA8ODsaw\nYcNUvqagoABnz57FxYsXYW5ujgEDBsDDwwON9DAgGNW/EEJqup49gX//W3/n0yjB3L9/H0uWLMGR\nI0fw999/4/fff0f//v2RkZGB+fPnY+bMmejWrZtaxzp69KjGwdrb26NPnz6wtrYGAPj4+ODSpUsY\nP348UlNTFfulpqZCUkE3+6CgIMX/vb291Xpai47mbckJIaSm6tIFuHIFKCwE6lVy94+OjkZ0dHT1\nTsjUlJSUxGxtbZmlpSWTyWRMJBKx48ePK7a3b9+eTZ8+Xd3DqcXb25tdvHhR8fOTJ09Yly5dWG5u\nLisoKGADBw5kv/32G2OMMU9PT3bhwgVWXFzMhg4dyg4fPqzymBq8ZSVt2jB27VqVXkoIIQbDxYWx\n+HjNX1eVe6faBT6LFy+GWCxGfHw89uzZU2a7j48Pzp49W71s94/w8HDY29vjwoUL8PX1xdChQwEA\nVlZWmDt3Lrp164bOnTvDw8NDsS0kJARTp06Fi4sLnJ2dtVrBf/8+H6JfKtXaIQkhRBAeHkBcnH7O\npXYR2bFjxzBr1iw4ODggMzOzzPZWrVohJSVFK0H5+fnBr5zJ7seNG4dx48aVWe/h4YH4+HitnP91\np07x+pcqdAsihBCD0rUrcPEin4hM19R+gnn69ClatmxZ7vb8/HwUFhZqJShDc/Ik1b8QQmoHfT7B\nqJ1g7OzscP369XK3x8TEwNnZWStBGZroaKBfP6GjIISQ6itd0a9raieYUaNGYdu2bYiPjy8zhMzB\ngwfx008/YfTo0VoPUGgpKcCzZ0DbtkJHQggh1WdpCdjZATdu6P5cao+mnJ2djZ49e+Lu3bvo06cP\nfv/9d8hkMmRnZyM2NhadOnXCuXPnYG5uruuYq0XT4Q4iIvj0yIcO6TAoQgjRo3HjgIEDgcmT1X+N\nTkdTbtSoEc6fP4+pU6fizz//BMD7sty6dQsffPABoqOjDT65VMXly0DnzkJHQQgh2qOvehi1n2Be\nvnwJ039mqmGMISMjA4wx2NraKgbGlMvlFXZwNASaZuFhw4CAAGDUKN3FRAgh+nT6NDB/PnDhgvqv\n0ekTzIQJE5RO1LRpUzRr1kyRXNLS0tC/f3+NTl4T0BMMIaS26dwZiI/XfUW/2gkmPDwcs8sZSvjh\nw4fo37+/yv4xNVlGBq/gb91a6EgIIUR7GjYEHByAhATdnkftBPP9999jw4YN+PLLL5XWZ2RkYMCA\nAUhPT8eRI0e0HqCQSp5eqIMlIaS28fDgHS51Se2e/JMmTYJcLseiRYsgkUgwYcIEPH78GAMHDkRK\nSgqOHDkCDw8PXcaqd5cuUfEYIaR26tqVV/RPmaK7c2g0mvKiRYsgl8vx3nvvwdjYGKtXr0ZSUhKi\noqLg5eWlqxgFc/kyr+QnhJDaxsMD+PFH3Z5D7VZkJYqLi/HOO+8gIiIC9evXx6FDh9C3b19dxad1\nmrSEcHUFwsNpkjFCSO3z7BnQtCkfyNfIqPL9q9KKrNwEs3PnzjI99kvk5eVh/vz58PPzK9NybOLE\niRoFoG/qXqSnT4EWLYDs7MrnTSCEkJrIwYE3WXZ0rHxfrSYYcRWmbhSJRCgqKtL4dfqk7kU6cwb4\n5BMgJkYPQRFCiAD69wcWLeK9+itTlQRT7nfzEydOaHSg2ubyZT4oHCGE1FbOzkBionoJpirKTTDq\nTCNcm126xOevJoSQ2srZGbh9W3fH17wcTA/279+Pdu3awcjICJcuXVKsf/HiBfz9/dGhQwdIpVKl\nPjlxcXFwd3eHi4sL5syZU+0YqAc/IaS203WCKfcJpqSSf/z48RCLxRVW+pemjUp+d3d3hIeHIzAw\nUGn93r17AQBXr15FXl4epFIpxo4dCwcHB8yYMQPbtm2Dp6cnfHx8EBUVVeVpk1+84I+N7u7VfiuE\nEGKwBEswkydPhkgkwpgxY2BiYoLJao7rrI0E4+bmpnJ9ixYt8Pz5cxQVFeH58+cwMTGBpaUl0tLS\nkJOTA09PT0UMERERVU4w167xC29mVuW3QAghBs/JCUhKAoqLgSq066pUpZX8xsbGSj9XRJ0nnOoY\nPHgwwsLC0KJFC+Tm5uKbb76BlZUVbt++DTs7O8V+EokEcrm8yueh4jFCSF1gYQFYWwNyOWBvr/3j\nq13JX1ml/+bNm7Fu3TokqDl6mkwmQ3p6epn1wcHBGFZO9/ndu3cjLy8PaWlpePz4Md58800MGDBA\nrfOVFhQUpPi/t7d3mfd2/ToVjxFC6oaSYrLXE0x0dDSio6OrdWytdSHMyMjAzZs31d7/6NGjGp/j\n/Pnz8PPzg5GREWxtbdGrVy/ExcWhd+/eSE1NVeyXmppa4bw0pROMKnfuAHW8ER0hpI4oSTD9+imv\nf/3L99KlSzU+tkG2IiutdMceNzc3RVHd8+fPceHCBbi5uaF58+awtLRETEwMGGMICwvDiBEjqnzO\nO3eANm2qHTohhBg8XVb0G2SCCQ8Ph729PS5cuABfX18MHToUABAYGIj8/Hy4u7vD09MTU6ZMQfv2\n7QEAISEhmDp1KlxcXODs7FzlCv7iYiA5mRIMIaRuKOlsqQsGOcqWn58f/Pz8yqw3NTXF7t27Vb7G\nw8MD8fHx1T53WhrQqBHQoEG1D0UIIQavzj3BCCkpiZ5eCCF1h5MTrxbQcJgxtVT4BPPVV1+p3fT4\n/PnzOm+mrA937vALTgghdYGlJS+xSUsDWrbU7rErTDDz5s3T7tlqAEowhJC6pqSYTK8JRtMRlWvD\nE0xSEjB4sNBREEKI/pQkmD59tHvcChNMXRxRmZ5gCCF1ja4q+qmS/zXUB4YQUtdQgtGDp0+B3Fyg\neXOhIyGEEP1xcaEEo3MlTZRrQVUSIYSozcmJd7bUdlNlSjClUB8YQkhd1LgxYGoK/P23do9LCaYU\nquAnhNRVuqiHoQRTClXwE0LqqpIe/dpECaaUpCR6giGE1E3NmlERmU5RERkhpK6ytQUyMrR7TEow\n/ygoAFJTgVathI6EEEL0jxKMDqWk8P4vpqZCR0IIIfpHCUaHqHiMEFKX1ZkEM2/ePLRt2xYdO3bE\nyJEjkZ2drdi2cuVKuLi4wM3NDUeOHFGsj4uLg7u7O1xcXDBnzhyNz0l9YAghdVmdSTCDBg3C9evX\nceXKFbi6umLlypUAgISEBOzbtw8JCQmIiorCzJkzwf7pejpjxgxs27YNiYmJSExMRFRUlEbnpCcY\nQkhdVmcSjEwmg1jMQ/Py8kJqaioAIDIyEv7+/jA2NoajoyOcnZ0RExODtLQ05OTkwNPTEwAwceJE\nREREaHRO6gNDCKnLGjUC8vKAly+1d0yDTDClbd++HT4+PgCABw8ewM7OTrHNzs4Ocrm8zHqJRAK5\nXK7ReagPDCGkLhOJABsbIDNTe8escD4YXZLJZEhPTy+zPjg4GMOGDQMArFixAiYmJhg7dqxWzx0U\nFKT4v7e3N7y9vXH/PuDgoNXTEEJIjVJSTCaRANHR0YiOjq7W8QRLMEePHq1w+44dO/Dbb7/h+PHj\ninUSiQQpKSmKn1NTU2FnZweJRKIoRitZL5FIyj126QQDAIWFQHY2YG2t4ZsghJBapHQ9TMmX7xJL\nly7V+HgGWUQWFRWFNWvWIDIyEmZmZor1w4cPx969e5Gfn4/k5GQkJibC09MTzZs3h6WlJWJiYsAY\nQ1hYGEaMGKH2+Z484eWP9QRLt4QQIjxtV/Qb5C31ww8/RH5+PmQyGQCgR48eCAkJgVQqxejRoyGV\nSlGvXj2EhIRA9M/kLSEhIQgICEBeXh58fHwwZMgQtc+XmcnLHgkhpC7TdoIRMabtKWYMm0gkwutv\n+cwZYOFC4OxZgYIihBADsGwZkJ8PfPFF2W2q7p2VMcgiMn2jJxhCCNH+EwwlGFCCIYQQQPvNlCnB\ngBIMIYQA9ASjE5RgCCGEEoxOUIIhhBBKMDpBCYYQQoAmTYCsLKCoSDvHowQDSjCEEAIARkaAlRXw\n6JF2jkcJBpRgCCGkhDaLySjBgBIMIYSUsLXVXlPlOp9g8vOB3Fw+FhkhhNR1Njb0BKM1jx7xiq1/\nhjQjhJA6jYrItIiKxwgh5BVKMFpECYYQQl6hBKNFlGAIIeQVSjBaRAmGEEJeoQSjRZRgCCHklVqf\nYObNm4e2bduiY8eOGDlyJLKzswEAR48eRdeuXdGhQwd07doVJ0+eVLwmLi4O7u7ucHFxwZw5c9Q+\nFyUYQgh5pdb3gxk0aBCuX7+OK1euwNXVFStXrgQA2Nra4tdff8XVq1exc+dOTJgwQfGaGTNmYNu2\nbUhMTERiYiKioqLUOhclGEIIeaVkThhtzHVskAlGJpNBLOaheXl5ITU1FQDQqVMnNG/eHAAglUqR\nl5eHgoICpKWlIScnB56engCAiRMnIiIiQq1zUYIhhJBXTE0BMzPgn4KjajHIBFPa9u3b4ePjU2b9\nwYMH4eHhAWNjY8jlctjZ2Sm2SSQSyOVytY5PCYYQQpRpqx6mXvUPUTUymQzp6ell1gcHB2PYsGEA\ngBUrVsDExARjx45V2uf69etYsGABjh49WqVzBwUFKf6fmuoNGxvvKh2HEEJqI1tb4Pffo/HDD9HV\nOo6IMW2UtGnfjh07sHXrVhw/fhxmZmaK9ampqRgwYAB27NiBHj16AADS0tLQv39/3LhxAwDw448/\n4tSpU/juu+/KHFckEqH0W7awAB4+BBo00PEbIoSQGmLYMGDqVODtt1+te/3eqQ6DLCKLiorCmjVr\nEBkZqZRcsrKy4Ovri1WrVimSCwC0aNEClpaWiImJAWMMYWFhGDFiRKXnyc3lE+tYWOjkbRBCSI2k\nrSIyg0wwH374IZ49ewaZTIbOnTtj5syZAIANGzbgzp07WLp0KTp37ozOnTsj85/2dCEhIZg6dSpc\nXFzg7OyMIUOGVHqekvoXGuiSEEJe0VZTZYMtItOV0o95ly4BU6YAf/0lcFCEEGJA1q4F0tKAr756\nta7WFJHpC7UgI4SQsmp1EZm+UIIhhJCytDXpGCUYSjCEEKJEKuUtyaqLEgwlGEIIUdK6NfBP26pq\noQRDCYYQQnSCEgwlGEII0QlKMJRgCCFEJyjBUIIhhBCdoARDCYYQQnSizvbkZ4zPe5CdDZibCx0V\nIYQYNurJr4GcHMDEhJILIYToSp1NMFQ8RgghukUJhhBCiE7U2QTz/DnvrUoIIUQ36mwlPyGEEPXV\nmkr+efPmoW3btujYsSNGjhyJ7Oxspe33799HgwYN8FWpyQri4uLg7u4OFxcXzJkzR98hE0IIeY1B\nJphBgwbh+vXruHLlClxdXbFy5Uql7XPnzoWvr6/SuhkzZmDbtm1ITExEYmIioqKi9BlyjRQdHS10\nCAaDrsUrdC1eoWtRPQaZYGQyGcRiHpqXlxdSU1MV2yIiItCmTRtIpVLFurS0NOTk5MDT0xMAMHHi\nREREROg36BqI/nheoWvxCl2LV+haVI9BJpjStm/fDh8fHwDAs2fPsHr1agQFBSntI5fLYWdnp/hZ\nIpFALpfrM0xCCCGvqSfUiWUyGdLT08usDw4OxrB/ZrpZsWIFTExMMHbsWABAUFAQPv74Y9SvX58q\n6gkhxNAxAxUaGsp69uzJ8vLyFOvefPNN5ujoyBwdHZmVlRWztrZmGzduZGlpaczNzU2x3549e1hg\nYKDK4zo5OTEAtNBCCy20aLA4OTlpfB8X7AmmIlFRUVizZg1OnToFMzMzxfrTp08r/r906VI0bNgQ\nM/+Zds3S0hIxMTHw9PREWFgYZs+erfLYt2/f1m3whBBCABhoHcyHH36IZ8+eQSaToXPnzookUpGQ\nkBBMnToVLi4ucHZ2xpAhQ/QQKSGEkPLUuY6WhBBC9MMgn2C0ISoqCm5ubnBxccGqVatU7jN79my4\nuLigY8eOuHz5sp4j1J/KrsUPP/yAjh07okOHDujVqxeuXr0qQJS6p87vBAD8+eefqFevHv773//q\nMTr9UudaREdHo3Pnzmjfvj28vb31G6AeVXYtMjMzMWTIEHTq1Ant27fHjh079B+knkyZMgXNmjWD\nu7t7uftodN/UuNamBigsLGROTk4sOTmZ5efns44dO7KEhASlfQ4dOsSGDh3KGGPswoULzMvLS4hQ\ndU6da3H+/HmWlZXFGGPs8OHDtfJaqHMdSvbr168f8/X1ZQcOHBAgUt1T51o8efKESaVSlpKSwhhj\nLCMjQ4hQdU6da7FkyRK2YMECxhi/DtbW1qygoECIcHXu9OnT7NKlS6x9+/Yqt2t636yVTzCxsbFw\ndpsyDPYAAAUdSURBVHaGo6MjjI2NMWbMGERGRirt8/PPP2PSpEkAeGfOrKwsPHz4UIhwdUqda9Gj\nRw80atQIQNmOrbWFOtcBAL799lu88847sLW1FSBK/VDnWuzZswejRo1S9C+zqaVDj6tzLVq0aIGn\nT58CAJ4+fYomTZqgXj2DbB9VbW+++SYaN25c7nZN75u1MsHI5XLY29srfrazsyvT8VLVPrXxxqrO\ntSht27Ztio6ttYm6vxORkZGYMWMGAD64X22kzrVITEzE48eP0a9fP3Tt2hVhYWH6DlMv1LkW06ZN\nw/Xr19GyZUt07NgR69at03eYBkPT+2atTMPq3hjYa+0bauMNRZP3dPLkSWzfvh3nzp3TYUTCUOc6\nfPTRR/jyyy9LTatdO9u/qHMtCgoKcOnSJRw/fhy5ubno0aMHunfvDhcXFz1EqD/qXIvg4GB06tQJ\n0dHRuHPnDmQyGa5cuYKGDRvqIULDo8l9s1YmGIlEgpSUFMXPKSkpSkPJqNonNTUVEolEbzHqizrX\nAgCuXr2KadOmISoqqsJH5JpKnesQFxeHMWPGAOAVu4cPH4axsTGGDx+u11h1TZ1rYW9vDxsbG5ib\nm8Pc3Bx9+vTBlStXal2CUedanD9/HosXLwYAODk5oXXr1vjf//6Hrl276jVWQ6DxfVOrNUQGoqCg\ngLVp04YlJyezly9fVlrJ/8cff9TKim3G1LsW9+7dY05OTuyPP/4QKErdU+c6lBYQEMAOHjyoxwj1\nR51rcePGDTZgwABWWFjInj9/ztq3b8+uX78uUMS6o861+Pjjj1lQUBBjjLH09HQmkUjYo0ePhAhX\nL5KTk9Wq5Ffnvlkrn2Dq1auHDRs2YPDgwSgqKsJ7772Htm3bYvPmzQCAwMBA+Pj44LfffoOzszMs\nLCwQGhoqcNS6oc61WLZsGZ48eaKoezA2NkZsbKyQYWudOtehrlDnWri5uWHIkCHo0KEDxGIxpk2b\npjSCeW2hzrVYtGgRJk+ejI4dO6K4uBirV6+GtbW1wJHrhr+/P06dOoXMzEzY29tj6dKlKCgoAFC1\n+yZ1tCSEEKITtbIVGSGEEOFRgiGEEKITlGAIIYToBCUYQgghOkEJhhBCiE5QgiGEEKITlGAIIYTo\nBCUYQgghOkEJhhBCiE5QgiGEEKITtXIsMkJqgs2bNyM+Ph4ikQhz587Fzp07UVRUhJs3b2LLli21\nclRrUrdQgiFEAHK5HMXFxVi0aBHs7Oxgbm6O1atXAwDef/99LFy4EN99953AURJSPVRERogALly4\nAF9fX9y4cQP169dHcHCwYpulpSVOnjwpYHSEaAclGEIEMGrUKDg4OODMmTPo3bu30hzv8fHxsLS0\nFDA6QrSDEgwhAjpz5gz69eun+Dk/Px9nz55VWkdITUUJhhCBFBQUICYmBt7e3op1kZGRYIzhgw8+\nEC4wQrSEEgwhArl06RJyc3Px999/AwAePXqETz/9FJs2bUKrVq0Ejo6Q6qNWZIQI5MyZM/D29sbV\nq1dx+fJl3Lp1C5s2bcLgwYOFDo0QraAEQ4hAzpw5gx49emDx4sVCh0KITlARGSECYIzh3Llz6N69\nu9ChEKIzlGAIEcC1a9fw+PFjeHp6Ch0KITpDCYYQPTtw4ADGjBkDkUiEadOm4erVq0KHRIhOiBhj\nTOggCCGE1D70BEMIIUQnKMEQQgjRCUowhBBCdIISDCGEEJ2gBEMIIUQnKMEQQgjRCUowhBBCdIIS\nDCGEEJ34f7SxB5mAeQ9BAAAAAElFTkSuQmCC\n",
       "text": [
        "<matplotlib.figure.Figure at 0xcaaaff0>"
       ]
      }
     ],
     "prompt_number": 3
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Note that our estimator $\\hat{p}$ (red circle) is not equal to the true value of $p$ (green square), but it is at the maximum of the likelihood function. This may sound disturbing, but keep in mind this estimate is a function of the random data; and since that data can change, the ultimate estimate can likewise change. I invite you to run this notebook a few times to observe  this. Remember that the estimator is a *function* of the data and is thus also a *random variable*, just like the data is. \n",
      "\n",
      "Let's write some code to empirically examine the behavior of the maximum likelihood estimator using a simulation of multiple trials. All we're doing here is combining the last few blocks of code."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def estimator_gen(niter=10,ns=100):\n",
      "    'generate data to estimate distribution of maximum likelihood estimator'\n",
      "    out=[]\n",
      "    x=sympy.symbols('x',real=True)\n",
      "    L=   p**x*(1-p)**(1-x)\n",
      "    for i in range(niter):\n",
      "        xs = sample(ns) # generate some samples from the experiment\n",
      "        J=np.prod([L.subs(x,i) for i in xs]) # objective function to maximize\n",
      "        logJ=sympy.expand_log(sympy.log(J)) \n",
      "        sol=sympy.solve(sympy.diff(logJ,p),p)[0]\n",
      "        out.append(float(sol.evalf()))\n",
      "    return out if len(out)>1 else out[0] # return scalar if list contains only 1 term\n",
      "    \n",
      "etries = estimator_gen(100) # this may take awhile, depending on how much data you want to generate\n",
      "hist(etries) # histogram of maximum likelihood estimator\n",
      "title('$\\mu=%3.3f,\\sigma=%3.3f$'%(mean(etries),std(etries)),fontsize=18)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 4,
       "text": [
        "<matplotlib.text.Text at 0xcee5cd0>"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAXYAAAEOCAYAAACdLzzJAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAG59JREFUeJzt3X10FNXhxvFnQ4IQSCRAEoRg0CASSIhR8CU1slQCBqUC\nWitSrQna1kIrVYvVn9UQLUQ5rYocrXJOKag1VbBiEaOpZiMFedWD1KixNqkRBSGQmJeSF5jfH5Gl\nSyDszk7eLt/POXvO7t25c+/eTB6Gm5m7LsuyLAEAjBHS2R0AADiLYAcAwxDsAGAYgh0ADEOwA4Bh\nCHYAMAzBDgCGIdgBwDChnd0BmGXdunXKz89XYmKiPvzwQ02ePFk33nhjQPvYsmWLcnNztXbtWp/y\n559/Xm+++abi4+P1n//8RzNmzNDVV1/tff+dd97RLbfcorS0NA0YMEDV1dUqLS3V0qVLNWbMGEc+\nn5PsjlWg9eyOJ7oxC3DIhg0brIEDB1oHDhywLMuyamtrrbi4OGv16tV+76Ours4aMWKENWHCBJ/y\nxx57zIqPj7dqa2u928XExFgbN270blNUVGQNHjzYOv30062BAwdaP/jBD6xPPvnEgU/mPLtjFWi9\nYMYT3RdTMXDMggULNH36dPXr10+S1KdPH91www168MEH/d7H4sWLlZCQIOt/Vrqor6/X/fffr8zM\nTPXp00eSFB4erssuu0yPPvqodzuXy6VFixapqqpKe/fuVX5+vkaMGOHQp3OW3bEKtF4w44nui2CH\nIxoaGlRUVKSkpCSf8qSkJO3YsUOVlZUn3UdhYaFSUlIUGxvrU15SUqKamhrFxMT4lMfFxamwsFCH\nDx/2llndYOkju2MVaD0nxhPdE8FuuGuuuUbLli3zKZsyZYpeeuklR9spKytTc3OzIiMjfcqPvC4r\nK2uz/oEDB1RcXKxp06a1CufTTjtNUuvQbmpqUnV1tT7//HNvWWlpqe666y4tWLBAP/zhD/WXv/zF\nr/43NjZqwYIFio+PV0hIiM+jV69eqqqq8ms//rA7VoHUc2o80T3xx1OD1dXV6W9/+5vmzp3rLaus\nrNQbb7yh3NxcR9vav3+/JHn/a39E3759ve225ZFHHtG999573PeSkpIUFxenr776yqf8n//8pyRp\n3759GjZsmKSWs9GXX35ZLpdLNTU1Gj58uE477TRNmzbthG03NjYqMzNTPXv2VH5+vlwul7KysjR+\n/Hjde++96tOnj3fqwwl2xyqQek6NJ7ongt1gGzduVI8ePfSd73zHW7Z+/Xr17dtXF1xwQavtf/Sj\nH+nrr7/2a9/R0dFauXKl93VoaMuh1KNHD5/tGhsbJUnNzc0n3Nfq1as1efJkRUREeMtcLpfP8yef\nfFK33XabKisrNWDAAG3cuFFNTU0+bZ5//vlasWKFt25ERITGjx+v++67r81gv//++1VXV6c333zT\nu685c+Zo+fLlOvPMM09Yz+542R0rf+s5NZ7ovgh2g73zzjtKS0tTz549fcrS09N9ftGPWLFihe22\njszXHjs/W1NTI0knPOP96quvVFJSot/85jc+5cdOE1x11VUaNGiQFi9erKioKCUlJSktLU2bNm3S\nWWedJUk+QXZEeHi4SkpKVFVVddw+VFdXa8mSJVq1apVPoDU0NHiD7kTsjpfdsfKnnpPjie6LYDdY\ncXGxJk2a1Kps5syZjrc1ePBghYeHa8+ePT7lR6YHTnR1yrp16/TRRx8pKyvLW1ZUVKTGxkZlZWXp\ne9/7nqZPny5JGjt2rMaOHevdbsWKFRo3bpz69eunmpoajRkzRtOnT9fvf/977zbffPONXC6XwsLC\njtv++vXrdejQIU2cONGnfMOGDUpLSwtgBPxnd6z8qffqq686Mp7o3gh2Qx08eFBbt27VHXfc4S2r\nrq7Wjh079PTTT2vbtm3at2+frrjiCu/7wUzF9OzZUxkZGSopKfHZbvv27UpNTVV0dPRx9zN79mzN\nnj3bp2zChAlyuVxavny5t+z2229XUVGRPvjgA0kt0w/FxcV67LHHJEkhISFqaGjQueee67Ov0tJS\nXXTRRa3mpY/473//q4EDB/r8r2bXrl0qLCzUli1b2hwDu+Nld6z8qefUeKKb67Qr6NGuPB6P5XK5\nrBdffNGyLMs6fPiwNW/ePCskJMSqq6uzHn74Yau+vt7RNl9//XUrOjraqq6utizLsvbu3WtFRUVZ\nr7zyinebdevWWQMHDrQKCwtPuJ9LL73USk9P9ynLysqybrnlFu/ru+++25oxY4bPNvfcc49VWlrq\nfb1582ard+/e1rZt27xlb775pjVw4EBr7dq1lmVZ1p49e6yoqCirsrLSsizLamhosCZPnmwtXbo0\n0I8fELtj5U+9Y9kdT3RfnLEbqri4WOeee66ee+45ffzxx2psbNS8efNUU1Oj//u//9Po0aPVu3dv\nR9u84oor9Mgjj+jWW29VSkqKduzYoSVLlrS6Tb25ufm4fyD861//qj/84Q/atGmTXC6XJk+erJ/+\n9KeaPn26HnroIeXk5OgXv/iFqqurFR8fr/z8fJ/6OTk5Wrhwofbs2aPTTjtNX3/9td59912lpKT4\nbNfY2KiXX35ZV155pWJiYpSfn6+5c+fqnHPO0a5duzR37lxdddVVjo7NseyOlb/1pODHE92Xy7JO\nfEdHdna2XnvtNcXExGjnzp2SWtadmDt3rpqamhQaGqonn3xS48aN67AOwz+XX365JkyYoPvuu6+z\nu9IlPfDAA1qwYEFndwNoF23eoJSVlaWCggKfsvnz5+vBBx/U+++/r9zcXM2fP79dO4jANTU1adOm\nTT6XOeKogwcPckkfjNZmsKenpysqKsqn7IwzzlB1dbUkqaqqSkOGDGm/3sGWrVu3qqmpSRdffHFn\nd6VLeu655/T973+/s7sBtJuA59jz8vJ06aWX6q677tLhw4f17rvvtke/EISysjJlZmY6Podugqqq\nKn399ddKTEzs7K4A7abNOXZJKi8v19SpU71z7BMnTtScOXM0ffp0vfTSS3rmmWdUWFjYIZ0FAPjh\nZJfNlJWVWUlJSd7XERER3ueHDx+2IiMjj1svISHBksSDBw8ePAJ4JCQkBHmxo4312IcPH67i4mJJ\n0ttvv33Cu+Q+++wzWZbFw7L0wAMPdHofusqDsWAsGIu2H5999lmgsdxKm3PsM2fOVHFxsfbt26eh\nQ4cqNzdXzzzzjObMmaOGhgb17t1bzzzzTNCdAAA4p81gf+GFF45bvnnz5nbpDAAgeHzRRgdwu92d\n3YUug7E4irE4irFw1kmvirG9Y5dL7bRrADCWE9nJGTsAGIZgBwDDEOwAYBiCHQAMQ7ADgGEIdgAw\nDMEOOCAysr9cLle7PyIj+3f2R0U3wHXsgANcLpda1nBq95b4vTIc17EDAFoh2AHAMAQ7ABiGYAcA\nwxDsAGAYgh0ADNNmsGdnZys2NlbJyck+5U888YQSExOVlJSku+++u107CAAITJvfoJSVlaWf//zn\nuummm7xlRUVFevXVV/XBBx8oLCxMe/fubfdOAgD81+YZe3p6uqKionzKnnrqKd1zzz0KCwuTJEVH\nR7df7wAAAQt4jv3TTz/VO++8o4svvlhut1vbtm1rj34BAGxqcyrmeJqbm3XgwAFt2rRJW7du1XXX\nXad///vfx902JyfH+9ztdvO9hgBwDI/HI4/H4+g+T7pWTHl5uaZOnaqdO3dKkjIzM/XrX/9a48eP\nlyQNHz5cmzdv1oABA3x3zFoxOIWwVgyc0ilrxUybNk1vv/22JKm0tFSNjY2tQh0A0HnanIqZOXOm\niouLVVlZqaFDhyo3N1fZ2dnKzs5WcnKyevbsqZUrV3ZUXwEAfmDZXsABTMXAKSzbCwBohWAHAMMQ\n7ABgGIIdAAxDsAOAYQh2ADAMwQ4AhiHYAcAwBDuA44qM7C+Xy9Uhj8jI/p39cY3CnaeAA0y887Tj\nPpPEHbVHcecpAKAVgh0ADEOwA4BhCHYAMAzBDgCGIdgBwDBtBnt2drZiY2OVnJzc6r3f/e53CgkJ\n0f79+9utcwCAwLUZ7FlZWSooKGhVXlFRocLCQsXHx7dbxwAA9rQZ7Onp6YqKimpVfscdd+iRRx5p\nt04BAOwLeI59zZo1iouL05gxY9qjPwCAIIUGsnF9fb0WLlyowsJCb1lbt77m5OR4n7vdbrnd7oA7\nCAAm83g88ng8ju7zpGvFlJeXa+rUqdq5c6d27typiRMnKjw8XJL0xRdfaMiQIdqyZYtiYmJ8d8xa\nMTiFsFZM0K2RF99yIjsDOmNPTk7Wnj17vK/POussbd++Xf37szIbAHQVbc6xz5w5U2lpaSotLdXQ\noUO1fPlyn/db/kUHAHQlLNsLOICpmKBbIy++xbK9AIBWCHYAMAzBDgCGIdgBwDAEOwAYhmAHAMMQ\n7ABgGIIdAAwT0JICQHcSGdlfNTUHOrsbQIfjzlMYq6PvnOTO06BaIy++xZ2nAIBWCHYAMAzBDgCG\nIdgBwDAEOwAYhmAHAMOcNNizs7MVGxur5ORkb9mvfvUrJSYmKiUlRTNmzFB1dXW7dhIA4L+TBntW\nVpYKCgp8yiZNmqQPP/xQO3bs0IgRI7Ro0aJ26yAAIDAnDfb09HRFRUX5lGVkZCgkpKXqRRddpC++\n+KJ9egcACFjQc+x//OMfNWXKFCf6AgBwQFBrxfz2t79Vz549dcMNNxz3/ZycHO9zt9stt9sdTHMA\nYByPxyOPx+PoPv1aK6a8vFxTp07Vzp07vWV/+tOftGzZMr311lvq1atX6x2zVgw6GWvFBNkSa8V0\nCiey09YZe0FBgRYvXqzi4uLjhjoAoPOc9Ix95syZKi4u1r59+xQbG6sFCxZo0aJFamxsVP/+/SVJ\nl1xyiZ588knfHXPGjk7GGXuQLXHG3imcyE6W7YWxCPYgWyLYOwXL9gIAWiHYAcAwBDsAGIZgBwDD\nEOwAYBiCHQAMQ7ADgGGCWisGQEcL/fb6cuDECHagW2lWx950he6IqRgAMAzBDgCGIdgBwDAEOwAY\nhmAHAMMQ7ABgmDaDPTs7W7GxsUpOTvaW7d+/XxkZGRoxYoQmTZqkqqqqdu8kAMB/bQZ7VlaWCgoK\nfMry8vKUkZGh0tJSXX755crLy2vXDgIAAnPSb1A69ousR44cqeLiYsXGxmr37t1yu936+OOPW++Y\nb1BCJzP1G5TM+0wtbZEXLTrlG5T27Nmj2NhYSVJsbKz27NkTVAcAAM4K6o+nLpeLdSsAoIsJeK2Y\nI1MwgwYN0ldffaWYmJgTbpuTk+N97na75Xa77fQR7Swysr9qag50UGthkpo6qC2g6/N4PPJ4PI7u\nM+A59vnz52vAgAG6++67lZeXp6qqquP+AZU59u7DzLloU9sy8TO1tEVetHAiO9sM9pkzZ6q4uFj7\n9u1TbGyscnNzdfXVV+u6667T559/rmHDhunFF19Uv3792qVz6BgEe3dqy8TP1NIWedGi3YM9qB0T\n7N0Gwd6d2jLxM7W0RV606JSrYgAAXRvBDgCGIdgBwDAEOwAYhmAHAMMQ7ABgGIIdAAxDsAOAYQh2\nAKeUyMj+3gUM2/sRGdm/Uz4jd56CO0+7VVsmfqaWtjoqLzr6eA/0c3HnKQCgFYIdAAxDsAOAYQh2\nADAMwQ4AhiHYAcAwtoN90aJFGj16tJKTk3XDDTeooaHByX4BAGyyFezl5eVatmyZ3nvvPe3cuVOH\nDh1Sfn6+030DANgQaqdSZGSkwsLCVF9frx49eqi+vl5Dhgxxum8AABtsnbH3799fd955p84880wN\nHjxY/fr108SJE53uGwDABlvB/tlnn+mxxx5TeXm5vvzyS9XW1ur55593um8AABtsTcVs27ZNaWlp\nGjBggCRpxowZ2rhxo2bNmuWzXU5Ojve52+2W2+223VEAMJHH45HH43F0n7YWAduxY4dmzZqlrVu3\nqlevXrr55pt14YUXas6cOUd3zCJg3QaLgHWntkz8TC1tsQjYtzU6axGwlJQU3XTTTRo7dqzGjBkj\nSfrxj38cVEcAAM5g2V5wxt6t2jLxM7W0xRn7tzVYthcAcCyCHQAMQ7ADgGEIdgAwDMEOAIYh2AHA\nMAQ7ABjG1pICp6rXXntNW7du7ZC2+vbtqzvvvPPba24BwH/coBSAceMytG3bIEkJ7d5WSMhCVVVV\nKiIiot3b4gal7tSWiZ9JksIkNXdQW5LpNyhxxh6wmyRltHsrPXr8rt3bALqOZnXsP1hmY44dAAxD\nsAOAYQh2ADAMwQ4AhiHYAcAwBDsAGMZ2sFdVVenaa69VYmKiRo0apU2bNjnZLwCATbavY7/99ts1\nZcoUrVq1Ss3Nzaqrq3OyXwAAm2wFe3V1tdavX68VK1a07CQ0VKeffrqjHQMA2GNrKqasrEzR0dHK\nysrS+eefr1tvvVX19fVO9w0AYIOtM/bm5ma99957Wrp0qcaNG6d58+YpLy9Pubm5Ptvl5OR4n7vd\nbrnd7mD6CgDG8Xg88ng8ju7T1iJgu3fv1iWXXKKysjJJ0j/+8Q/l5eVp7dq1R3ds7CJg89URa8WE\nhUWosvJLFgGjrU5qh7acaqszFgGzNRUzaNAgDR06VKWlpZKkv//97xo9enRQHQEAOMP2VTFPPPGE\nZs2apcbGRiUkJGj58uVO9gsAYJPtYE9JSemwL50AAPiPO08BwDAEOwAYhmAHAMMQ7ABgGIIdAAxD\nsAOAYQh2ADAMwQ4AhiHYAcAwBDsAGIZgBwDDEOwAYBiCHQAMQ7ADgGEIdgAwDMEOAIYJKtgPHTqk\n1NRUTZ061an+AACCFFSwP/744xo1atS3X4YMAOgKbAf7F198oXXr1umWW24J+hu1AQDOsR3sv/zl\nL7V48WKFhDBNDwBdia0vs167dq1iYmKUmpoqj8dzwu1ycnK8z91ut9xut53mAMBYHo+nzRy1w2XZ\nmEe599579eyzzyo0NFQHDx7UN998o2uuuUYrV648umOXy7gpmnHjMrRt23xJGe3eVlhYhCorv1RE\nRES7t9XyN5KO+lnRVvdoh7acaivQHHQiO23NoyxcuFAVFRUqKytTfn6+vvvd7/qEOgCg8zgyQc5V\nMQDQddiaY/9f48eP1/jx453oCwDAAVzSAgCGIdgBwDAEOwAYhmAHAMMQ7ABgGIIdAAxDsAOAYQh2\nADAMwQ4AhiHYAcAwBDsAGIZgBwDDEOwAYBiCHQAMQ7ADgGEIdgAwjO1gr6io0IQJEzR69GglJSVp\nyZIlTvYLAGCT7W9QCgsL06OPPqrzzjtPtbW1uuCCC5SRkaHExEQn+wcACJDtM/ZBgwbpvPPOkyT1\n7dtXiYmJ+vLLLx3rGADAHkfm2MvLy/X+++/roosucmJ3AIAgBP1l1rW1tbr22mv1+OOPq2/fvj7v\n5eTkeJ+73W653e5gmwMAo3g8Hnk8Hkf36bIsy7JbuampSVdddZUyMzM1b9483x27XApi113SuHEZ\n2rZtvqSMdm8rLCxClZVfKiIiot3bcrlckjrqZ0Vb3aMd2nKqrUBz0InstD0VY1mWZs+erVGjRrUK\ndQBA57Ed7Bs2bNBzzz2noqIipaamKjU1VQUFBU72DQBgg+059ksvvVSHDx92si8AAAdw5ykAGIZg\nBwDDEOwAYBiCHQAMQ7ADgGEIdgAwDMEOAIYJeq0YtJ/Bg+NVW3ugs7sBoJsh2LuwllDvqPVHAJiC\nqRgAMAzBDgCGIdgBwDAEOwAYhmAHAMMQ7ABgGNvBXlBQoJEjR+qcc87Rww8/7GSfAABBsBXshw4d\n0ty5c1VQUKCSkhK98MIL+uijj5zum0E8nd0BdEmezu4ADGUr2Lds2aLhw4dr2LBhCgsL0/XXX681\na9Y43TeDeDq7A+iSPJ3dARjKVrDv2rVLQ4cO9b6Oi4vTrl27HOsUAMA+W0sKuFyn5i3oYWEhCg//\njUJDlwRU7+DBT9Sr1/aA6tTW/jeg7QHgCFvBPmTIEFVUVHhfV1RUKC4uzmebhISEU/YfgONpbPzU\nZs2OGsOO/FnR1lELOqgdu2gr6JYCzMGEhITg27QsK+BVppqbm3Xuuefqrbfe0uDBg3XhhRfqhRde\nUGJiYtAdAgAEx9YZe2hoqJYuXarJkyfr0KFDmj17NqEOAF2ErTN2AEDXZeuqmJPdnLRmzRqlpKQo\nNTVVF1xwgd5++23ve8OGDdOYMWOUmpqqCy+80H7Puwh/b9TaunWrQkNDtXr16oDrdhfBjMWpdlx4\nPB6dfvrpSk1NVWpqqh566CG/63Y3gY7Fgw8+6H3PpOPCn5+rx+NRamqqkpKS5Ha7A6rrwwpQc3Oz\nlZCQYJWVlVmNjY1WSkqKVVJS4rNNbW2t9/kHH3xgJSQkeF8PGzbMqqysDLTZLsmfsTiy3YQJE6wr\nr7zSWrVqVUB1u4tgxsKyTr3joqioyJo6daqtut1JMGNhWeYcF/6Mw4EDB6xRo0ZZFRUVlmVZ1t69\ne/2ue6yAz9j9uTmpT58+3ue1tbUaOHDgsf+YBNpsl+TvjVpPPPGErr32WkVHRwdct7sIZiyOONWO\ni+N93lP1uGjrZ2/CceHPOPz5z3/WNddc473C8Ehu2jkmAg52f29OeuWVV5SYmKjMzEwtWXL0um+X\ny6WJEydq7NixWrZsWaDNdyn+jMWuXbu0Zs0a3XbbbZKOXvpk2k1ewYzFkeen0nHhcrm0ceNGpaSk\naMqUKSopKfG7bncSzFgcec+E48Kfcfj000+1f/9+TZgwQWPHjtWzzz7rd91jBXxVjL/XZE6bNk3T\npk3T+vXrdeONN+qTTz6RJG3YsEFnnHGG9u7dq4yMDI0cOVLp6emBdqNL8Gcs5s2bp7y8PLlcLlmW\n5T37MO0a/2DGQjr1jovzzz9fFRUVCg8P1+uvv65p06aptLS0A3rXsYIdC1OOC3/GoampSe+9957e\neust1dfX65JLLtHFF19sKysCPmP35+ak/5Wenq7m5mZVVlZKks444wxJUnR0tKZPn64tW7YE2oUu\nw5+x2L59u66//nqdddZZWr16tX72s5/p1VdfDXgcu7pgxkI69Y6LiIgIhYeHS5IyMzPV1NSk/fv3\nKy4u7pQ7Lk40FpI5x4U/4zB06FBNmjRJvXv31oABA3TZZZdpx44d9rIi0D8CNDU1WWeffbZVVlZm\nNTQ0HHci/1//+pd1+PBhy7Isa/v27dbZZ59tWZZl1dXVWd98841lWS1/YE1LS7PeeOONQLvQZfgz\nFv/r5ptvtlavXm2rblcXzFicisfF7t27vb8jmzdvtuLj4/2u250EMxYmHRf+jMNHH31kXX755VZz\nc7NVV1dnJSUlWR9++KGtYyLgqZgT3Zz09NNPS5J+8pOfaPXq1Vq5cqXCwsLUt29f5efnS5J2796t\nGTNmSGq5e3XWrFmaNGlSoF3oMvwZi0DrdlfBjMWpeFysWrVKTz31lEJDQxUeHu79HTkVj4sTjYVJ\nx4U/4zBy5EhdccUVGjNmjEJCQnTrrbdq1KhRkhTwMcENSgBgGL4aDwAMQ7ADgGEIdgAwDMEOAIYh\n2AHAMAQ7ABiGYAcAwxDsAGCY/wec6io5ZiGzMQAAAABJRU5ErkJggg==\n",
       "text": [
        "<matplotlib.figure.Figure at 0xcc03290>"
       ]
      }
     ],
     "prompt_number": 4
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Note that the mean of the estimator ($\\mu$) is pretty close to the true value, but looks can be deceiving. The only way to know for sure is to check if the estimator is unbiased, namely, if\n",
      "\n",
      "$$ \\mathbb{E}(\\hat{p}) = p $$"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Because this problem is simple, we can solve for this in general noting that since $x=0$ or $x=1$, the terms in the product of $\\mathcal{L}$ above are either $p$, if $x_i=1$ or $1-p$ if $x_i=0$. This means that we can write\n",
      "\n",
      "$$ \\mathcal{L}(p|\\mathbf{x})= p^{\\sum_{i=1}^n x_i}(1-p)^{n-\\sum_{i=1}^n x_i} $$\n",
      "\n",
      "with corresponding log as\n",
      "\n",
      "$$ J=\\log(\\mathcal{L}(p|\\mathbf{x})) =  \\log(p)  \\sum_{i=1}^n x_i +   \\log(1-p) \\left(n-\\sum_{i=1}^n x_i\\right)$$ \n",
      "\n",
      "Taking the derivative of this gives:\n",
      "\n",
      "$$  \\frac{dJ}{dp} = \\frac{1}{p}\\sum_{i=1}^n x_i + \\frac{(n-\\sum_{i=1}^n x_i)}{p-1} $$\n",
      "\n",
      "and solving this leads to\n",
      "\n",
      "$$  \\hat{p} = \\frac{1}{ n} \\sum_{i=1}^n x_i $$\n",
      "\n",
      "This is our *estimator* for $p$. Up til now, we have been using `sympy` to solve for this based on the data $x_i$ but now we have it generally and don't have to solve for it again. To check if this estimator is biased, we compute its expectation:\n",
      "\n",
      "$$ \\mathbb{E}\\left(\\hat{p}\\right) =\\frac{1}{n}\\sum_i^n \\mathbb{E}(x_i) = \\frac{1}{n} n \\mathbb{E}(x_i) $$\n",
      "\n",
      "by linearity of the expectation and where\n",
      "\n",
      "$$\\mathbb{E}(x_i)  = p$$\n",
      "\n",
      "Therefore,\n",
      "\n",
      "$$ \\mathbb{E}\\left(\\hat{p}\\right) =p $$\n",
      "\n",
      "This means that the esimator is unbiased. This is good news. We almost always want our estimators to be unbiased. Similarly, \n",
      "\n",
      "$$ \\mathbb{E}\\left(\\hat{p}^2\\right) = \\frac{1}{n^2} \\mathbb{E}\\left[\\left(  \\sum_{i=1}^n x_i \\right)^2 \\right]$$\n",
      "\n",
      "and where\n",
      "\n",
      "$$ \\mathbb{E}\\left(x_i^2\\right) =p$$\n",
      "\n",
      "and by the independence assumption,\n",
      "\n",
      "$$ \\mathbb{E}\\left(x_i x_j\\right) =\\mathbb{E}(x_i)\\mathbb{E}( x_j) =p^2$$\n",
      "\n",
      "Thus,\n",
      "\n",
      "$$ \\mathbb{E}\\left(\\hat{p}^2\\right) =\\left(\\frac{1}{n^2}\\right) n \n",
      "\\left[\n",
      "p+(n-1)p^2\n",
      "\\right]\n",
      "$$\n",
      "\n",
      "So, the variance of the estimator, $\\hat{p}$ is the following:\n",
      "\n",
      "$$ \\sigma_\\hat{p}^2 = \\mathbb{E}\\left(\\hat{p}^2\\right)- \\mathbb{E}\\left(\\hat{p}\\right)^2  = \\frac{p(1-p)}{n} $$\n",
      "\n",
      "Note that the $n$ in the denominator means that the variance asymptotically goes to zero as $n$ increases (i.e. we consider more and more samples). This is good news also because it means that more and more coin flips leads to a better estimate of the underlying $p$.\n",
      "\n",
      "Unfortunately, this formula for the variance is practically useless because we have to know $p$ to compute it and $p$ is the parameter we are trying to estimate in the first place! But, looking at $ \\sigma_\\hat{p}^2 $, we can immediately notice that if $p=0$, then there is no estimator variance because the outcomes are guaranteed to be tails. Also, the maximum of this variance, for whatever $n$, happens at $p=1/2$. This is our worst case scenario and the only way to compensate is with more samples (i.e. larger $n$). \n"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "All we have computed is the mean and variance of the estimator. In general, this is insufficient to characterize the underlying probability density of $\\hat{p}$, except if we somehow knew that  $\\hat{p}$ were normally distributed. This is where the powerful [*central limit theorem*](http://mathworld.wolfram.com/CentralLimitTheorem.html) comes in. The form of the estimator, which is just a mean estimator, implies that we can apply this theorem and conclude that  $\\hat{p}$ is normally distributed. However, there's a wrinkle here: the theorem tells us that  $\\hat{p}$ is asymptotically normal, it doesn't quantify how many samples $n$ we need to approach this asymptotic paradise. In our simulation this is no problem since we can generate as much data as we like, but in the real world, with a costly experiment, each sample may be precious. In the following, we won't apply this theorem and instead proceed analytically.\n"
     ]
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "Probability Density for the Estimator"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "To write out the full density for $\\hat{p}$, we first have to ask what is the probability that the estimator will equal a specific value and the tally up all the ways that could happen with their corresponding probabilities. For example, what is the probability that\n",
      "\n",
      "$$  \\hat{p} = \\frac{1}{n}\\sum_{i=1}^n x_i  = 0 $$\n",
      "\n",
      "This can only happen one way: when $x_i=0 \\hspace{0.5em} \\forall i$. The probability of this happening can be computed from the density\n",
      "\n",
      "$$ f(\\mathbf{x},p)= \\prod_{i=1}^n \\left(p^{x_i} (1-p)^{1-x_i}  \\right) $$\n",
      "\n",
      "$$ f\\left(\\sum_{i=1}^n x_i  = 0,p\\right)= \\left(1-p\\right)^n  $$\n",
      "\n",
      "Likewise, if $\\lbrace x_i \\rbrace$ has one $i^{th}$  value equal to one, then\n",
      "\n",
      "$$ f\\left(\\sum_{i=1}^n x_i  = 1,p\\right)= n p \\prod_{i=1}^{n-1} \\left(1-p\\right)$$\n",
      "\n",
      "where the $n$ comes from the $n$ ways to pick one value equal to one from the $n$ elements $x_i$. Continuing this way, we can construct the entire density as\n",
      "\n",
      "$$ f\\left(\\sum_{i=1}^n x_i  = k,p\\right)= \\binom{n}{k} p^k  (1-p)^{n-k}  $$\n",
      "\n",
      "where the term on the left is the binomial coefficient of $n$ things taken $k$ at a time. This is the binomial distribution and it's not the density for $\\hat{p}$, but rather for $n\\hat{p}$. We'll leave this as-is because it's easier to work with below. We just have to remember to keep track of the $n$ factor."
     ]
    },
    {
     "cell_type": "heading",
     "level": 4,
     "metadata": {},
     "source": [
      "Confidence Intervals"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Now that we have the full density for $\\hat{p}$, we are ready to ask some meaningful questions. For example,\n",
      "\n",
      "$$ \\mathbb{P}\\left( | \\hat{p}-p | \\le \\epsilon p \\right) $$\n",
      "\n",
      "Or, in words, what is the probability we can get within $\\epsilon$ percent of the true value of $p$. Rewriting,\n",
      "\n",
      "$$ \\mathbb{P}\\left(  p - \\epsilon p \\lt \\hat{p} \\lt p + \\epsilon p \\right) = \\mathbb{P}\\left(  n p - n \\epsilon p \\lt \\sum_{i=1}^n x_i   \\lt n p + n \\epsilon p \\right)$$\n",
      "\n",
      "Let's plug in some live numbers here for our worst case scenario where $p=1/2$. Then, if $\\epsilon = 1/100$, we have\n",
      "\n",
      "$$ \\mathbb{P}\\left( \\frac{99 n}{200} \\lt \\sum_{i=1}^n x_i   \\lt \\frac{101 n}{200} \\right)$$\n",
      "\n",
      "Since the sum in integer-valued, we need $n> 100$ to even compute this. Thus, if $n=101$ we have\n",
      "\n",
      "$$ \\mathbb{P}\\left( \\frac{9999}{200} \\lt \\sum_{i=1}^{101} x_i   \\lt \\frac{10201}{200} \\right) = f\\left(\\sum_{i=1}^{101} x_i  = 50,p\\right)= \\binom{101}{50} (1/2)^{50}  (1-1/2)^{101-50} = 0.079$$\n",
      "\n",
      "This means that in the worst-case scenario for $p=1/2$, given $n=101$ trials, we will only get within 1% of the actual $p=1/2$ about 8% of the time. If you feel disappointed, that only means you've been paying attention. What if the coin was really heavy and it was costly to repeat this 101 times? Then, we would be within 1% of the actual value only 8% of the time. Those odds are terrible.\n",
      "\n",
      "Let's come at this another way: given I could only flip the coin 100 times, how close could I come to the true underlying value with high probability (say, 95%)? In this case we are seeking to solve for $\\epsilon$. Plugging in gives,\n",
      "\n",
      "$$ \\mathbb{P}\\left(  50 - 50 \\epsilon \\lt \\sum_{i=1}^{100} x_i   \\lt 50 + 50 \\epsilon  \\right) = 0.95$$\n",
      "\n",
      "which we have to solve for $\\epsilon$. Fortunately, all the tools we need to solve for this are already in `scipy`."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import scipy.stats\n",
      "\n",
      "b=scipy.stats.binom(100,.5) # n=100, p = 0.5, distribution of the estimator \\hat{p}\n",
      "\n",
      "f,ax= subplots()\n",
      "ax.stem(arange(0,101),b.pmf(arange(0,101))) # heres the density of the sum of x_i\n",
      "\n",
      "g = lambda i:b.pmf(arange(-i,i)+50).sum() # symmetric sum the probability around the mean\n",
      "print 'this is pretty close to 0.95:%r'%g(10)\n",
      "ax.vlines( [50+10,50-10],0 ,ax.get_ylim()[1] ,color='r',lw=3.)\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "this is pretty close to 0.95:0.95395593307062954\n"
       ]
      },
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 5,
       "text": [
        "<matplotlib.collections.LineCollection at 0xd079630>"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEACAYAAAC08h1NAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XtQG9fdN/CvYnhxQxo78Ti4lsgQS7IlfJFxcTRPJrGV\nOmYD2DQXtyVJEyYmKY9TLDuTPzr1ZYzb+kIzz9tAaCbYfcsbJ30x08sMFCkyoYlMLpZJalInhdaY\ngY4kx0ybRGkdx2CU8/5BkBFIKwFCEuz3M6PxXs6ufnuQflrv2T1HJYQQICIixbgu0QEQEVF8MfET\nESkMEz8RkcIw8RMRKQwTPxGRwjDxExEpTMTE73A4YDAYoNfrUVlZGbKM1WqFXq+HyWRCR0dHYHlV\nVRVWrlyJFStWoKqqKnZRExHRpMkmfr/fj/LycjgcDnR2dqK+vh5dXV1BZex2O86fP4/u7m4cOXIE\n27ZtAwB8+OGH+NWvfoV3330Xf/nLX9Dc3Iyenp7pOxIiIoqKbOJvb2+HTqdDVlYWUlNTUVxcjMbG\nxqAyTU1NKCkpAQCYzWb4fD5cvHgRXV1dMJvNmDt3LubMmYP169fjD3/4w/QdCRERRUU28Xu9XmRm\nZgbmNRoNvF5vxDIXLlzAypUr8eabb+KTTz7B5cuXYbPZ4PF4Yhw+ERFNVIrcSpVKFdVOQvX6YDAY\n8KMf/Qh5eXlIT09HTk4OrruObclERIkmm/jVajXcbndg3u12Q6PRyJbxeDxQq9UAgK1bt2Lr1q0A\ngF27duHWW28d9x46nY7X/omIJkir1eL8+fOT21jIuHr1qliyZIno7e0VAwMDwmQyic7OzqAyNptN\n5OfnCyGEOHXqlDCbzYF1/f39Qggh/vGPfwiDwSA+++yzce8RIQRF2bdvX6JDSBqsi2tYF9ewLq6Z\nSu6UPeNPSUlBTU0NJEmC3+9HaWkpjEYjamtrAQBlZWUoKCiA3W6HTqdDeno66urqAttv2bIFH3/8\nMVJTU/HCCy/gxhtvnNyvExERxYxs4geA/Px85OfnBy0rKysLmq+pqQm5bVtb2xRCIyKi6cDW1iRi\nsVgSHULSYF1cw7q4hnURG6qvrhUlLgCVKuRdQUREFN5UcifP+ImIFIaJn4hIYZj4iYgUhomfiEhh\nmPiJiBSGiZ+ISGGY+ImIFIaJn4hIYZj4iYgUhomfiEhhmPiJiBSGiZ+ISGGY+ImIFIaJn4hIYSIm\nfofDAYPBAL1ej8rKypBlrFYr9Ho9TCYTOjo6AssPHTqE5cuXY+XKlXj44YcxMDAQu8iJiGhSZBO/\n3+9HeXk5HA4HOjs7UV9fj66urqAydrsd58+fR3d3N44cOYJt27YBAPr6+nD06FGcOXMGH3zwAfx+\nP44fPz59R0JERFGRHXqxvb0dOp0OWVlZAIDi4mI0NjbCaDQGyjQ1NaGkpAQAYDab4fP50N/fjxtv\nvBGpqam4fPky5syZg8uXL0OtVk/fkRDFkM3WhurqFgwMpCAtbQhWax4KC9clOiyimJBN/F6vF5mZ\nmYF5jUaD06dPRyzj9XqxZs0aPPPMM7j11lvxta99DZIk4Z577olx+ESxZ7O1YceOE+jpORBY1tOz\nGwCY/GlWkE38KpUqqp2EGv6rp6cHzz33HPr6+jBv3jx85zvfwW9+8xs88sgj48pWVFQEpi0WC8fV\nVLpRn4eg6Tiprm4JSvoA0NNzAM8/vzf+iT/BdUHJw+l0wul0xmRfsolfrVbD7XYH5t1uNzQajWwZ\nj8cDtVoNp9OJO+64AwsWLAAAPPDAA3jnnXciJn4i7N9/bToBn42BgdBfiytX5sQ5EiS8Lih5jD0p\n3j/6szFBso27ubm56O7uRl9fHwYHB9HQ0ICioqKgMkVFRTh27BgAwOVyYf78+cjIyMCyZcvgcrnw\nxRdfQAiB1tZWZGdnTzpQoulks7VBkvbAYqnAhx92jVnbBmAPzp7thSTtgc3WlogQiWJG9ow/JSUF\nNTU1kCQJfr8fpaWlMBqNqK2tBQCUlZWhoKAAdrsdOp0O6enpqKurAwCsXr0ajz32GHJzc3Hddddh\nzZo1+MEPfjD9R0Q0QeOv6bchJeW/MTT04lfzJwAcwKefAi0tvN5PM59KhLpAH88AVKqQbQSkYKPb\nluLw2ZCkPWhp+dmYpW1YsOAFfPxx6FuQJWkvHI6fTnts8a4Lmjmmkjv55C4pXuhr+uuwYoUh7DYJ\nud5PFCNM/KR4aWlDIZfPnesPu43cOqJkx8RPime15kGr3R20TKvdhe3bN341HX4d0UzExE+KV1i4\nDlVVEiRpL4Dh6/dVVfcGGm/l1hHNRGzcpeSTwAZNlSr4LUfPj10Xt4BG8HtCo7Bxl4iIosbET0Sk\nMEz8REQKI/vkLtFsFqrrZUC+0ZbdNdNswMRPihS+62X5JM7ummk24F09lHzicCdL6G4axr/l2Lt6\nQu9rGrtv4F09FAbv6iGaoHBdL08Gu2+gmYaJnxQpXDcNk8HuG2imYeInRQrXTUMk7L6BZgMmflKk\ncN00RMLuG2g2YOMuJZ84N2jKdcsQbl3cum9g4y6FMa2Nuw6HAwaDAXq9HpWVlSHLWK1W6PV6mEwm\ndHR0AAD+/ve/IycnJ/CaN28eqqurJxUkERHFjuwZv9/vx7Jly9Da2gq1Wo21a9eivr4eRqMxUMZu\nt6OmpgZ2ux2nT5/Gjh074HK5gvbz5ZdfQq1Wo729HZmZmcEB8IyfxuIZf3AAI/g9oVGm7Yy/vb0d\nOp0OWVlZSE1NRXFxMRobG4PKNDU1oaSkBABgNpvh8/nQ398fVKa1tRVarXZc0icioviTTfxerzco\nWWs0Gni93ohlPB5PUJnjx4/j4YcfjkW8REQ0RbJPsajCPao4xtj/bozebnBwEH/84x/Dtg8AQEVF\nRWDaYrHAYrFE9b5ERErhdDrhdDpjsi/ZxK9Wq+F2uwPzbrcbGo1GtozH44FarQ7Mv/rqq/jmN7+J\nhQsXhn2f0YmfiIjGG3tSvH///knvS/ZST25uLrq7u9HX14fBwUE0NDSgqKgoqExRURGOHTsGAHC5\nXJg/fz4yMjIC6+vr6/HQQw9NOkAiIoot2TP+lJQU1NTUQJIk+P1+lJaWwmg0ora2FgBQVlaGgoIC\n2O126HQ6pKeno66uLrD9559/jtbWVhw9enR6j4KIiKLGB7go+fB2zuAARvB7QqOwd04iIooaB2Ih\nRZnMqFsT3R/77qFkx8RPijHZUbcmvj+OyEXJjdf4KflM03XtaEbdmsg1/ry80PuL6YhcvMZPYfAa\nP1EUYjnqltz+OCIXJTsmflKMWI66Jbc/jshFyY6JnxRjsqNuTXR/HJGLkh0TPynGZEfdmuj+2LBL\nyY6Nu5R84tCgGe2DWZMpF/NAR/B7QqOwcZeIiKLGxE9EpDBM/ERECsPET0SkMEz8REQKw8RPRKQw\nTPxERAoTMfE7HA4YDAbo9fqwA6ZbrVbo9XqYTCZ0dHQElvt8PmzZsgVGoxHZ2dlwuVyxi5yIiCZF\nNvH7/X6Ul5fD4XCgs7MT9fX16OrqCipjt9tx/vx5dHd348iRI9i2bVtg3Y4dO1BQUICuri6cPXsW\nRqNxeo6CiIiiJpv429vbodPpkJWVhdTUVBQXF6OxsTGoTFNTE0pKSgAAZrMZPp8P/f39+Oyzz/Dm\nm29i69atAIbH7503b940HQYREUVLNvF7vV5kZmYG5jUaDbxeb8QyHo8Hvb29WLhwIR5//HGsWbMG\nTz75JC5fvhzj8ImIaKJkOyhXje4nRMbY/iJUKhWGhoZw5swZ1NTUYO3atdi5cycOHz6Mn/zkJ+O2\nr6ioCExbLBZYLJao3peISCmcTiecTmdM9iWb+NVqNdxud2De7XZDo9HIlvF4PFCr1RBCQKPRYO3a\ntQCALVu24PDhwyHfZ3TiJyKi8caeFO/fv3/S+5K91JObm4vu7m709fVhcHAQDQ0NKCoqCipTVFSE\nY8eOAQBcLhfmz5+PjIwMLFq0CJmZmTh37hwAoLW1FcuXL590oESTYbO1QZL2wGKpgCTtgc3WNive\ni2gqZM/4U1JSUFNTA0mS4Pf7UVpaCqPRiNraWgBAWVkZCgoKYLfbodPpkJ6ejrq6usD2zz//PB55\n5BEMDg5Cq9UGrSOabnKDoU92gHU5HHidZgr2x0/JJ0Z90IcbXF2S9uLEiZ/GvD/+0DFMceB19sdP\nYbA/fqIQkmEwdA68TsmIiZ9mrWQYDJ0Dr1MyYuKnWSveg6Fz4HWaKZj4adaK92DoHHidZgo27lLy\nmYYGzVgMoj6ZcjEJfAS/JzQKG3eJiChqTPxERArDxE9EpDBM/ERECsPET0SkMEz8REQKw8RPRKQw\nTPxERArDxE9EpDBM/EREChMx8TscDhgMBuj1elRWVoYsY7VaodfrYTKZ0NHREVielZWFVatWIScn\nB7fffnvsoiYiokmTHYHL7/ejvLwcra2tUKvVWLt2LYqKimA0GgNl7HY7zp8/j+7ubpw+fRrbtm2D\ny+UCMNyXhNPpxM033zy9R0FERFGTPeNvb2+HTqdDVlYWUlNTUVxcjMbGxqAyTU1NKCkpAQCYzWb4\nfD709/cH1rMDNiKi5CKb+L1eLzIzMwPzGo0GXq836jIqlQr33HMPcnNzcfTo0VjGTUREkyR7qUcV\nbiDRMcKd1b/11ltYvHgx/vnPf2Ljxo0wGAy46667Jh4lERHFjGziV6vVcLvdgXm32w2NRiNbxuPx\nQK1WAwAWL14MAFi4cCHuv/9+tLe3h0z8FRUVgWmLxQKLxTLhAyEims2cTiecTmdM9iU7EMvQ0BCW\nLVuGP/3pT1i8eDFuv/121NfXj2vcrampgd1uh8vlws6dO+FyuXD58mX4/X58/etfx+eff468vDzs\n27cPeXl5wQFwIBYaa4qDj9hsbaiubsHAQArS0oZgteZh06Z1cR2Ipbl5fAyTGo2LA7FQGFPJnbJn\n/CkpKaipqYEkSfD7/SgtLYXRaERtbS0AoKysDAUFBbDb7dDpdEhPT0ddXR0A4OLFi3jggQcADP+A\nPPLII+OSPlGs2Wxt2LHjBHp6DgSW9fTsBhDfIRBDxwAOxUhJgUMvUvKZwlmuJO1BS8vPQq6L5xl/\n6Nj2wuH46YSOh2f8FA6HXiT6ysCA7H9iE+rKlTmJDoEIABM/zTJpaUOJDiGsuXP9iQ6BCAATP80y\nVmsetNrdQcu02l1xjyNUDNu3b4x7HEShMPHTrFJYuA5VVRIkaS+A4evqVVX3xj2OUDGwYZeSBRt3\nKfnEqEFzOhptp1puUgcxgt8TGoWNu0REFDUmfiIihWHiJyJSGCZ+IiKFYeInIlIYJn4iIoVh4ici\nUhgmfiIihWHiJyJSGCZ+IiKFYeInIlKYiInf4XDAYDBAr9ejsrIyZBmr1Qq9Xg+TyYSOjo6gdX6/\nHzk5Odi8eXNsIiYioimRTfx+vx/l5eVwOBzo7OxEfX09urq6gsrY7XacP38e3d3dOHLkCLZt2xa0\nvqqqCtnZ2VCFG5aIiIjiSjbxt7e3Q6fTISsrC6mpqSguLkZjY2NQmaamJpSUlAAAzGYzfD4f+vv7\nAQAejwd2ux1PPPEEe+AkIkoSsonf6/UiMzMzMK/RaOD1eqMu8/TTT+PZZ5/FddexKYGIKFnIDlAa\n7eWZsWfzQgg0NzfjlltuQU5ODpxO56QDJIqGzdaG6uoWDAykIC1tCFZrHoDkGfgkVHwcmIUSRTbx\nq9VquN3uwLzb7YZGo5Et4/F4oFar8fvf/x5NTU2w2+24cuUK/v3vf+Oxxx7DsWPHxr1PRUVFYNpi\nscBisUzycEiJbLY27NhxAj09BwLLenp2I1kSf/j4wORPUXM6nbE7iRYyrl69KpYsWSJ6e3vFwMCA\nMJlMorOzM6iMzWYT+fn5QgghTp06Jcxm87j9OJ1OsWnTppDvESEEUqLhsaaGX1HIy9sdtEmozcfu\nKty66SgXLj5J2hP54CZYF6QcU8mdsmf8KSkpqKmpgSRJ8Pv9KC0thdFoRG1tLQCgrKwMBQUFsNvt\n0Ol0SE9PR11dXch98a4emi4DA7If44QLF9+VK3PiHAnRsIjfmPz8fOTn5wctKysrC5qvqamR3cf6\n9euxfv36SYRHFFla2lCiQ5AVLr65c/1xjoRoGG+3oRnPas2DVrs7aJlWuytB0YwXLr7t2zcmKCJS\nOiZ+mvEKC9ehqkqCJO0FAEjSXlRV3ZvgqK4JFx8bdilRVF81EiQuAJWKD3dRsNHtQRP8bKhU1zYJ\nN50s5aI+oBH8ntAoU8mdPOMnIlIYJn4iIoVh4iciUhgmfiIihWHiJyJSGCZ+IiKFYeInIlIYJn4i\nIoVh4iciUhgmfiIihWHiJyJSGCZ+IiKFYeInIlKYiInf4XDAYDBAr9ejsrIyZBmr1Qq9Xg+TyYSO\njg4AwJUrV2A2m7F69WpkZ2fjxz/+cWwjJyKiSZEdgcvv96O8vBytra1Qq9VYu3YtioqKYDQaA2Xs\ndjvOnz+P7u5unD59Gtu2bYPL5cLcuXPxxhtv4Prrr8fQ0BDuvPNOvPXWW7jzzjun/aBo9rPZ2lBd\n3YKBgRSkpQ3Bas2bMf3bz+TYaXaQTfzt7e3Q6XTIysoCABQXF6OxsTEo8Tc1NaGkpAQAYDab4fP5\n0N/fj4yMDFx//fUAgMHBQfj9ftx8883TdBikJDZbG3bsOIGengOBZT09IyNcJX8CDRc7kz/Fi+yl\nHq/Xi8zMzMC8RqOB1+uNWMbj8QAY/h/D6tWrkZGRgbvvvhvZ2dmxjJ0Uqrq6JShxAkBPzwE8//xr\nCYpoYmZy7DQ7yJ7xq0aP/iNj7CgwI9vNmTMH77//Pj777DNIkgSn0wmLxTJu+4qKisC0xWIJWYZo\nxMBA6I/tlStz4hxJ7Mzk2Ck+nE4nnE5nTPYlm/jVajXcbndg3u12Q6PRyJbxeDxQq9VBZebNm4fC\nwkK89957ERM/USRpaUMhl8+d649zJLEzk2On+Bh7Urx///5J70v2Uk9ubi66u7vR19eHwcFBNDQ0\noKioKKhMUVERjh07BgBwuVyYP38+MjIy8K9//Qs+nw8A8MUXX+C1115DTk7OpAMlGmG15kGr3R20\nTKvdhe3bNyYooomZybHTLCEisNvtYunSpUKr1YqDBw8KIYR48cUXxYsvvhgo88Mf/lBotVqxatUq\n8ec//1kIIcTZs2dFTk6OMJlMYuXKleLnP/95yP1HEQIpzfCw4sOvMJqbTwpJ2iMAISRpj2huPhnY\ndPRuQk0nuly42EOKoi5ImaaSO1Vf7SBhpjJSPM1So9uWInw2VKrgIqPnw00nYznZAxzB7wmNMpXc\nySd3iYgUhomfiEhhmPiJiBSGiZ+ISGGY+ImIFIaJn4hIYZj4iYgUhomfiEhhmPiJiBSGiZ+ISGGY\n+ImIFIaJn4hIYWT74ydKJqHGqp0JQy1GwjF4Kd6Y+GlGCD/O7sxOkHLjBzP503ThpR6aEcKNszvT\nzfTxg2lmYuKnGSHcOLsz3WwcP5iSX1SJ3+FwwGAwQK/Xo7KyMmQZq9UKvV4Pk8mEjo4OAMNj9N59\n991Yvnw5VqxYgerq6thFTooSbpzdmW42jh9MM0CkIbqGhoaEVqsVvb29YnBwUJhMJtHZ2RlUxmaz\nifz8fCGEEC6XS5jNZiGEEB999JHo6OgQQgjxn//8RyxdunTctlGEQEoTYrjB5uaTQqvdFbRKq/1x\nUg2pOJly4Y4rMBwjh16kMKaSOyOe8be3t0On0yErKwupqakoLi5GY2NjUJmmpiaUlJQAAMxmM3w+\nH/r7+7Fo0SKsXr0aAHDDDTfAaDTiwoULMf/xotmvsHAdqqokSNJeAIAk7UVV1b0Jjmrqwh0XG3Zp\nOkW8cOr1epGZmRmY12g0OH36dMQyHo8HGRkZgWV9fX3o6OiA2WyORdykQIWF61BYuA4qFeBw/DTR\n4cTMbD0uSl4RE79q9GDPMsSYQX9Hb3fp0iVs2bIFVVVVuOGGG8ZtW1FREZi2WCywWCxRvScRkVI4\nnU44nc6Y7Cti4ler1XC73YF5t9sNjUYjW8bj8UCtVgMArl69igcffBDf//73cd9994V8j9GJn4iI\nxht7Urx///5J7yviNf7c3Fx0d3ejr68Pg4ODaGhoQFFRUVCZoqIiHDt2DADgcrkwf/58ZGRkQAiB\n0tJSZGdnY+fOnZMOkoiIYifiGX9KSgpqamogSRL8fj9KS0thNBpRW1sLACgrK0NBQQHsdjt0Oh3S\n09NRV1cHAHj77bfxyiuvYNWqVcjJyQEAHDp0CPfeO/Mb5YiIZiqVGHtxPt4BqFTj2gdI4Ua3K4X4\nbKhU1xaPnpZbl+zlQs1HUxekXFPJnXxyl4hIYZj4iYgUZnZ2gEKzxmztijmccceLG1CIS4kOi2YZ\nJn5KWjbcMCu7YpYz7njRDcDO5E8xxcZdSj5fNWhKyEUL3g1ZJNGNsdPVuBuKhLVw4D027lIQNu7S\nrDSA9ESHkBSusB4oxpj4KWml4fNEh5AU5rIeKMaY+ClpWfE3aLW7g5ZptbsSFE18jDtefBfb8bcE\nRUOzFRM/Ja1CXJqVXTHLGXe8eJUNuxRzbNyl5BPiadVka4yd7id3A9N8cpfCYOMuERFFjYmfiEhh\nmPiJiBSGiZ+ISGGY+ImIFIaJn5KKzdYGCbmwYD0k5MJma0t0SAnDuqDpElXidzgcMBgM0Ov1qKys\nDFnGarVCr9fDZDKho6MjsHzr1q3IyMjAypUrYxMxzVo2Wxt27DiBFryLk3CiBe9ix44Tik14rAua\nNiKCoaEhodVqRW9vrxgcHBQmk0l0dnYGlbHZbCI/P18IIYTL5RJmszmwrq2tTZw5c0asWLEi5P6j\nCIEUIi9vtxi+WT34JUl7xOiPSbhpuXXJXi7UunB1QSTE1HJnxDP+9vZ26HQ6ZGVlITU1FcXFxWhs\nbAwq09TUhJKSEgCA2WyGz+fDxYsXAQB33XUXbrrpppj/YNHsMzAQupfwK1fmxDmS5MW6oFiImPi9\nXi8yMzMD8xqNBl6vd8JliCJJSxsKuXzuXH+cI0lerAuKhYgDsajCdRI+hhjz6HC02wFARUVFYNpi\nscBisUS9Lc0eVmseenp2Bw1EotXuwvbt9+LEiQQGliBabei6IGVyOp1wOp0x2VfExK9Wq+F2uwPz\nbrcbGo1GtozH44FarY46iNGJn5SrsHB4ZK3nN63FCbwLCWuxvep/AsuVpqpKYl1QwNiT4v37909+\nZ5EaAa5evSqWLFkient7xcDAQMTG3VOnTgU17gohRG9vLxt3KXqA+Or/kKMXRZyeyeXCrgtRF0RC\nTHPjbkpKCmpqaiBJErKzs/G9730PRqMRtbW1qK2tBQAUFBRgyZIl0Ol0KCsrwwsvvBDY/qGHHsId\nd9yBc+fOITMzE3V1dZP/lSIioiljt8yUcDZbG6qrWzAwkIK0tCFYW6qwCf+BwLV+ipOtG+V4dsus\ngoCACrbmk8H1ZM3jpR8Fm0rujHiNn2g6jTy0NboRswfdCYwoOdlww/h66hkerYvJnyaKZ/yUUJK0\nBy0tPwu5jmf8CJzx52EtWvDuuDqSpL1wOH46bjnNfhyIhWascA9tUbABpIdczge6aDKY+Cmhwj20\nRcHS8HnI5XygiyaDiZ8SymrNg1a7O2iZFt9NUDTJy4q/ja8n7S5s374xQRHRTMbETwlVWLgOVVUS\nJGkvgOFr1lV4NcFRJZ9CXBpfT1X3smGXJoWNu5QQ427htOZh06Z1425hZONu6Lpobh5ff/wRUBbe\nzkkzSshbOHt2A2DiihZv7aSp4Bk/xZ3sLZw84w+eDlMXofDWTmXh7Zw0o/AWzunBWzspWkz8FBc2\nWxskaQ8slgp8+GFXosOZZdoA7MHZs72QpD0cnpEi4qkXTbvx1/TbkJLy3xgaejFQRqvdhZ6eg4kJ\ncAYK7qv/BIAD+PRToKWF1/spMp7x07Srrm4JaogE1mFo6GEsWFAM4NqtiRS90bd2AgeC1vX0HMDz\nz78W/6BoxmDjLk2bkVs2T5/24LPP/u+49evXV+DkyYqoGzSTqTE20Y27o8uN14abbnoBq1YZeKvn\nLMbbOSnpBF/e2ROyDLsbmA5tAE7g00+P4+TJ4SW89ENjRbzU43A4YDAYoNfrUVlZGbKM1WqFXq+H\nyWRCR0fHhLal2WWkEffRR18YdXknDwC7G5guwV05tCD40k8benpUePTR/8OGX7pGbniuoaEhodVq\nRW9vrxgcHIw47KLL5QoMuxjNtl9dZprMyGGz0htvvJHoEKLS3HxS5OXtFuvX7xM5OaUiJ2dbYHrR\noqfF8IWIfV/9O/I6KYA9AhBCkvaI5uaTQojhdSOChxt8I2hluCELoxq+cAaUC7suirpobj4pJGm4\nbm+66bEx9b4r6G/wta99V6xYsUPk5e0W+/b9MvB3zMvbHfibJLOZ8h2Jh6nkTtlLPe3t7dDpdMjK\nygIAFBcXo7GxEUajMVCmqakJJSUlAACz2Qyfz4eLFy+it7c34rYjJGkP/uu/FuPUqQsYGEjBv//t\nAfC/cOONtwRNp6UNxa1cPN9rZPqTTzqxbNmdSVcXY6c/+uhGXLz4vzFyWeHaGeYeACMPZo3tdXMd\nRp7Mje4hI2cUZZTCKbu2sHAdCgvXQaUC1q7NREvL6LXX7qQCTuCLLxrw4YfAhx+24fXX/9+oO6va\n8Oabv4RW+wekpl5Coj9n4crNhO9IvMpNidyvwm9/+1vxxBNPBOZffvllUV5eHlRm06ZN4u233w7M\nb9iwQbz33nvid7/7XcRtR361gJMiJaVs1JnhrhDTIo7l4vleo6f3JWFdjC23O8x0qLP8XUHzWu2P\ng85Wh//+IaYBAezjGf9XExOpi+bmk0KrDa730K/Rf7tk/JyFKzcTviPxKjfmAzQBsltGk7w3bdok\n3nrrrcD85BJ/uGQyNrHEq1w832v09L4krItIyT3cPkY+pMGXd8Z+VkfPB6bBxD+Vuhh96Sf8K9zf\nMVk+ZzP5OxKvcmM+QBMgu+WpU6eEJEmB+YMHD4rDhw8HlSkrKxP19fWB+WXLlomLFy9Gta0QQgBa\nMZz8+eIx9iVmAAAFMklEQVSLL774iv6lnXTil73Gn5ubi+7ubvT19WHx4sVoaGhAfX19UJmioiLU\n1NSguLgYLpcL8+fPR0ZGBhYsWBBxWwAQ4rxcCEREFGOyiT8lJQU1NTWQJAl+vx+lpaUwGo2ora0F\nAJSVlaGgoAB2ux06nQ7p6emoq6uT3ZaIiBIr4U/uEhFRfCW0rx4lP+Dldrtx9913Y/ny5VixYgWq\nq6sBAJ988gk2btyIpUuXIi8vDz6fL8GRxoff70dOTg42b94MQLn1AAA+nw9btmyB0WhEdnY2Tp8+\nrcj6OHToEJYvX46VK1fi4YcfxsDAgGLqYevWrcjIyMDKlSsDy+SO/dChQ9Dr9TAYDGgJvp83pIQl\nfr/fj/LycjgcDnR2dqK+vh5dXcrprjc1NRW/+MUv8Ne//hUulwu//OUv0dXVhcOHD2Pjxo04d+4c\nNmzYgMOHDyc61LioqqpCdnY2VF91PqPUegCAHTt2oKCgAF1dXTh79iwMBoPi6qOvrw9Hjx7FmTNn\n8MEHH8Dv9+P48eOKqYfHH38cDocjaFm4Y+/s7ERDQwM6OzvhcDjw1FNP4csvv5R/g0k3C0/RO++8\nE3TXz6FDh8ShQ4cSFU7Cffvb3xavvfZa4K4oIYT46KOPxLJlyxIc2fRzu91iw4YN4vXXXxebNm0S\nQghF1oMQQvh8PnHbbbeNW660+vj444/F0qVLxSeffCKuXr0qNm3aJFpaWhRVD729vWLFihWB+XDH\nPvaOSUmSxKlTp2T3nbAzfq/Xi8zMzMC8RqOB1+tNVDgJ1dfXh46ODpjNZvT39yMjIwMAkJGRgf7+\n/gRHN/2efvppPPvss7juumsfRyXWAwD09vZi4cKFePzxx7FmzRo8+eST+PzzzxVXHzfffDOeeeYZ\n3HrrrVi8eDHmz5+PjRs3Kq4eRgt37BcuXIBGowmUiyaXJizxq0L3J6s4ly5dwoMPPoiqqip8/etf\nD1qnUqlmfT01NzfjlltuQU5ODkSY+wyUUA8jhoaGcObMGTz11FM4c+YM0tPTx13OUEJ99PT04Lnn\nnkNfXx8uXLiAS5cu4ZVXXgkqo4R6CCfSsUeql4QlfrVaDbfbHZh3u91Bv1pKcPXqVTz44IN49NFH\ncd999wEY/iW/ePEiAOCjjz7CLbfcksgQp90777yDpqYm3HbbbXjooYfw+uuv49FHH1VcPYzQaDTQ\naDRYu3YtAGDLli04c+YMFi1apKj6eO+993DHHXdgwYIFSElJwQMPPIBTp04prh5GC/edGJtLPR4P\n1Gq17L4SlvhHPxw2ODiIhoYGFBUVJSqcuBNCoLS0FNnZ2di5c2dgeVFREV566SUAwEsvvRT4QZit\nDh48CLfbjd7eXhw/fhzf+ta38PLLLyuuHkYsWrQImZmZOHfuHACgtbUVy5cvx+bNmxVVHwaDAS6X\nC1988QWEEGhtbUV2drbi6mG0cN+JoqIiHD9+HIODg+jt7UV3dzduv/12+Z3FukFiIux2u1i6dKnQ\narXi4MGDiQwl7t58802hUqmEyWQSq1evFqtXrxavvvqq+Pjjj8WGDRuEXq8XGzduFJ9++mmiQ40b\np9MpNm/eLIQQiq6H999/X+Tm5opVq1aJ+++/X/h8PkXWR2VlpcjOzhYrVqwQjz32mBgcHFRMPRQX\nF4tvfOMbIjU1VWg0GvHrX/9a9tgPHDggtFqtWLZsmXA4HBH3zwe4iIgUhoOtExEpDBM/EZHCMPET\nESkMEz8RkcIw8RMRKQwTPxGRwjDxExEpDBM/EZHC/H+94x7arATSjgAAAABJRU5ErkJggg==\n",
       "text": [
        "<matplotlib.figure.Figure at 0x560e610>"
       ]
      }
     ],
     "prompt_number": 5
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The two vertical lines in the plot show how far out from the mean we have to go to accumulate 95% of the probability. Now, we can solve this as\n",
      "\n",
      "$$ 50 + 50 \\epsilon = 60 $$\n",
      "\n",
      "which makes $\\epsilon=1/5$ or 20%. So, flipping 100 times means I can only get within 20% of the real $p$ 95% of the time in the worst case scenario (i.e. $p=1/2$).\n"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "b=scipy.stats.bernoulli(.5) # coin distribution\n",
      "xs = b.rvs(100) # flip it 100 times\n",
      "phat = mean(xs) # estimated p\n",
      "\n",
      "print abs(phat-0.5) < 0.5*0.20 # did I make it w/in interval 95% of the time?"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "True\n"
       ]
      }
     ],
     "prompt_number": 6
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Let's keep doing this and see if we can get within this interval 95% of the time."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "out=[]\n",
      "b=scipy.stats.bernoulli(.5) # coin distribution\n",
      "for i in range(500): # number of tries\n",
      "    xs = b.rvs(100) # flip it 100 times\n",
      "    phat = mean(xs) # estimated p\n",
      "    out.append(abs(phat-0.5) < 0.5*0.20 ) # within 20% \n",
      "\n",
      "print 'Percentage of tries within 20 interval = %3.2f'%(100*sum(out)/float(len(out) ))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Percentage of tries within 20 interval = 96.20\n"
       ]
      }
     ],
     "prompt_number": 7
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Well, that seems to work. Now we have a way to get at the quality of the estimator, $\\hat{p}$. "
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Summary"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "In this section, we explored the concept of maximum likelihood estimation using a coin flipping experiment both analytically and numerically with the  scientific Python tool chain. There are two key points to remember. First, maximum likelihood estimation produces a function of the data that is itself a random variable, with its own statistics and distribution. Second, it's worth considering how to analytically derive the density function of the estimator rather than relying on canned packages to compute confidence intervals wherever possible. This is especially true when data is hard to come by and the approximations made in the central limit theorem are therefore harder to justify."
     ]
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "References"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "This [IPython notebook](www.ipython.org) is available for [download](https://github.com/unpingco/Python-for-Signal-Processing/blob/master/Maximum_likelihood.ipynb). I urge you to experiment with the calculations for different parameters. As always, corrections and comments are welcome!"
     ]
    }
   ],
   "metadata": {}
  }
 ]
}