2.1-binary-variables.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Populating the interactive namespace from numpy and matplotlib\n"
     ]
    }
   ],
   "source": [
    "%pylab inline\n",
    "from IPython.display import display, Math, Latex\n",
    "maths = lambda s: display(Math(s))\n",
    "latex = lambda s: display(Latex(s))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "class Coin(object):\n",
    "    \"\"\"\n",
    "    A simple coin class, that we can toss to get heads or tails, \n",
    "    it could be represented by a binary random variable x,\n",
    "    where x=1 (heads) or x=0 (tails)\n",
    "\n",
    "    Can have a biased coin, by changing mu, the probability x=1 (heads)\n",
    "    \"\"\"\n",
    "    def __init__(self, mu=0.5):\n",
    "        self.mu = mu\n",
    "        \n",
    "    def toss(self, n=1):\n",
    "        return int_(random_sample(size=n) < self.mu)\n",
    "    \n",
    "    def pmf(self, x):\n",
    "        # p(x=0|mu) = (1-mu)\n",
    "        # p(x=1|mu) = mu\n",
    "        # ie, the Bernoulli distribution\n",
    "        # which can be written as:\n",
    "        return self.mu**x * (1-self.mu)**(1-x)\n",
    "\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([1])"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Now we can toss a coin\n",
    "Coin().toss()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0, 1, 1, 0, 1])"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Now we can toss a coin 5 times\n",
    "Coin().toss(5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Fair Coin   --  Heads: 490 Tails: 510\n",
      "Biased Coin --  Heads: 593 Tails: 407\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "Text(0.5, 0, 'Biased Coin')"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeYAAAD8CAYAAACiqQeGAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAXPElEQVR4nO3df5TV9X3n8ed7BYMWE36IyjKmkIatGjKMZlCjrkXtatSqtMTVrFvBNUtOjjE22awaT3M2Pdk9abN7QuOamkM1FaubxGoaSeqJVYxLEkVFRYISF7SJDlgZUVErRMH3/nG/gwPcYS4z9858gOfjnHvu9/v5fube97365nW/3/vjG5mJJEkqw78a7gIkSdK7DGZJkgpiMEuSVBCDWZKkghjMkiQVZMRwFwBw8MEH5+TJk4e7DKl4jz766EuZOWG469gV+1nq3656uYhgnjx5MsuWLRvuMqTiRcSvh7uG/tjPUv921cseypYkqSAGsyRJBTGYJUkqSBHvMUs93n77bbq6uti8efNwlzKsRo0aRVtbGyNHjhzuUqQBsZdrBtLLBrOK0tXVxUEHHcTkyZOJiOEuZ1hkJhs2bKCrq4spU6YMdznSgNjLA+9lD2WrKJs3b2b8+PH7bCMDRATjx4/f5/c0tGezlwfeywazirMvN3IPnwPtDfz/eGDPgcEsSVJBDGaVLaK5lwaMHj16u/WbbrqJz3zmM015ODNnzvTHN7RvspcbZjBLklQQg1naDd3d3cyePZsZM2YwY8YMfv7znwPw8MMPc8IJJ3D00Udzwgkn8PTTTwOwadMmLrzwQtrb27ngggvYtGkTAFu3bmXu3LlMmzaND3/4w8yfP3/YHpO0Lyq5l/26lLSDTZs20dHRsW395Zdf5txzzwXgiiuu4HOf+xwnnXQSzz33HGeccQarVq3iiCOOYMmSJYwYMYJ7772Xa665hjvuuIPrr7+eAw88kBUrVrBixQqOOeYYAJYvX87atWtZuXIlAK+++urQP1BpL7en9nJDwRwRY4AbgGlAAv8JeBr4HjAZ+BXw7zPzlah9BO0bwFnAm8DczHxs0JVKQ+SAAw5g+fLl29Zvuummbe8l3XvvvTz11FPbtr322mu8/vrrbNy4kTlz5rB69WoigrfffhuAJUuW8NnPfhaA9vZ22tvbAfjABz7As88+y+WXX87ZZ5/N6aefPiSPzV7WvmRP7eVGD2V/A/hxZh4BTAdWAVcDizNzKrC4Wgc4E5haXeYB1w+6SqkQ77zzDg8++CDLly/f9kr5oIMO4ktf+hKnnHIKK1eu5Ic//OF231us93WJsWPH8sQTTzBz5ky++c1v8slPfnKoHoK9LFF2L/cbzBHxXuBk4EaAzHwrM18FzgMWVtMWArOq5fOAm7NmKTAmIiYOulKpAKeffjrXXXfdtvWeV+MbN25k0qRJQO1VeY+TTz6ZW2+9FYCVK1eyYsUKAF566SXeeecdZs+ezVe+8hUee6z1O6L2svSuknu5kT3mDwDdwN9ExOMRcUNE/BZwaGa+AFBdH1LNnwQ83+vvu6qx7UTEvIhYFhHLuru7B/UgtBfLbO5lkK699lqWLVtGe3s7Rx11FN/61rcAuPLKK/niF7/IiSeeyNatW7fN//SnP80bb7xBe3s7X/va1zj22GMBWLt2LTNnzqSjo4O5c+fy1a9+ddC1NaAlvQz2sxpgLzcssp8HGBGdwFLgxMx8KCK+AbwGXJ6ZY3rNeyUzx0bEPwBfzcyfVeOLgSsz89G+7qOzszP9bqcAVq1axZFHHjncZRSh3nMREY9mZudAbm8oehnsZ9XYy+/a3V5uZI+5C+jKzIeq9duBY4AXew5rVdfre80/vNfftwHrGn4EklrFXt4dzf5BjBb+kIb2Lv0Gc2b+M/B8RPxuNXQa8BSwCJhTjc0B7qyWFwEXR83xwMaew2SSho+9LO0ZGv0e8+XArRGxP/AscAm1UL8tIi4FngPOr+beRe3rFWuofcXikqZWLGkw7GWpcA0Fc2YuB+odCz+tztwELhtkXZJawF6WyudPckqSVBCDWZKkghjMKtpQf8B1w4YNdHR00NHRwWGHHcakSZO2rb/11lt1/+aMM87g9ddfZ8uWLYwZM6buHGlfZy83zpNYSL2MHz9+2y8AffnLX2b06NF84Qtf2OXf3H333QBs2bKl5fVJasye3MvuMUsNOuecc/jIRz7Chz70IW644YZt421tbTudUWbt2rWcdNJJdHR0MG3aNB544IGhLldSH0rvZfeYpQYtXLiQcePG8eabb9LZ2cns2bMZO3Zs3bm33HIL55xzDldddRVbt27ddu5WScOv9F42mKUGzZ8/n0WLFgHQ1dXFM888Q2dn/V/HnDFjBp/61KfYvHkzs2bNYvr06UNZqqRdKL2XPZQtNeDee+9lyZIlLF26lCeeeIL29vbtTge3o1NPPZX777+fiRMnctFFF207K42k4bUn9LLBLDVg48aNjBs3jgMOOIAnn3ySRx55ZJfzf/3rX3PYYYcxb9485s6dy+OPPz5ElUralT2hlz2UraI14exuTXH22WezYMECpk+fzhFHHMFxxx23y/mLFy/m61//OiNHjmT06NHccsstQ1SpVCZ7uXH9nvZxKHiaOPXwVHHvavZpH4fKXtPPpZzZqYB/owfCXn5XK077KEmShojBLElSQQxmFaeEt1eGm8+B9gb+fzyw58BgVlFGjRrFhg0b9umGzkw2bNjAqFGjhrsUacDs5YH3sp/KVlHa2tro6uqiu7t7uEsZVqNGjaKtrW24y5AGzF6uGUgvG8wqysiRI5kyZcpwlyFpkOzlgfNQtiRJBTGYJUkqiMEsSVJBDGZJkgpiMEuSVBCDWZKkghjMkiQVpKFgjohfRcQvImJ5RCyrxsZFxD0Rsbq6HluNR0RcGxFrImJFRBzTygegckWUcdH27GepbLuzx3xKZnb0Ok3V1cDizJwKLK7WAc4EplaXecD1zSpWUtPYz1KhBnMo+zxgYbW8EJjVa/zmrFkKjImIiYO4H0mtZz9LhWg0mBP4x4h4NCLmVWOHZuYLANX1IdX4JOD5Xn/bVY1JKoP9LBWs0d/KPjEz10XEIcA9EfHLXcyt967eTqcXqf5BmAfw/ve/v8EyJDWB/SwVrKE95sxcV12vB/4eOBZ4seeQVnW9vpreBRze68/bgHV1bnNBZnZmZueECRMG/ggk7Rb7WSpbv8EcEb8VEQf1LAOnAyuBRcCcatoc4M5qeRFwcfVpzuOBjT2HyCQNL/tZKl8jh7IPBf4+at87GQH8n8z8cUQ8AtwWEZcCzwHnV/PvAs4C1gBvApc0veqSlfL9nH345OTaJftZKly/wZyZzwLT64xvAE6rM57AZU2pTlJT2c9S+fzlL0mSCmIwS5JUEINZkqSCGMySJBXEYJYkqSAGsyRJBTGYJUkqiMEsSVJBDGZJkgpiMEuSVBCDWZKkghjMkiQVxGCWJKkgBrMkSQUxmCVJKojBLElSQQxmSZIKYjBLklQQg1mSpIIYzJIkFcRgliSpIAazJEkFMZglSSqIwSxJUkEMZkmSCtJwMEfEfhHxeET8qFqfEhEPRcTqiPheROxfjb+nWl9TbZ/cmtIlDYS9LJVtd/aYrwBW9Vr/C2B+Zk4FXgEurcYvBV7JzA8C86t5ksphL2u3RZRx2Rc0FMwR0QacDdxQrQdwKnB7NWUhMKtaPq9ap9p+WjVf0jCzl6XyNbrH/JfAlcA71fp44NXM3FKtdwGTquVJwPMA1faN1fztRMS8iFgWEcu6u7sHWL6k3dT0Xgb7WWqmfoM5Iv4AWJ+Zj/YerjM1G9j27kDmgszszMzOCRMmNFSspIFrVS+D/Sw104gG5pwInBsRZwGjgPdSe9U9JiJGVK+k24B11fwu4HCgKyJGAO8DXm565ZJ2l70s7QH63WPOzC9mZltmTgYuBO7LzIuAnwAfr6bNAe6slhdV61Tb78vMuq+yJQ0de1naMwzme8xXAZ+PiDXU3ne6sRq/ERhfjX8euHpwJUpqMXtZKkgjh7K3ycz7gfur5WeBY+vM2Qyc34TaJLWIvSyVy1/+kiSpIAazJEkFMZglSSqIwSxJUkEMZkmSCmIwS5JUEINZkqSCGMySJBXEYJYkqSAGsyRJBTGYJUkqiMEsSVJBDGZJkgpiMEuSVBCDWZKkghjMkiQVxGCWJKkgBrMkSQUxmCVJKojBLElSQQxmSZIKYjBLklQQg1mSpIIYzJIkFaTfYI6IURHxcEQ8ERFPRsSfVeNTIuKhiFgdEd+LiP2r8fdU62uq7ZNb+xAkNcp+lsrXyB7zb4BTM3M60AF8LCKOB/4CmJ+ZU4FXgEur+ZcCr2TmB4H51TxJZbCfpcL1G8xZ80a1OrK6JHAqcHs1vhCYVS2fV61TbT8tIqJpFUsaMPtZKl9D7zFHxH4RsRxYD9wDPAO8mplbqildwKRqeRLwPEC1fSMwvs5tzouIZRGxrLu7e3CPQlLD7GepbA0Fc2ZuzcwOoA04Fjiy3rTqut6r6dxpIHNBZnZmZueECRMarVfSINnPUtl261PZmfkqcD9wPDAmIkZUm9qAddVyF3A4QLX9fcDLzShWUvPYz1KZGvlU9oSIGFMtHwD8PrAK+Anw8WraHODOanlRtU61/b7M3OkVtqShZz9L5RvR/xQmAgsjYj9qQX5bZv4oIp4CvhsR/x14HLixmn8j8LcRsYbaK+sLW1C3pIGxn6XC9RvMmbkCOLrO+LPU3p/acXwzcH5TqpPUVPazVD5/+UuSpIIYzJIkFcRgliSpIAazJEkFMZglSSqIwSxJUkEMZkmSCmIwS5JUEINZkqSCGMySJBXEYJYkqSAGsyRJBTGYJUkqiMEsSVJBDGZJkgpiMEuSVBCDWZKkghjMkiQVxGCWJKkgBrMkSQUxmCVJKojBLElSQQxmSZIKYjBLklQQg1mSpIL0G8wRcXhE/CQiVkXEkxFxRTU+LiLuiYjV1fXYajwi4tqIWBMRKyLimFY/CEn9s5elPUMje8xbgP+SmUcCxwOXRcRRwNXA4sycCiyu1gHOBKZWl3nA9U2vWtJA2MvSHqDfYM7MFzLzsWr5dWAVMAk4D1hYTVsIzKqWzwNuzpqlwJiImNj0yiXtFntZ2jPs1nvMETEZOBp4CDg0M1+AWsMDh1TTJgHP9/qzrmpsx9uaFxHLImJZd3f37lcuacCa2cvV7dnPUpM0HMwRMRq4A/iTzHxtV1PrjOVOA5kLMrMzMzsnTJjQaBmSBqnZvQz2s9RMDQVzRIyk1si3Zub3q+EXew5rVdfrq/Eu4PBef94GrGtOuZIGw16WytfIp7IDuBFYlZlf77VpETCnWp4D3Nlr/OLqE53HAxt7DpNJGj72srRnGNHAnBOBPwZ+ERHLq7FrgD8HbouIS4HngPOrbXcBZwFrgDeBS5pasaSBspelPUC/wZyZP6P+e00Ap9WZn8Blg6xLUpPZy9KewV/+kiSpIAazJEkFMZglSSqIwSxJUkEMZkmSCmIwS5JUEINZkqSCGMySJBXEYJYkqSAGsyRJBTGYJUkqiMEsSVJBDGZJkgpiMEuSVBCDWZKkghjMkiQVxGCWJKkgBrMkSQUxmCVJKojBLElSQQxmSZIKYjBLklQQg1mSpIIYzJIkFaTfYI6Ib0fE+ohY2WtsXETcExGrq+ux1XhExLURsSYiVkTEMa0sXtLusZ+l8jWyx3wT8LEdxq4GFmfmVGBxtQ5wJjC1uswDrm9OmZKa5CbsZ6lo/QZzZi4BXt5h+DxgYbW8EJjVa/zmrFkKjImIic0qVtLg2M9S+Qb6HvOhmfkCQHV9SDU+CXi+17yuamwnETEvIpZFxLLu7u4BliGpCexnqSDN/vBX1BnLehMzc0FmdmZm54QJE5pchqQmsJ+lYTDQYH6x55BWdb2+Gu8CDu81rw1YN/DyJA0B+1kqyECDeREwp1qeA9zZa/zi6tOcxwMbew6RSSqW/SwVZER/EyLiO8BM4OCI6AL+G/DnwG0RcSnwHHB+Nf0u4CxgDfAmcEkLapY0QPazVL5+gzkzP9HHptPqzE3gssEWJak17GepfP7ylyRJBTGYJUkqiMEsSVJBDGZJkgpiMEuSVBCDWZKkghjMkiQVxGCWJKkgBrMkSQUxmCVJKojBLElSQQxmSZIKYjBLklQQg1mSpIIYzJIkFcRgliSpIAazJEkFMZglSSqIwSxJUkEMZkmSCmIwS5JUEINZkqSCGMySJBXEYJYkqSAGsyRJBWlJMEfExyLi6YhYExFXt+I+JA0N+1kaWk0P5ojYD/gmcCZwFPCJiDiq2fcjqfXsZ2notWKP+VhgTWY+m5lvAd8FzmvB/UhqPftZGmIjWnCbk4Dne613AcftOCki5gHzqtXfRMTKFtSyuw4GXhruImhGHRFl1FFADc15Kop4LgB+d4jvb0/t51L+e9nL27Of39VnL7cimOs9bbnTQOYCYAFARCzLzM4W1LJbrKO8OkqoobQ6hvou64wV388l1GAd1tFfDX1ta8Wh7C7g8F7rbcC6FtyPpNazn6Uh1opgfgSYGhFTImJ/4EJgUQvuR1Lr2c/SEGv6oezM3BIRnwHuBvYDvp2ZT/bzZwuaXccAWcf2SqijhBpgH61jD+7nEmoA69iRdbyrzxoic6e3iyRJ0jDxl78kSSqIwSxJUkEMZkmSCmIwS5JUEINZkqSCGMySJBXEYN4LRMTWiFje6zJ5F3P/dUTc3uDtHhYR342IZyLiqYi4KyL+TTNuW9qX9erZJyLisYg4oRpvaQ9FxMyI+FEf246NiCXVKT5/GRE3RMSBu7itcz0NaGv4Pea9QES8kZmjB3kbIzJzS6/1AB4AFmbmt6qxDuCgzPzpoAqW9nG9ezYizgCuyczfG4L7nQl8ITP/YIfxQ4GHgQsz88Gq/2cDP83MF1tdl7bnHvNeKiImR8RPq1fjvV+RT+45809EzI2Iv4uIHwL/uMNNnAK83RPKAJm5PDN/GjX/MyJWRsQvIuKCPm77+xHx44hYHRFfG5IHLu153gu8Ajv1UF89PLHas11e9eC/rcZPj4gHq7l/FxE9wf+xag/4Z8Af9VHDZdRehD8IkDW3Z+aLETEuIn4QESsiYmlEtFe3OzcirquWb4qIayPigYh4NiI+3rqna+/XirNLaegdEBHLq+V/ysw/BNYD/y4zN0fEVOA7QL2zqXwUaM/Ml3cYnwY82sf9/RHQAUyndvq0RyJiSZ15HcDRwG+ApyPif2fm83XmSfuanp4dBUwETq0zp68e/g/A3Zn5PyJiP+DAiDgY+FPg9zPzXyLiKuDz1Qviv65ufw3wvT7qmQYs7GPbnwGPZ+asiDgVuJlab+9oInAScAS131P3ba0BMpj3Dpsyc8dGGQlcVx1+3gr09d7wPXVCuT8nAd/JzK3AixHxf4EZwIod5i3OzI0AEfEU8Ntsf25faV+1rWcj4qPAzRExbYc5ffXwI8C3I2Ik8IPMXB4RvwccBfy8dhSa/YEHqYXkP2Xm6uq+buHd82Y36iRqh7XJzPsiYnxEvK/OvB9k5jvAU9WhcQ2Qh7L3Xp8DXqS2V9tJrVHr+Zc+xp8EPtLHtkZPVf6bXstb8YWgtJPq8PHBwIQdNtXt4cxcApwMrAX+NiIuptaT92RmR3U5KjMv7bmLBsrY3X6vd5u9+73RfyNUh8G893of8EL1CvaPqZ0ZaHfcB7wnIv5zz0BEzKhemS8BLoiI/SJiArV/JB5uUt3SPiUijqDWnxt22FS3hyPit4H1mfnXwI3AMcBS4MSI+GA158DqGxS/BKZExO9Ut/mJPsq4DpgTEcf1qus/RsRh1Pr9ompsJvBSZr42uEetXXEPZu/1V8AdEXE+8BP63jOuKzMzIv4Q+MvqKxGbgV8Bf0KtUT8KPEHtlfOVmfnPsYuvaUnaTu/PhQQwJzO3Voehe/TVwzOB/xoRbwNvABdnZndEzAW+ExHvqeb9aWb+v4iYB/xDRLwE/Iza+8nbqT7kdSHwvyLiEOAdan3+feDLwN9ExArgTWBOU54B9cmvS0mSVBAPZUuSVBCDWZKkghjMkiQVxGCWJKkgBrMkSQUxmCVJKojBLElSQf4/K2rG220gseUAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 576x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Toss a fair coin and a biased coin 1000 times\n",
    "faircoin = Coin(.5)\n",
    "biasedcoin = Coin(.6)\n",
    "\n",
    "nsamples = 1000\n",
    "\n",
    "ftosses = faircoin.toss(nsamples)\n",
    "btosses = biasedcoin.toss(nsamples)\n",
    "\n",
    "bheads = sum(btosses == 1)\n",
    "btails = sum(btosses == 0)\n",
    "\n",
    "fheads = sum(ftosses == 1)\n",
    "ftails = sum(ftosses == 0)\n",
    "\n",
    "print(\"Fair Coin   -- \", \"Heads:\", fheads, \"Tails:\", ftails)\n",
    "print(\"Biased Coin -- \", \"Heads:\", bheads, \"Tails:\", btails)\n",
    "\n",
    "figure(figsize=(8,4))\n",
    "subplot(1,2,1)\n",
    "b = bar([0,1],[fheads,ftails], color=['r','b'], width=.5)\n",
    "b[0].set_label(\"Heads\")\n",
    "b[1].set_label(\"Tails\")\n",
    "gca().axes.xaxis.set_ticklabels([])\n",
    "legend()\n",
    "xlim(-1,2); ylim(0,650)\n",
    "xlabel(\"Fair Coin\")\n",
    "\n",
    "subplot(1,2,2)\n",
    "b = bar([0,1],[bheads,btails], color=['r','b'], width=.5)\n",
    "b[0].set_label(\"Heads\")\n",
    "b[1].set_label(\"Tails\")\n",
    "gca().axes.xaxis.set_ticklabels([])\n",
    "legend()\n",
    "xlim(-1,2); ylim(0,650)\n",
    "xlabel(\"Biased Coin\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Consider $p(x_n|\\mu) = \\mu^{x_n}(1-\\mu)^{1-x_n}$ for each sample $x_n$ in our dataset D"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.6 0.4 0.6 0.6 0.4 0.4 0.4 0.6 0.6 0.4 0.4 0.6 0.4 0.6 0.4 0.4 0.6 0.4\n",
      " 0.6 0.4 0.6 0.6 0.4 0.6 0.6 0.4 0.6 0.6 0.6 0.6 0.6 0.4 0.6 0.4 0.4 0.6\n",
      " 0.4 0.6 0.4 0.6 0.6 0.6 0.6 0.4 0.6 0.6 0.6 0.6 0.6 0.6]\n"
     ]
    }
   ],
   "source": [
    "# Show first 50\n",
    "print(biasedcoin.pmf(btosses)[:50])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can find the *likelihood* of all these independent samples by multiplying\n",
    "\n",
    "$p(D|\\mu) = \\prod_{n=1}^{N} p(x_n|\\mu)$\n",
    "\n",
    "And taking the log to give the log likelihood, which is nicer numerically:\n",
    "\n",
    "$\\log(p(D|\\mu)) = \\sum_{n=1}^{N} \\log(p(x_n|\\mu))$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Add likelihood and loglikelihood to our coin class\n",
    "Coin.likelihood = lambda self, data: product(self.pmf(data))\n",
    "Coin.loglikelihood = lambda self, data: sum(log(self.pmf(data)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "3.0346453801326247e-294\n",
      "-675.8499227660134\n"
     ]
    }
   ],
   "source": [
    "print(biasedcoin.likelihood(btosses))\n",
    "print(biasedcoin.loglikelihood(btosses))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, in order to estimate $\\mu$ from our samples we want to maximise the log likelihood -- or minimise the negative log likelihood. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x112b9d7b8>"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYsAAAEQCAYAAABBQVgLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nOzdd3zU9f3A8dc7OyGDbCABElZYslFQxI24UcFRB7ZW21psa1ttra3aWrurtb+qFasV68CtqCiiiIiyJewVRiBk7z3v8/vje4EAGZfkvrlL8n4+Hve43Hfd+xvCve+zxRiDUkop1RofTweglFLK+2myUEop1SZNFkoppdqkyUIppVSbNFkopZRqk5+nA7BDTEyMSUpK8nQYSinVrWzatCnfGBPb3L4emSySkpLYuHGjp8NQSqluRUTSW9qn1VBKKaXapMlCKaVUmzRZKKWUalOPbLNQSnUPdXV1ZGRkUF1d7elQepWgoCASExPx9/d3+RxNFkopj8nIyCAsLIykpCRExNPh9ArGGAoKCsjIyCA5Odnl87QaSinlMdXV1URHR2ui6EIiQnR0dLtLc5oslFIepYmi63Xkd67JoqnqElj5J8jY5OlIlFLKq2iyaMo4YOUf4cg6T0eilOoivr6+TJgw4djj0KFDLR6bmZnJ3LlzXbpudnY2N9xwA0OHDmX06NFceuml7N271y3X9gRt4G4qMALEFyrzPR2JUqqLBAcHk5qa6tKxAwYM4M033zxle319PX5+xz9OjTFcffXVzJ8/n8WLFwOQmppKTk4OI0aMaNe1vYWWLJry8YGQKKgs8HQkSikPOnToEGeffTaTJk1i0qRJfP3118e2jx07FoAXXniBefPmccUVVzBr1qwTzv/888/x9/fn+9///rFtEyZM4Oyzz8YYw7333svYsWM57bTTeO2115q99jXXXMPs2bMZPnw49913X1fcdqu0ZHGykBhNFkp5wG/f38HOzFK3XnP0gHAeumJMq8dUVVUxYcIEAJKTk3nnnXeIi4tj+fLlBAUFsW/fPm688cZm55tbs2YNW7duJSoq6oTt27dvZ/Lkyc2+39tvv01qaipbtmwhPz+fqVOnMnPmzFOOS01NZfPmzQQGBpKSksLdd9/NwIEDXb11t9NkcbKQaKjQZKFUb9FcNVRdXR0LFiwgNTUVX1/fFtsaLrroolMSRVtWr17NjTfeiK+vL/Hx8Zxzzjls2LCBcePGnXDcBRdcQEREBACjR48mPT1dk4VXCYmCvD2ejkKpXqetEkBXevzxx4mPj2fLli04HA6CgoKaPa5Pnz7Nbh8zZkyL7Q/GGJdiCAwMPPazr68v9fX1Lp1nF22zOFkfrYZSqrcrKSmhf//++Pj48L///Y+GhoZ2nX/++edTU1PDs88+e2zbhg0b+OKLL5g5cyavvfYaDQ0N5OXlsWrVKk4//XR334LbabI4WUg0VBWCw+HpSJRSHnLXXXexaNEipk2bxt69e1ssQbRERHjnnXdYvnw5Q4cOZcyYMTz88MMMGDCAq6++mnHjxjF+/HjOP/98/vKXv9CvXz+b7sR9xNUiUXcyZcoU0+HFj9Y+DR//Eu47aFVJKaVss2vXLkaNGuXpMHql5n73IrLJGDOlueO1ZHGykBjruULHWiilVCNNFidrLE1ou4VSSh2jyeJkfZwlC00WSil1jCaLk4VEW8865YdSSh2jyeJkx5KFliyUUqqRJouT+QeDfx+oLPR0JEop5TU0WTQnJFp7QynVSzROUT5+/PgTJg20e8rwlStXcvnllze7b/369cycOZOUlBRGjhzJd7/7XSorK1u81pIlS/jTn/5kV6iATvfRPJ15Vqleo+ncUMuWLeP+++/niy++8NiU4Tk5OcybN4/Fixczffp0jDG89dZblJWVERIS0uw5V155JVdeeaWtcWnJojk65YdSvVJpaSmRkZHAiVOGtzRleVZWFjNnzmTChAmMHTuWL7/8EoBPPvmE6dOnM2nSJObNm0d5eTkAH3/8MSNHjmTGjBm8/fbbzcbw5JNPMn/+fKZPnw5Yo8Hnzp1LfHw8hYWFzJkzh3HjxjFt2jS2bt0KWFOaL1iwAIDbbruNH/3oR5x55pkMGTLEbQlPSxbNCYmG/JZXtFJK2eCjX0L2Nvdes99pcEnr1TONU5RXV1eTlZXFihUrTjmmpSnLX3nlFS6++GIeeOABGhoaqKysJD8/n9///vd8+umn9OnThz//+c889thj3Hfffdxxxx2sWLGCYcOGcf311zcbz/bt25k/f36z+x566CEmTpzIu+++y4oVK7j11lubXbgpKyuL1atXs3v3bq688kq3VKdpsmhOSLQ2cCvVSzSthlqzZg233nor27dvP+GYlqYsnzp1Kt/5zneoq6tjzpw5TJgwgS+++IKdO3dy1llnAVBbW8v06dPZvXs3ycnJDB8+HICbb76ZhQsXtivW1atX89ZbbwHWZIUFBQWUlJScctycOXPw8fFh9OjR5OTktO8X0gLbkoWIDAReBPoBDmChMeYJEXkYuAPIcx76K2PMUuc59wO3Aw3Aj4wxy5zbZwNPAL7Af4wx9rbkhERDbTnUVYN/81MTK6XcrI0SQFeYPn06+fn55OXlnbC9pSnLZ86cyapVq/jwww+55ZZbuPfee4mMjOSiiy7i1VdfPeEaqampiEibMYwZM4ZNmzZx1VVXnbKvubn8mrtm0+nN3TX/n51tFvXAz4wxo4BpwA9FZLRz3+PGmAnOR2OiGA3cAIwBZgNPiYiviPgCTwKXAKOBG5tcxx461kKpXmn37t00NDQQHR19wvaWpixPT08nLi6OO+64g9tvv51vvvmGadOm8dVXX5GWlgZAZWUle/fuZeTIkRw8eJD9+/cDnJJMGi1YsIBFixaxbt26Y9teeuklsrOzmTlzJi+//DJg9aaKiYkhPDzc7b+H5thWsjDGZAFZzp/LRGQXkNDKKVcBi40xNcBBEUkDGid5TzPGHAAQkcXOY3faFfsJySKitZCVUt1d02VVjTEsWrQIX1/fE4656667uPbaa3njjTc477zzjk1ZvnLlSv7617/i7+9PaGgoL774IrGxsbzwwgvceOON1NTUAPD73/+eESNGsHDhQi677DJiYmKYMWPGKdVdAPHx8SxevJif//zn5Obm4uPjw8yZM7nmmmt4+OGH+fa3v824ceMICQlh0aJFNv92juuSKcpFJAlYBYwFfgrcBpQCG7FKH0Ui8i9grTHmJec5zwEfOS8x2xjzXef2W4AzjDELTnqPO4E7AQYNGjQ5PT294wGnfw3/vQRueQeGnt/x6yilWqVTlHuO101RLiKhwFvAT4wxpcDTwFBgAlbJ4++NhzZzumll+4kbjFlojJlijJkSGxvbuaCPlSy0kVsppcDm3lAi4o+VKF42xrwNYIzJabL/WeAD58sMoOlq5IlApvPnlrbbI0RnnlVKqaZsK1mI1UT/HLDLGPNYk+39mxx2NdBYabcEuEFEAkUkGRgOrAc2AMNFJFlEArAawZfYFTcAwX0B0Sk/lOoCPXG1Tm/Xkd+5nSWLs4BbgG0i0jhq5FdYvZkmYFUlHQK+B2CM2SEir2M1XNcDPzTGNACIyAJgGVbX2eeNMTtsjBt8fCE4UksWStksKCiIgoICoqOjXepWqjrPGENBQcGx7r+usrM31Gqab29Y2so5jwKPNrN9aWvn2UKn/FDKdomJiWRkZJwyrkHZKygoiMTExHadoyO4WxISrclCKZv5+/uTnJzs6TCUC3QiwZZoslBKqWM0WbREk4VSSh2jyaIljclCe2oopZQmixaFRIOjHqpPndFRKaV6G00WLemjA/OUUqqRJouW6MyzSil1jCaLlmiyUEp1Mz9evJnbX9hgy7U1WbSkMVnolB9KqW4ivaCS2gaHLdfWZNESLVkopbqZvLIaYkMD2z6wAzRZtCSgD/iHQLl71q9VSik7GWPIK68hNkyTRdcSgcgkKDzo6UiUUqpNpdX11NY7NFl4RGQyFGmyUEp5v7wyawlXTRaeEJUMRYfAYU+DkVJKuUt+uTNZaJuFB0QmQX01lGd7OhKllGpVY8kiRksW9jPGUFJVR0VNvbUhyjl1srZbKKW83LFqKJtKFi2uZyEiP23txKZLpfYU2aXVTP/jCn4/Zyw3TxtstVmA1W6RdJZng1NKqVbkldfg7ytEBPvbcv3WFj8Kcz6nAFM5vu71FcAqW6LxsJjQQEQg15mh6TsIxFdLFkopr5dfVkNMaCA+PvYsT9tisjDG/BZARD4BJhljypyvHwbesCUaD/P39SEqJOBYcQ5ff+g7UHtEKaW8Xl65lSzs4kqbxSCgtsnrWiDJlmi8QGxYIHll1cc3RCZryUIp5fXyyuwbkAeurcH9P2C9iLzjfD0HWGRbRB4WFx50vBoKrEbuHe+0fIJSSnmBvLIaxg6IsO36bSYLY8yjIvIRcDZggG8bYzbbFpGHxYYGsje77PiGyGSoKrIewZGeC0wppVrgcBgKKmptLVm42nW2AXA0efRYceGB5JfX4HA4l1PV7rNKKS9XVFlLg8MQExpg23u0mSxE5MfAy0AMEAe8JCJ32xaRh8WFBVLvMBRVOptpmnafVUopL5TXOHo7LMi293ClzeJ24AxjTAWAiPwZWAP8n21ReVCc85edW1ZDdGigNYobtGShlPJads8LBa5VQwlWNVSjBue2Hqnxl32skTswFPrEaclCKeW1uiJZuFKy+C+wztkbSoCrgOdsi8jD4py/7LyTe0QVHvJMQEop1YZj80LZ2GbhSm+ox0RkJTDDualH94aKC28sWZw01uLQlx6KSCmlWpdfXkOQvw+hga58/++Y9vSGMvSC3lAhAX6EBvqRW3pSyaI0E+qqWz5RKaU8pHFAnoh9LQTaG6oZcWGBJ1ZDRSYDBorTPRaTUkq1JK/cvrW3G2lvqGbEnJwsooZYz4UHITbFM0EppVQL8spqSIruY+t7aG+oZsSFBZ7YZhGlYy2UUt4rv9ze0dvQ/t5QYM0N1WN7Q4E11iK3LPf4hpBoCAiDwgOeC0oppZpR1+Cg0OapPsD13lBfAGdhlSh6dG8osHpEVdY2UF5Tb/UuEIHYEZC7y9OhKaXUCQrKrdkmPJ4snFKBrMbjRWSQMeawbVF5WONYi9zSakJjQ62N/cfDtrfAGCt5KKWUFzg+xsLeZOFKb6i7gRxgOfAB8KHzuceKbW5gXr9xUFMCRYc8E5RSSjUjr9xqX/WGksWPgRRjTIGtkXiRpvNDHdN/vPWcteV4g7dSSnlYfpmzGsrTJQvgCFBiaxReJu7k+aEA4kaDjx9kb/VQVEopdarjM856qGQhIj91/ngAWCkiHwLHPj2NMY/ZGpkH9Q3xJ8DX58Tus/5BEDvKKlkopZSXyCurISzQjyB/X1vfp7WSRZjzcRirvSKgybawti4sIgNF5HMR2SUiO5wjwRGRKBFZLiL7nM+Rzu0iIv8UkTQR2Soik5pca77z+H0iMr/jt+saEXGuxV1z4o7+4yAz1WrkVkopL2D32tuNWixZGGN+28lr1wM/M8Z8IyJhwCYRWQ7cBnxmjPmTiPwS+CXwC+ASYLjzcQbwNHCGiEQBDwFTsOan2iQiS4wxRZ2Mr1WnjOIGq90i9WUoy4LwAXa+vVJKuSSvrIaYLkgWLZYsROQfzuf3RWTJyY+2LmyMyTLGfOP8uQzYBSRgTXG+yHnYIqxBfji3v2gsa4G+ItIfuBhYbowpdCaI5cDsDt1tO8SFBZ44mSCc2MitlFJeILOkivhw+1bIa9Rab6j/OZ//1tk3EZEkYCKwDog3xmSBlVBEJM55WAJWY3qjDOe2lraf/B53AncCDBo0qLMhExcWyMZDhSdujB8LCGRthZRLOv0eSinVGVW1DRwtrmLu5ETb36u1aqhNzucvOvMGIhIKvAX8xBhT2soUus3tMK1sP3GDMQuBhQBTpkzpdKNCXFgQRZV11NY7CPBzFsACQyFmuJYslFJeYX9eOcbAsLhQ29+rtd5Q22jmQxnrw9sYY8a1dXER8cdKFC8bY952bs4Rkf7OUkV/oHESpgxgYJPTE4FM5/ZzT9q+sq337qzGBqP88hoG9A0+vqPfODi81u63V0qpNu3PKwc8nCyAyztzYbGKEM8Bu07qZrsEmA/8yfn8XpPtC0RkMVYDd4kzoSwD/tDYawqYBdzfmdhc0XSsxQnJov942P4mVBRAn2i7w1BKqRal5ZbjI5AcY+/05NB6NdSxlX5EZDAw3BjzqYgEt3ZeE2cBtwDbRCTVue1XWEnidRG5Hatb7jznvqXApUAaUAl82xlHoYg8AmxwHvc7Y8xJjQnud2x51dKTVsdrbOTO3gJDz7c7DKWUalFabjmDokII9LN3jAW48KEvIndgNRxHAUOxqoH+DVzQ2nnGmNW0vO7FKecaYwzwwxau9TzwfFuxulPjlB+NoyOP6e+sfcvSZKGU8qy03PIuqYIC16b7+CFWKaEUwBizD2t51R4tJjQAEcg5uftscCT0HWQNzlNKKQ+pb3BwqKCCoV6ULGqMMbWNL0TEj+YbvnsUP18fEiODOeBsQDpB4lSrkVtHciulPCS9sJK6BsOwWO9JFl+IyK+AYBG5CHgDeN/esLzD8Lgw0nKbSRZJM6A8Gwr2d31QSikFxz6bvKka6pdAHrAN+B6w1BjzgK1ReYnh8aEcyKugvsFx4o6kmdbzoVVdH5RSSnE8WXhTNdREY8yzxph5xpi5xphnReQK2yPzAiPiwqhtcHCooPLEHdFDIaw/HFrtmcCUUr3e/txy4sMDCQ/yP75x+9uw9XVb3s+VZPGsiJzW+EJEbgR+bUs0XmZEvDW57r6cshN3iFhVUQe/1HYLpZRHpOU10xNqw3Ow0Z6Oo64ki7nAIhEZ5exGexfWwLgeb1hcKCKwr9l2i7OhIhfy93Z9YEqpXs0Yw/7c8lMbt0uOQIQ980S1Oc7CGHNARG4A3sWa0G+WMabKlmi8THCAL4mRwew9uWQBVskC4NCXEJvStYEppXq1rJJqKmobTixZOBqgNLPrk0Uzc0NFAb7AOhHBlbmheoIRcWHsy2mmZBE1BMITrKqoqd/t+sCUUr3W8Z5QTdahK88FR51HShadmhuqpxgeH8aqfXnUNTjw921SaydiVUWlfWq1W7Q8m65SSrlVs91mSzKs54iBzZzRea21WRQ554cqa+HRK4yID6WuwZBeUHHqzqQZUJkPebu7PjClVK+VlldORLA/MaEBxzeWOJf98UDJ4hWs0sUmTl1XwgBDbInIyxzvEVV+YpEPIPls6/nglxA3qosjU0r1Vo1zQp2wPtCxkoU9yaLFkoUx5nLnc7IxZojzufHRKxIFwNBYq0fU3ubaLSKTIGKQ1citlFJdwBjDvpyyZnpCZUBgOARF2PK+rTVwT2rtxMb1tXu64ABfBkaGsDe3hZq35Jmw631oqANf/+aPUUopN0kvqKSoso5xA09KCiUZtpUqoPVqqL+3ss8AvWZ+7hHxoacOzGs08lJIfckazT30vK4NTCnV62xKLwJg8uDIE3fYOMYCWl/8SD/5nIbHh/HF3mZ6RAEMOQ/8gmHPUk0WSinbbUwvIizQjxEnt6GWZEDiFNve15UR3L1eqz2iAkJg2AWw+0Od+kMpZbtv0ouYODgSH58mjdu1FVBVaGvJQpOFC4Y7M3izjdwAIy+D0qOQpQsiKaXsU1JVx97cMqacUgV11Hq2aYwFaLJwyfEeUS20W4yYDeJjlS6UUsommw8XYUwL7RXgsQZuoMVeUSVAujGm3v0heZ/gAF+SovuwI7O0+QNComDQmVayOL9XTMirlPKATelF+AhMGNj3xB02j7EA10oWTwFrgYXAs8AaYDGwV0R6xeyzYGXyTelFmJbaJUZeBrk7ofBA1wamlOo1NqUXMap/OH0CT/qeX5Jh1W6E9bftvV1JFoewFkCaYoyZDEwEtgMXAn+xLTIvc3pSFIUVtexvbk1usLrQAuxe2nVBKaV6jfoGB6lHik9trwArWYT1t3WslyvJYqQxZkfjC2PMTqzk0au+Qk9NjgJg/cGi5g+ITIL402D3B10XlFKq19idXUZlbQOTmk0W9o6xANeSxR4ReVpEznE+nsKqggoE6myNzoskRYcQExrIhkOFLR806go4vPZ4/aFSSrnJRudnz5SkqFN32jx6G1xLFrcBacBPgHuAA85tdUCvGYUmIpyeHMn6g60ki3HXAQa2vtZlcSmleodNh4vpFx7EgIigE3c4HFbXfU8nC+eqeP8HPIi19vYTxphKY4zDGNNCBX7PNDUpiqPFVRwtbmGhwKhkq1dU6qs6QE8p5VbfpBcxOSnyxJlmASryoKHW1jEW4EKyEJFzgX3Av7B6Ru0VkZm2RuWlpjqLfxtaK11MuBEK9kHGxi6KSinV06UXVHC0uIqpLTVug+dLFlgTCs4yxpxjjJkJXAw8bmtUXmpU/3DCAv1Y31q7xeg51lxRW17pusCUUj3aZ7tyAThvZNypO7tgQB64liz8jTF7Gl8YY/YCvXIubl8fYdLgyNZLFkHhMOpy2P4W1FV3XXBKqR5rxe5chsWFMji6z6k7vahksVFEnhORc52PZ7FWz+uVTk+OYl9uOUUVtS0fNP5GqC6BvR91XWBKqR6prLqOdQcLuGBUM6UKsJJFQCgE9W1+v5u4kix+AOwAfgT8GNgJfN/OoLzZsXaL1qqihpwLYQOshm6llOqEL/flU9dguGBkfPMHNI6xOLnh281c6Q1VY4x5zBhzjTHmamPM48aYGluj8mLjEiMI8PVpPVn4+ML46yHt0+OzQSqlVAd8tiuXiGB/Jg1qoeRQcgTCE2yPo8VkISLbRGRrSw/bI/NSQf6+TBjUl9VpBa0fOPk2wMDG57oiLKVUD9TgMKzck8u5KbH4nbzwGkB9LeTuhrhRtsfS2qyzl9v+7t3UhaPi+MPS3RwprGRgVEjzB0UmQcqlsPG/MPNe8A/u0hiVUt3floxiCipquWBUC1VQebugoQYGTLQ9lhZLFsaY9NYetkfmxS4a3Q+AT3fltH7gGd+zVq/a9kYXRKWU6mlW7MrF10c4Z3hs8wdkbraePZksVMuSY/owPC6UT3a0kSySzob4sbD23zqiWynVbp/uymHK4EgiQloYrZC5GQIjIGqI7bFosuigWWPiWX+okOLKVrrQilili9wdcGh11wWnlOr2jhRWsju7rOUus2AliwETbO8JBS4mCxEJFpEUu4PpTi4a3Y8Gh2HF7tzWDzxtHgRHwbp/d01gSqke4d3NVk/KS09rYUGjumrI2dklVVDg2txQVwCpwMfO1xNEZIndgXm7cQkRxIcHtl0V5R9s9Yza/SEU7O+S2JRS3Zsxhrc3H2XakCgSI1voRJO7Axx1kNDcytfu50rJ4mHgdKAYwBiTCiS1dZKIPC8iuSKyvcm2h0XkqIikOh+XNtl3v4ikicgeEbm4yfbZzm1pIvJL12/NXj4+wkWj41m1L4/quobWDz7je+AXCKsf65rglFLd2uYjxRzMr+CaSa1M4dGFjdvgWrKoN8aUdODaLwCzm9n+uDFmgvOxFEBERgM3AGOc5zwlIr4i4gs8CVwCjAZudB7rFS4a3Y/K2ga+Sstv/cCwflbpYstiKDrUFaEppbqxtzZlEOTv03IVFFjJIiTa9qnJG7mSLLaLyLcAXxEZLiL/B3zd1knGmFVAK8OcT3AVsNg5Wvwg1mJLpzsfacaYA8aYWmCx81ivMH1INGGBfm1XRQGc9RMQX/jy7/YHppTqtmrqG3h/Syazx/QjNLCVoXBHN1ulii5o3AbXksXdWN/4a4BXgBKsVfM6aoFzFPjzItI4OXsCcKTJMRnObS1tP4WI3CkiG0VkY15eXifCc12Anw/njYxj2c5saurbqIoK7w+TboXUV6CoVw9TUUq1YsWuXEqr61uvgqqttAbkdVEVFLiWLFKMMQ8YY6Y6H782xnR07u2ngaHABCALa60MgOZSo2ll+6kbjVlojJlijJkSG9vCABYbXDs5keLKOj7d2UavKIAZ94D4wOpeuRyIUsoFb32TQXx4IGcNi2n5oOxtYBxelyweE5HdIvKIiIzpzJsZY3KMMQ3GGAfwLFY1E1glhqYVb4lAZivbvcaMYTEMiAjitY1H2j44IgEm3gKbX4Liw/YHp5TqVvLLa1i5J485ExPw9WmleqmLG7fBtVlnzwPOBfKAhc4JBn/dkTcTkaatNVcDjT2llgA3iEigiCQDw4H1wAZguIgki0gAViO4V3Xb9fUR5k5O5Mt9eWS2tDZ3U2f/1CpdrPi9/cEppbqVxesPU+8wzJvcxkJGmZshNB7CWmkAdzOXBuUZY7KNMf/EWsciFXiwrXNE5FVgDZAiIhkicjvwl8bZbIHzgHuc198BvI61VsbHwA+dJZB6YAGwDNgFvO481qvMmzIQY+DNTRltHxyRCNPvgq2vwdFv7A9OKdUt1NY7eHFNOmcPj2FYXFjrB2d+06WN2+DaoLxRzvER24F/YfWEanP9PmPMjcaY/sYYf2NMojHmOWPMLcaY04wx44wxVxpjspoc/6gxZqgxJsUY81GT7UuNMSOc+x7t4H3aamBUCGcOjeb1jUdwOFyYA2rGTyEkBj75jc4ZpZQC4IOtmeSW1XD7jOTWDyzNgvy9MPD01o9zM1dKFv8FioBZxphzjDFPG2NcaM3tXa6fOpCMoirWHmhjnQuw1uk+735IXw17ltofnFLKqxljeG71QYbFhXLOiDY66KQtt56HX9z6cW7mSpvFNGPME8YYr2pY9jYXj+lHWJCfaw3dAJNug5gRsPxBaKizNTallHdbd7CQHZmlfOesZKStqqW9y6yV8eI71d+o3VpbKe915/PJK+Y1tjmoJoL8fbl6YgIfbc8mr8yFVWd9/eCiR6AgDdY9Y3+ASimv9fzqg0SG+HPNpDaWR62vgQMrYfhFXdpeAa2XLH7sfL4cuKLJo/G1OsltZyZR1+Dgha8PunbCiIutouTnf4BiF0skSqkeJb2gguW7crjpjMEE+fu2fvDhNVBb3uVVUND6SnmNjc93NbNK3l1dE173MiQ2lNlj+vG/NemU19S3fYIIXPpXa8ccV+wAACAASURBVHDNR7+wP0CllNd58vM0/H19uGX64LYP3vsJ+AbAkHPsD+wkrjRwX9TMtkvcHUhPcefMIZRW17N4vYuD7iIHw7m/hD0fwq4P7A1OKeVVDuSV89Y3R7n5jMHEhwe1fcK+TyBpBgT0sT+4k7TWZvEDEdmGNU6iaZvFQUDbLFowcVAkZyRH8dzqg9Q1OFw7afoPIW4MfHQf1JTZG6BSyms88dk+Anx9+MG5Q9s+uPAAFOzzSBUUtF6yeAWrbWIJJ7ZZTDbG3NwFsXVb3z93KFkl1by/xcUOZL7+cMU/oDQTPn3Y1tiUUt5hb04ZS7ZkMv/MJGLDAl044RPreXhzlT32a63NosQYc8g5uC4dqMKaxC9URAZ1WYTd0LkjYkmJD+OZLw64NkgPrAE2034AG/4DaZ/ZG6BSyuMeX76XPgF+fG/mENdO2PcJRA+DaBdKITZwaVlVEdkHHAS+AA4BH7V6Ui8nItx13lD25JTxwbastk9odMGDEJMC7y2AqiL7AlRKedT2oyV8tD2b22ckE9knoO0Tqkvh0GqPVUGBaw3cvwemAXuNMcnABcBXtkbVA1wxbgCj+ofzt2V7qK13se3CPxiueQYqcmHpvfYGqJTyCGMMv/tgJ5Eh/tx+dhtTezTa8TY01MDYa+0NrhWuJIs6Y0wB4CMiPsaYz7HWo1Ct8PERfjE7hcOFlbzqas8osCYHm3kfbHsDtr9lX4BKKY94f2sW6w8W8vOLUwgP8nftpM0vQewoSJhkb3CtcCVZFItIKLAKeFlEngBcGESgzhkRy/Qh0fzzs32UVbdjSo+zfwoJU+D9n1g9IJRSPUJFTT2PfriTsQnh3DDVxabf3N2QsQEm3tzlo7abciVZXIXVuH0P1vTh+9ER3C4REX55yUgKKmp59ksXR3WD1Ttq7vPWuhevz4e6ji5MqJTyJv/6PI2c0hp+e+XY1hc3air1JfDxg3HX2xtcG1yZSLCicW0JY8wiY8w/ndVSygXjB/blstP6858vD5BT2o4P/cjBcPW/IXsrLLvfvgCVUl3iQF45//nyANdOSmTy4EjXTmqogy2LYcRsCO265aKb40pvqDIRKT3pcURE3hERF/t89W73zU6h3mH43fs723diyiVw5o9g4/Ow7U17glNK2c7hMDzwznaC/Hz5xSUprp+47xOoyLOWY/Ywl9bgBu4FErAWPfo51vrZi4Hn7Qut5xgc3YcF5w3jw21ZfL6nnUuBXPAgDJpudafNTLUnQKWUrV5el86aAwX86rJRxIW5MK1Ho80vWcunDrvQvuBc5EqymG2MecYYU2aMKTXGLAQuNca8BrhYllLfO2cIQ2L78Jt3t1NV2+D6ib7+cN2LEBINi78FZTn2BamUcrvDBZX88aPdnD08hhumDnT9xNIsa+2K8TdaSxp4mCvJwiEi14mIj/NxXZN9uiaoiwL9fHl0zmlkFFXxfyv2te/k0Di48RVroN5rN1tz2iulvJ7DYbj3zS34ivDna8e1vbBRU2v+ZT1Pvs2W2NrLlWRxE3ALkAvkOH++WUSCgQU2xtbjTB8azbWTElm46gC7skrbd3L/8TDnachYD0t+pGt3K9UNvLjmEOsOFvKby0czoG+w6ydWFsLG/8JpcyHKxYF7NnOlN9QBY8wVxpgYY0ys8+c0Y0yVMWZ1VwTZkzxw2Sj6hgRwz2upVNe1ozoKYMwcOO8B2LoYPvudPQEqpdxi+9ES/vDRbs4fGce8KYntO3ndM1BXATPusSe4DnClN9QIEflMRLY7X48TkV/bH1rPFNUngL/OHcfu7DL+tmxP+y8w816YNB9WPwbrFro/QKVUp5VV17HglW+ICgngb/PGt6/6qaYM1v0bRl4OcaPsC7KdXKmGeha4H6gDMMZsBW6wM6ie7ryRcdwybTD/WX2Q1fvy23eyCFz2GKRcaq1/sfM9e4JUSnWIMYb7397GkaIq/u9bE4lyZaLApja9ANXFMOOntsTXUa4kixBjzPqTtul0H530q0tHMTS2Dz97I5Wiitr2nezrB9c+B4lT4a3vQtqn9gSplGq3l9Yd5oOtWfx8VgpTk6Lad3JdNXz9L0g+BxIn2xNgB7mSLPJFZCjOnk8iMhdox7zbqjnBAb48ccNECitquef1VNfXvWgUEAI3vQ6xKbD4Jjj4pT2BKqVctu5AAb97fwfnpcS6vk5FU+sXQnk2zPy5+4PrJFeSxQ+BZ4CRInIU+AnwA1uj6iXGJkTw0BVjWLknj3981s7utADBkXDLuxCZBK9cD4fXuj1GpZRrDhdU8v2XNjEwKoR/3DARH1fnfmpUnger/mqtWZE8054gO8HV3lAXArHASGPMDGPMIdsj6yVuOmMQ8yYn8s/P9vHpzg4MuOsTA7cugfD+8NJcSF/j/iCVUq0qq67j9kUbcBh4bv5UIoJdnHq8qRWPQF0lXPyo+wN0A1d6QwWKyLeAHwP3iMiDIvKg/aH1DiLCI3PGMjYhnHteS+VAXnn7LxIWD/Pfh7B+8NI1sP9z9weqlGpWXYODu1/dzMH8Cp6+aRLJMX3af5GsrfDNi3D6nRAz3P1BuoEr1VDvYU1TXg9UNHkoNwny9+XfN0/G38+H77ywgYLyDozQDh8A314KkclWldTeZe4PVCl1AofDcN+bW1m5J49H5ozlzGEx7b+IMfDx/Va18jn3uT9IN3ElWSQaY643xvzFGPP3xoftkfUyiZEhPHvrFLJKqvnuixvbN39Uo9A4uO0DiB9tzSO15TX3B6qUAqwusr//cBfvbD7Kz2eN4MbTXVzM6GTb34L01XD+A1bC8FKuJIuvReQ02yNRTB4cyRM3TCT1SDE/XryZhvb2kAIIiYJb37Nmqn3nTvjqnzo1iFI2eGrlfp7/6iDfPiuJH543rGMXKc+FpffCgEkw6Ta3xuduriSLGcAmEdkjIltFZJuIbLU7sN5q9th+PHj5aD7ZmcOD723HdOSDPigCbn4LxlwNy38Dy34FDof7g1Wql3rmi/38ddke5kwYwG8uG92+EdqNjLGWTq6tsBY684KZZVvjSnSX2B6FOsG3z0omt6yGp1fuJ9DPl99cPqr9f4x+gXDt8xDWH9Y+BUXpcM1CCAy1J2ileol/f7GfP320m8vH9edv88a3v4tso62vwZ4P4aJHrPFSXq7NZGGMSe+KQNSJ7rs4heq6Bp7/6iCB/j7cd3FK+xOGjw9c/Aer0fvjX8Dzs+HGV6FvO+bUV0od89TKNP7y8R6uGD+Ax68bj5+vK5UzzSjNtKbrGTgNpv/QvUHapIN3quwmIjx4+WhuOmMQT6/cz+PL93asSkoEzrgTbnoDitPh2fN18J5S7WSM4c8f7+YvH+/hqgmdTBQNdfDm7dbznKfAx9e9wdpEk4UXExEeuWos108ZyD9XpPHIB7s6ljDAWpbxu59a1VAvXGZNgawN30q1qb7BwS/e2srTK/fzrTMG8dh1EzqeKACWPwSHv4YrnoDooe4L1Gbe3aKi8PER/njNaYQE+vL8Vwcpr6njj9eMw7cj9aSxKXDH5/DO960icMZGuOIfENCBQURK9QJVtQ3c/epmPt2Vw48vGM5PLhzescbsRtvfhrVPwunfg3HXtX28F9Fk0Q34+FhVUuFB/jzx2T5Kq+r5xw0TCPLvQPE1uC/c8Aqs/juseBSytsC8/0L8GPcHrlQ3ll1SzR0vbmR7ZgmPXDWGW6Ynde6CubvhvQUw8AyY9Xu3xNiVtBqqmxAR7rloBA9ePpplO7O5YeFa8jsy0hushu+Z98Kt71rz5j97Pmx4TqullHLamlHMlf9azYG8cp6bP6XziaIsB16ZZ5Xi570Afu1c48IL2JYsROR5EcltXGHPuS1KRJaLyD7nc6Rzu4jIP0UkzTmWY1KTc+Y7j98nIvPtire7+M6MZJ6+aTK7s0u5+qmvSMst6/jFhpwL318Ng8+CD38Kr90MFe1cjEmpHua91KNc98wa/H19eOuuMzl/ZHznLlhTBi/PhYoC+NZr1tQ83ZCdJYsXgNknbfsl8JkxZjjwmfM1WGM5hjsfdwJPg5VcgIeAM4DTgYcaE0xvNntsP167czpVtQ6ufvLrjs1W2yg0Dm560+rrve8TeGoa7PnYfcEq1U3U1jt48L3t/HhxKuMS+vLegrMY2S+8cxetr4XXb4WcHXDdIkiY1PY5Xsq2ZGGMWQUUnrT5KmCR8+dFwJwm2180lrVAXxHpD1wMLDfGFBpjioDlnJqAeqXxA60/5qSYPnz3xY08tnxv+xdQauTjA2f9yGr8Do2HV6+H934IVcXuDVopL3W0uIrrnlnDi2vSuXPmEF6+4wxiQgM7d9GGenj3B7B/BVz5Txh+kXuC9ZCubrOIN8ZkATif45zbE4AjTY7LcG5rafspROROEdkoIhvz8vLcHrg3SugbzBvfn85c53oYty/aQGF7l2htqt9YuGMFzLgHUl+xShm7l7ovYKW80AdbM5n9j1Wk5Zbz9E2T+NWlo/DvTNdYsBLFO3fC9jfhwt/CxJvdE6wHeUsDd3N90Uwr20/daMxCY8wUY8yU2NhYtwbnzYL8ffnr3HE8MmcsX6UVcMkTq/h6fyfaHfwC4cKH4bufQUg0LL4R3rgNyrLdFLFS3qGipp5739jCglc2MzQ2lKU/OptLTuvf+Qs31MPbd1izyV74W5jxk85f0wt0dbLIcVYv4XzOdW7PAJrOQZEIZLayXTUhItwybTBv33UmfQL8uOk/6/jbsj3UNXRi8sCESVa11Hm/tkoX/5oK6xaCowNTpyvlZb7en8/F/1jFm99kcPf5w3jj+9MZFB3S+QvXVcOb34Ydb8NFv+sxiQK6PlksARp7NM3HWlipcfutzl5R04ASZzXVMmCWiEQ6G7ZnObepZoxNiOD9u2cwd1Ii//o8jTlPfsXu7NKOX9AvAM65F+5aA4lT4KN7YeG5Ol2I6rYqaur5zbvb+daz6/D39eH1703nZ7NSOl/tBFBVZK1UuWsJXPxHOOvHnb+mF5EOTx/R1oVFXgXOBWKAHKxeTe8CrwODgMPAPGNMoVhDIv+F1XhdCXzbGLPReZ3vAL9yXvZRY8x/23rvKVOmmI0bN7r3hrqZj7dn88A72yitruMnF47gzplDOvcfwhjY8Q588msoPQpj51rfnCKabUJSyut8siObh5fsIKu0mu+clczPZ6UQHOCmeZlKMuCla6FgvzXd+Glz3XPdLiYim4wxU5rdZ1ey8CRNFpaC8hoefG8HH27LYlT/cP5w9VgmDupkz+PaClj9D/jqCRAfOHOB9Q0qMMw9QSvlZkeLq3h4yQ6W78whJT6MP1wzlsmDo9z3BofXWd1j66rghpcgeab7rt3FNFn0ch9vt75R5ZRVc8u0wfxsVgoRwf6du2hROnz2O6u3R0gMnPtLmDS/W45MVT1TVW0DC1cd4Okv0gD4yYUjuH1GsnuqnMAqbW/6Lyy9DyISrWl04ke759oeoslCUVZdx98/2cuiNYeIDAng57NSuH7qwI5NSNjU0U3wyYPWGsJ9B8O591sTpHWTaZdVz2OM4YOtWfzpo90cLa7istP6c/+lI0mMdEMDdqPaSmsyzs3/g2EXwbXPevX62a7SZKGO2X60hN++v4MNh4oYMyCcX182mulDozt3UWMg7TNY8TtrYsKYFDjnPmtZV00aqgt9vT+fP320m60ZJYzqH85DV4xm2pBO/n2fLHs7vHU75O2Gs38O5/2qx/yda7JQJ2j85vXHpbvILKnm3JRYfjF7JKP6d3JqA2OsniCf/xHydkH0cGvCwrHXev36wqp7Sz1SzGPL97Jqbx4DIoL42awU5kxM6HzJuSmHA9YvhOUPWrM3X/0MDD3Pfdf3AposVLOq6xpY9PUhnvw8jbKaeq4cP4AfXTCcobGdXKfb4YDd78MXf4Gc7dB3EJz5I5hwEwS4sSpA9XpbM4r5x6f7WLE7l8gQf35w7lBunZ7Usen7W1OwH5bcDelfwYjZcNWT0CfGve/hBTRZqFaVVNbx9Bf7WfT1IWrqG7hqQgILzh/mnqSx92NY/ThkrLcawk+/A6bcDqG9Z5S9ci9jDGsPFPLUyjS+3JdP3xB/7jh7CPPPTCI00M0l2IZ6WPsUfP4o+AbCxY9aU3d0ZgEkL6bJQrkkv7yGZ1cd4MU16VTXNzB7TD9+cO5QxiX27dyFjYHDa6zutns/tv7TjZsHZ/zAmo9KKRfUNzj4ZGcOz355gM2Hi4kJDeA7M5K5ZdpgwoI62buvOYdWw9J7IXcnpFwGl/0dwt0wHYgX02Sh2iW/vIYXvjrEojWHKKuu58yh0XznrGTOHxmHT2frgPP3Wd/UUl+F+ioYNN0qbYy6Enxt+A+vur2Sqjre2HiEF74+REZRFYOiQrhj5hDmTU50f3UTQPER+PRhq1t4xCCrNDHqih5bmmhKk4XqkLLqOl5Zd5gXvj5EVkk1SdEhzD8ziWsmJXZ+nEZlIaS+DBv+A0WHoE8cTLwJJt0KUUPcEr/q3nZklvDS2nTe3ZxJVV0DpydFcfvZyVw4Kt69DdeNqorgy8dg3TPW67N+bM3A3Iva2TRZqE6pa3Dw8fZsnlt9kNQjxQT7+3Ll+AHcNG0QpyVEdG4Be4cD0j6FTS9YVVSmAZLOtuqFR13Zq/6jKusLyvtbsnhtw2G2ZJQQ5O/DVeMTuGX6YMYmRNjzptWlsOFZ+OqfUF0C42+A8x6AvgPbPreH0WSh3GZbRgkvr0vnvVTr296o/uFcNyWRORMSiOzTydHbpZmw+WVIfckqbQSEwZirYNwN1tKvPt4yo75ypwaH4au0fN7ZfJSPt2dTVddASnwY100dyNxJiUSE2FQ9WV1idYVd86RVqhg+Cy54EPqdZs/7dQOaLJTblVTVsWRLJm9sPMLWjBL8fYVzRsQxZ+IALhwV37m6ZIcDDn9tJY5dS6C2HMITrPEaY6+B/hN6Rf1xT2aMYUtGCe9vyeSDrZnklNYQHuTHZeMGcP3UgYxP7GSJtTUlGbD2adi0CGrLrK6w59wHCZPteb9uRJOFstWurFLe/iaDJVus//ShgX5cOCqOS0/rz8wRsZ1LHLWVsGcpbH0d9n8GjnqrTWPM1VY1Vf/xmji6CYfDkJpRzLId2SzdlsWRwqpjXzKumZTA+SPj7GmwBqtH3pF1sP5Z2Pmu9XrM1dZywv3H2/Oe3ZAmC9UlGhyGtQcKWJKaybKd2RRX1hEa6Me5KbHMGtOPc1NiCe9MF8fKQtj9AWx/Gw6usto3IgbBqMsh5VKrZ5WOFPcq1XUNrDlQwGe7cli+M4ec0hr8fIQzh8Vwxbj+zBrTr/OdJVoNoMRasW7D85CzDQIjrI4UZ3wfIgfb977dlCYL1eXqGhys2V/A0m1ZfLorh/zyWvx9halJUZyXEsd5I2MZGhva8aqGigLY+xHsXAIHPoeGWgjqa9U7D58Fwy6AEDdOQ61cdqSwkpV781i1N4+v0vKprG0g2N+Xs4fHcMlp/Tg/Jd6+dgiwVnM8tNrqbbdzidVFO34sTP2uNcllQB/73rub02ShPKrBYUg9UsQnO3NYuTuPPTllAAyICGLG8BhmDI/lrKHRRIcGduwNaspg/+ew5yPYtwwqC6y1NhImw9ALYOj51s9a6rBFUUUtaw4U8FVaPl/vL+BgfgUAiZHBnJsSywWj4pk+JNq+KiawqpUyv4Ftb1klifJsqxRx2lyrZ92AiVpd6QJNFsqrHC2uYuWeXL7cm8/X+/Mpra4HYER8KNOGRHNGcjRTkiKJDw9q/8UdDZCZCvs+sR6ZmwEDgeFWj6rkmZB8NsSN0d5VHZRVUsXGQ0WsP1jI+oOFx5J/nwBfzhgSzdnDY5g5IpYhMX3sa6QGayqOjPWw633rUXIEfAOskuVpc62Ga/9g+96/B9JkobxWfYODrUdLWHuggLUHCtl4qJDK2gYABkYFM3lQJBMG9mXCoEhG9Q8j0K+d304rC632jQOfW8+FB6ztQX1h8JlWO8eg6VYjpy7cdIrK2np2ZJay5Ugxm48Uszm9iMySasBKDpOTojg9KZLpQ2MYlxjhvoWFWlKRb5Ui9y2DfcuhutiaPmbo+TD6Ski5pEesK+EpmixUt1HX4GBnZikbDhWyKb2Ibw4XkVNaA4C/r5DSL4yxAyIYmxDBqP7hpPQLa9/kcSUZcPBLa/bQ9K+hcL+13TcQEiZB4hRImGJVW0Uk9qqqi4LyGvZkl7Ezq5QdmaXsyCwhLbcch/MjIqFvMJMGRzJpUF8mDYpkzIBw/OxODjVl1rKlh1ZZSSJ7q7U9JMYqQYyYBcMu1GV93USTherWskqqSD1czJaMErYfLWF7ZgnFlXXH9g+ODmF4XBgj4kMZHh/KkJhQhsT2cW1yubJsOLLe6lZ5ZJ21eFNDrbWvT6w1pmPABOu531hrNcBunECMMeSV1XAgv4K03PJjjz05ZeSV1Rw7Lj48kLEDIhgzIJzxA/syLrEvsWEdbFNqj9Is57/FemvyyawtVq83H38YeAYMPReGnG/9m/SQBYe8iSYL1aMYYzhaXMXurDJ2ZZWyO7uMvTllHMyvoN5x/O85NiyQwVEhDIoOYXBUHwZGBZMYGUJCZDDxYYHNfyuur7W6WGZssto7slKtFdGMw9ofGA7xYyBuNMSNsp5jU7xqbYPK2nqOFlWRUVzF0aIqjhRVcrigksOFlRzKr6DCWc0HEBLgy/C4UIbFhTGqfxgj+4Uzsn8YMR3tbOAqY6Asy1p1LnsLHN1s/b7LMq39fkEwYBIknWW1NQ08XXsxdQFNFqpXqGtwcCi/ggP5FRzIq+BAXjnphdYHZXZp9QnH+oiVTPpFWIkjLjyQ2NAgYsICiO4TQFSfQKL6+BMRHECEXx0BBbshe5u1mFP2dsjdBTUlxy8YHAUxIyB6GEQPtR6RyRCZBEGdW4HQGENFbQPFlbUUV9ZRWFFLYUUtBRW15JXVkFdWQ25ZNTml1WSXVB/rMNAowNeHxKhgBkWFkBTdhyGxfUiOsR4DIoI7P5Nw68FDeQ4UpFlJN2+P9bvL3Wn1WmsUPcxKDgmTIPF0a8oNbUPqcposVK9XXddAZnEVGUXWI7ukiqySarJLq8kttT5si5pUbZ0sJMCX0EA/woL8CA3yJ8TPhwTfIpLNYRLqjzCg7jBxtYeJrskgrC7/hHOr/PpSGtSf0sB+lAb2pyQgjmK/OIr9YinyjaLARFLp8KGqroHK2gaqahuoqG2gvKaO8up6yqrrTygxNeXvK8SGBhIbFkh8eBD9I4KIjwgioW+w9YgMJi4syJ5ZWsFKBtXFUHLUag8qOQLF6dbcXkXpVoeC2vLjxweEWSWxuFFWQogfa1XvBdk0SaBqF00WSrmgtt7h/MZec+zbe2lVHcWVdZRU1VFeU09ZjfXhXVVbf+yDvabeYT3qGqh3GAIdFSQ4shkkOQySXAZLDgmSf+wRLLWnvHcJYRT79KXMty/lfn2p8oukLqAvDUEROIIi8QvpS0CfvgT0iSQsoi/h4X3p27cvEWHhiDvGjzgc1uC1uirrw72mzHpUl1iT7FUVWT3LKvOtHkkVeVZ7T3kO1J9YasM30Bod3Xews5TlLG3FjLDm+OrGbT49XWvJQkcpKeUU4OdDv4gg+kV0YHxHM4wxGAMOY/ARQQQErA/e0kwoPXrsAzeiLJuIilxrZHpFNlTugPLi420lrfHxA79gq9rGN8BaREp8rQZg8cF6V2OVAozDml/LUW815NfXQkPNqR/4zREfCIm2Gv77xFgNzmH9rEd4AkQMhIgEa20SHcPS42iyUMom4kwQPpz0TTokynq0taSsw2G1i1QVWWsu1JRaz7UV1rf/2grrQ76uynpuqHU+6qzBiabBej4ekDOJ+FkPXz+rIdk3wBq85h9iPQJCrIb8oHCrS2pwlDV2ITBck0AvpslCKW/l42N9SOsgM+UF9GuCUkqpNmmyUEop1SZNFkoppdqkyUIppVSbNFkopZRqkyYLpZRSbdJkoZRSqk2aLJRSSrWpR84NJSJ5QHo7T4sB8ts8qufR++5d9L57l/be92BjTGxzO3pksugIEdnY0gRaPZned++i9927uPO+tRpKKaVUmzRZKKWUapMmi+MWejoAD9H77l30vnsXt923tlkopZRqk5YslFJKtUmThVJKqTb1qmQhIrNFZI+IpInIL5vZHygirzn3rxORpK6P0v1cuO+fishOEdkqIp+JyGBPxGmHtu69yXFzRcSISI/oXunKfYvIdc5/9x0i8kpXx2gHF/7WB4nI5yKy2fn3fqkn4nQnEXleRHJFZHsL+0VE/un8nWwVkUkdeiNrneCe/wB8gf3AECAA2AKMPumYu4B/O3++AXjN03F30X2fB4Q4f/5BT7hvV+/deVwYsApYC0zxdNxd9G8+HNgMRDpfx3k67i6674XAD5w/jwYOeTpuN9z3TGASsL2F/ZcCH2Etxj4NWNeR9+lNJYvTgTRjzAFjTC2wGLjqpGOuAhY5f34TuEBETlpAudtp876NMZ8bYyqdL9cCiV0co11c+TcHeAT4C1DdlcHZyJX7vgN40hhTBGCMye3iGO3gyn0bINz5cwSQ2YXx2cIYswoobOWQq4AXjWUt0FdE+rf3fXpTskgAjjR5neHc1uwxxph6oASI7pLo7OPKfTd1O9a3kJ6gzXsXkYnAQGPMB10ZmM1c+TcfAYwQka9EZK2IzO6y6Ozjyn0/DNwsIhnAUuDurgnNo9r7GdAsP7eF4/2aKyGc3G/YlWO6G5fvSURuBqYA59gaUddp9d5FxAd4HLitqwLqIq78m/thVUWdi1WS/FJExhpjim2OzU6u3PeNwAvGmL+LyHTgf877dtgfnse45XOtN5UsMoCBTV4ncmoR9NgxIuKHVUxtrXjXHbhy34jIhcADwJXGmJouis1ubd17GDAWWCkih7Dqc5f0gEZuV//W3zPG1BljDgJ7sJJHd+bKr2xnbQAAAvxJREFUfd8OvA5gjFkDBGFNtteTufQZ0JbelCw2AMNFJFlEArAasJecdMwSYL7z57nACuNsIerG2rxvZ1XMM1iJoifUXTdq9d6NMSXGmBhjTJIxJgmrveZKY8xGz4TrNq78rb+L1bEBEYnBqpY60KVRup8r930YuABAREZhJYu8Lo2y6y0BbnX2ipoGlBhjstp7kV5TDWWMqReRBcAyrF4TzxtjdojI74CNxpglwHNYxdI0rBLFDZ6L2D1cvO+/AqHAG872/MPGmCs9FrSbuHjvPY6L970MmCUiO4EG4F5jTIHnou48F+/7Z8CzInIPVlXMbd39C6GIvIpVnRjjbIt5CPAHMMb8G6tt5lIgDagEvt2h9+nmvyellFJdoDdVQymllOogTRZKKaXapMlCKaVUmzRZKKWUapMmC6WUUm3SZKGUUqpNmiyUUkq1SZOFUjYSkdHOdTJ2tLA/QkQcItKTRs6rHkiThVL2apxnalML+ydhTfT2TdeEo1THaLJQyl6NyaKl+aYmO581WSivpslCKXu1VbLQZKG6BZ0bSimbiIgvUAoEAuFNViNsesxerKnBhzinClfKK2nJQin7jAZCgF0tJIpwYBhQpIlCeTtNFkrZx9XG7c1dE45SHafJQin7NCaL1Bb2X+p81vYK5fU0WShln8bG66qTd4hIBPAt58uWSh5KeQ1NFkrZwLmG+3jnyxtFJLDJvv7Aa0CCc9OeLg5PqXbTZKGUPcZire+cAZwOpIvIUhFZA/x/e3eI00AURQH0PhJsTT076CYIG2AD1Qi2Q2vYCivAYtEsoKYJ4lcM6NeknUwg5+gRz925fzLvfyZZJ/n+efa1qrbLjAnnERYwj9/vFW9JHpN8JblPcpfprveHJLtMR1Q3Sd4XmBHO5j8LmEFV7ZM8JXkeY7wsPQ9cSrOAeXRrPuBP0SzgyqrqNskh08vYaoxxXHgkuJhmAde3ybTi40NQ8F9oFgC0NAsAWsICgJawAKAlLABoCQsAWsICgJawAKB1AhqR10KkstBkAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "mus = arange(0.01,1,0.01)\n",
    "ll = zeros(mus.shape)\n",
    "\n",
    "for i,m in enumerate(mus):\n",
    "    # \"create\" a coin for each mu and see how likely the tosses would be\n",
    "    ll[i] = -Coin(m).loglikelihood(ftosses)\n",
    "plot(mus, ll, label=\"Fair Coin\")\n",
    "\n",
    "for i,m in enumerate(mus):\n",
    "    ll[i] = -Coin(m).loglikelihood(btosses)\n",
    "plot(mus, ll, label=\"Biased Coin\")\n",
    "    \n",
    "xlabel(\"$\\mu$\", size=22)\n",
    "ylabel(\"negative log likelihood\")\n",
    "legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From the graph it looks like the minimum $\\mu$ is close to the true one of the coin. And differentiating shows the maximum likelihood estimator is\n",
    "\n",
    "$\\mu_{ML} = \\frac{1}{N} \\sum_{n=1}^{N}x_n$\n",
    "\n",
    "i.e. the sample mean (which you probably knew already)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "Fair Coin: $\\mu_{ML} = 0.49$"
      ],
      "text/plain": [
       "<IPython.core.display.Latex object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "Biased Coin: $\\mu_{ML} = 0.593$"
      ],
      "text/plain": [
       "<IPython.core.display.Latex object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "muf = sum(ftosses)/float(nsamples)\n",
    "mub = sum(btosses)/float(nsamples)\n",
    "\n",
    "latex(\"Fair Coin: $\\mu_{ML} = %s$\" % muf)\n",
    "latex(\"Biased Coin: $\\mu_{ML} = %s$\" % mub)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now what if we want to estimate the *number* of heads? Call it $k$.\n",
    "\n",
    "For example, if we toss our biased coin 5 times, the probability of getting 3 heads *then* 2 tails is:\n",
    "\n",
    "$p(x_1=1)p(x_2=1)p(x_3=1)p(x_4=0)p(x_5=0) = 0.6 \\times 0.6 \\times 0.6 \\times 0.4 \\times 0.4$\n",
    "\n",
    "or $\\mu^k(1-\\mu)^{N-k}$\n",
    "\n",
    "But this is just one way to get 3 heads, there are many other sequences of 5 trials that would give us 3 heads. How many exactly? This is given by the [Binomial coefficient](https://en.wikipedia.org/wiki/Binomial_coefficient):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "@vectorize\n",
    "def C(N,k):\n",
    "    return math.factorial(N)//(math.factorial(N-k)*math.factorial(k))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "There is 1 way of getting 0 heads\n",
      "There are 5 ways of getting 1 head\n",
      "There are 10 ways of getting 2 heads\n",
      "There are 10 ways of getting 3 heads\n",
      "There are 5 ways of getting 4 heads\n",
      "There is 1 way of getting 5 heads\n"
     ]
    }
   ],
   "source": [
    "Ntoss = 5\n",
    "for k in range(Ntoss+1):\n",
    "    c =  C(Ntoss,k)\n",
    "    s = lambda n: \"s\" if (n-1) else \"\"\n",
    "    print(\"There\", \"are\" if (c-1) else \"is\", c, \"way\" + s(c), \"of getting\", k, \"head\" + s(k))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This leads us to the Binomial distribution\n",
    "\n",
    "$Bin(k|N,\\mu) = C(N,k)\\times\\mu^k(1-\\mu)^{N-k}$\n",
    "\n",
    "With mean $N\\mu$\n",
    "\n",
    "We can sample a Binomial distribution by tossing a coin N times and seeing how many heads we get. And if we repeat this step many times we can draw a histogram of how many times we get each number of heads."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 0, 'k (number of heads)')"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAEJCAYAAABv6GdPAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAUWElEQVR4nO3df/BldX3f8eeLH4IjNoh8IZvdxaW6/sBkXOwGaUgtolWETMEJtBCjaOls7GCrjZN2zWSMzpQpaat2bBKStSDYWBEjBhRqpSuE0gZwF5Yfy4puEGHDDrtGQdGUDPDuH/ezw2W5u9+73/v97t398HzM3Lnnfs7nnPM+y/L6nv18z/ncVBWSpL4cMO0CJEnzz3CXpA4Z7pLUIcNdkjpkuEtShwx3SerQrOGe5NAktyW5M8nGJB9r7Zcl+W6SDe21orUnyaeSbE5yV5LXL/RJSJKe7aAx+jwBnFJVjyc5GLg5yf9o636rqv50p/5vB5a31xuAi9u7JGkvmTXca/CU0+Pt48Httbsnn84APtu2uyXJ4UkWVdXWXW1w5JFH1rJly8avWpLE+vXrv19VM6PWjXPlTpIDgfXAK4A/qKpbk/wL4MIkHwHWAqur6glgMfDQ0OZbWtsuw33ZsmWsW7durJORJA0k+d6u1o31C9WqeqqqVgBLgBOS/DzwYeDVwC8CRwD/dsfxRu1iRFGrkqxLsm779u3jlCFJGtMe3S1TVY8CNwKnVtXWGngC+AxwQuu2BVg6tNkS4OER+1pTVSurauXMzMh/VUiS5micu2Vmkhzell8IvAX4VpJFrS3AmcA9bZNrgHe3u2ZOBB7b3Xi7JGn+jTPmvgi4vI27HwBcWVVfTfKNJDMMhmE2AO9r/a8DTgM2Az8F3jv/ZUuSdmecu2XuAo4f0X7KLvoXcMHkpUmS5sonVCWpQ4a7JHXIcJekDhnuktShsZ5QlfRsy1ZfO+0SxvLARadPuwRNiVfuktQhw12SOmS4S1KHDHdJ6pDhLkkdMtwlqUOGuyR1yHCXpA4Z7pLUIcNdkjpkuEtShwx3SeqQ4S5JHTLcJalDhrskdchwl6QOGe6S1KFZwz3JoUluS3Jnko1JPtbaj01ya5LvJPlCkhe09kPa581t/bKFPQVJ0s7GuXJ/Ajilql4HrABOTXIi8HvAJ6tqOfBD4PzW/3zgh1X1CuCTrZ8kaS+aNdxr4PH28eD2KuAU4E9b++XAmW35jPaZtv7NSTJvFUuSZjXWmHuSA5NsALYB1wN/CTxaVU+2LluAxW15MfAQQFv/GPDS+SxakrR7Y4V7VT1VVSuAJcAJwGtGdWvvo67Sa+eGJKuSrEuybvv27ePWK0kawx7dLVNVjwI3AicChyc5qK1aAjzclrcASwHa+p8BfjBiX2uqamVVrZyZmZlb9ZKkkca5W2YmyeFt+YXAW4BNwA3AWa3becDVbfma9pm2/htV9Zwrd0nSwjlo9i4sAi5PciCDHwZXVtVXk9wLXJHk3wF3AJe0/pcA/y3JZgZX7OcsQN2SpN2YNdyr6i7g+BHt9zMYf9+5/f8BZ89LdZKkOfEJVUnqkOEuSR0y3CWpQ4a7JHXIcJekDhnuktQhw12SOmS4S1KHDHdJ6pDhLkkdMtwlqUOGuyR1yHCXpA4Z7pLUIcNdkjpkuEtShwx3SeqQ4S5JHTLcJalDhrskdWjWL8iW1L9lq6+ddgljeeCi06ddwn7DK3dJ6tCs4Z5kaZIbkmxKsjHJB1r7R5P8VZIN7XXa0DYfTrI5yX1J3raQJyBJeq5xhmWeBD5UVbcneTGwPsn1bd0nq+o/DXdOchxwDvBa4OeA/5XklVX11HwWLknatVmv3Ktqa1Xd3pZ/DGwCFu9mkzOAK6rqiar6LrAZOGE+ipUkjWePxtyTLAOOB25tTe9PcleSS5O8pLUtBh4a2mwLu/9hIEmaZ2OHe5LDgC8BH6yqHwEXAy8HVgBbgY/v6Dpi8xqxv1VJ1iVZt3379j0uXJK0a2OFe5KDGQT756rqKoCqeqSqnqqqp4FP88zQyxZg6dDmS4CHd95nVa2pqpVVtXJmZmaSc5Ak7WScu2UCXAJsqqpPDLUvGur2DuCetnwNcE6SQ5IcCywHbpu/kiVJsxnnbpmTgHcBdyfZ0Np+Gzg3yQoGQy4PAL8BUFUbk1wJ3MvgTpsLvFNGkvauWcO9qm5m9Dj6dbvZ5kLgwgnqkiRNwCdUJalDhrskdchwl6QOGe6S1CHDXZI6ZLhLUocMd0nqkOEuSR3ya/a04PwKN2nv88pdkjpkuEtShwx3SeqQ4S5JHTLcJalDhrskdchwl6QOGe6S1CHDXZI6ZLhLUocMd0nqkOEuSR0y3CWpQ7OGe5KlSW5IsinJxiQfaO1HJLk+yXfa+0tae5J8KsnmJHclef1Cn4Qk6dnGuXJ/EvhQVb0GOBG4IMlxwGpgbVUtB9a2zwBvB5a31yrg4nmvWpK0W7OGe1Vtrarb2/KPgU3AYuAM4PLW7XLgzLZ8BvDZGrgFODzJonmvXJK0S3s05p5kGXA8cCtwdFVthcEPAOCo1m0x8NDQZltamyRpLxk73JMcBnwJ+GBV/Wh3XUe01Yj9rUqyLsm67du3j1uGJGkMY4V7koMZBPvnquqq1vzIjuGW9r6ttW8Blg5tvgR4eOd9VtWaqlpZVStnZmbmWr8kaYRx7pYJcAmwqao+MbTqGuC8tnwecPVQ+7vbXTMnAo/tGL6RJO0d43xB9knAu4C7k2xobb8NXARcmeR84EHg7LbuOuA0YDPwU+C981qxJGlWs4Z7Vd3M6HF0gDeP6F/ABRPWJUmagE+oSlKHDHdJ6pDhLkkdMtwlqUOGuyR1yHCXpA4Z7pLUIcNdkjpkuEtShwx3SeqQ4S5JHTLcJalDhrskdchwl6QOGe6S1CHDXZI6ZLhLUocMd0nqkOEuSR0y3CWpQ4a7JHVo1nBPcmmSbUnuGWr7aJK/SrKhvU4bWvfhJJuT3JfkbQtVuCRp18a5cr8MOHVE+yerakV7XQeQ5DjgHOC1bZs/THLgfBUrSRrPrOFeVTcBPxhzf2cAV1TVE1X1XWAzcMIE9UmS5mCSMff3J7mrDdu8pLUtBh4a6rOltUmS9qK5hvvFwMuBFcBW4OOtPSP61qgdJFmVZF2Sddu3b59jGZKkUeYU7lX1SFU9VVVPA5/mmaGXLcDSoa5LgId3sY81VbWyqlbOzMzMpQxJ0i7MKdyTLBr6+A5gx5001wDnJDkkybHAcuC2yUqUJO2pg2brkOTzwMnAkUm2AL8LnJxkBYMhlweA3wCoqo1JrgTuBZ4ELqiqpxamdEnSrswa7lV17ojmS3bT/0LgwkmKkiRNxidUJalDhrskdchwl6QOGe6S1CHDXZI6ZLhLUocMd0nqkOEuSR0y3CWpQ4a7JHXIcJekDhnuktQhw12SOmS4S1KHDHdJ6pDhLkkdMtwlqUOGuyR1yHCXpA4Z7pLUIcNdkjpkuEtSh2YN9ySXJtmW5J6htiOSXJ/kO+39Ja09ST6VZHOSu5K8fiGLlySNNs6V+2XAqTu1rQbWVtVyYG37DPB2YHl7rQIunp8yJUl7YtZwr6qbgB/s1HwGcHlbvhw4c6j9szVwC3B4kkXzVawkaTxzHXM/uqq2ArT3o1r7YuChoX5bWpskaS+a71+oZkRbjeyYrEqyLsm67du3z3MZkvT8Ntdwf2THcEt739batwBLh/otAR4etYOqWlNVK6tq5czMzBzLkCSNMtdwvwY4ry2fB1w91P7udtfMicBjO4ZvJEl7z0GzdUjyeeBk4MgkW4DfBS4CrkxyPvAgcHbrfh1wGrAZ+Cnw3gWoWZI0i1nDvarO3cWqN4/oW8AFkxYlSZqMT6hKUocMd0nqkOEuSR0y3CWpQ4a7JHXIcJekDs16K6Qk7Y+Wrb522iWM5YGLTl+Q/Rru+6Dn+19KSZNzWEaSOmS4S1KHDHdJ6pDhLkkdMtwlqUOGuyR1yHCXpA4Z7pLUIcNdkjpkuEtShwx3SeqQ4S5JHTLcJalDhrskdWiiKX+TPAD8GHgKeLKqViY5AvgCsAx4APgnVfXDycqUJO2J+bhyf1NVraiqle3zamBtVS0H1rbPkqS9aCGGZc4ALm/LlwNnLsAxJEm7MWm4F/D1JOuTrGptR1fVVoD2ftSEx5Ak7aFJv2bvpKp6OMlRwPVJvjXuhu2HwSqAY445ZsIyJEnDJrpyr6qH2/s24MvACcAjSRYBtPdtu9h2TVWtrKqVMzMzk5QhSdrJnMM9yYuSvHjHMvBW4B7gGuC81u084OpJi5Qk7ZlJhmWOBr6cZMd+/ntVfS3JN4Erk5wPPAicPXmZkqQ9Medwr6r7gdeNaP9r4M2TFCVJmoxPqEpShwx3SeqQ4S5JHTLcJalDhrskdchwl6QOGe6S1CHDXZI6ZLhLUocMd0nqkOEuSR0y3CWpQ4a7JHXIcJekDhnuktQhw12SOmS4S1KHDHdJ6tAk36G6T1i2+tpplzCWBy46fdolSHoe8cpdkjpkuEtShwx3SerQgoV7klOT3Jdkc5LVC3UcSdJzLUi4JzkQ+APg7cBxwLlJjluIY0mSnmuhrtxPADZX1f1V9bfAFcAZC3QsSdJOFircFwMPDX3e0tokSXtBqmr+d5qcDbytqv55+/wu4ISq+pdDfVYBq9rHVwH3zXshc3ck8P1pFzHPejun3s4H+jun3s4H9r1zellVzYxasVAPMW0Blg59XgI8PNyhqtYAaxbo+BNJsq6qVk67jvnU2zn1dj7Q3zn1dj6wf53TQg3LfBNYnuTYJC8AzgGuWaBjSZJ2siBX7lX1ZJL3A/8TOBC4tKo2LsSxJEnPtWBzy1TVdcB1C7X/BbZPDhdNqLdz6u18oL9z6u18YD86pwX5haokabqcfkCSOmS476S3aROSXJpkW5J7pl3LfEiyNMkNSTYl2ZjkA9OuaRJJDk1yW5I72/l8bNo1zZckBya5I8lXp13LpJI8kOTuJBuSrJt2PeNwWGZImzbh28A/YnA75zeBc6vq3qkWNoEkbwQeBz5bVT8/7XomlWQRsKiqbk/yYmA9cOb++t8oSYAXVdXjSQ4GbgY+UFW3TLm0iSX5TWAl8Heq6lemXc8kkjwArKyqfeke993yyv3Zups2oapuAn4w7TrmS1Vtrarb2/KPgU3sx08/18Dj7ePB7bXfX3ElWQKcDvzXadfyfGW4P5vTJuxHkiwDjgdunW4lk2nDFxuAbcD1VbVfn0/zn4F/Azw97ULmSQFfT7K+PV2/zzPcny0j2vb7q6geJTkM+BLwwar60bTrmURVPVVVKxg8yX1Ckv16+CzJrwDbqmr9tGuZRydV1esZzHR7QRvu3KcZ7s8267QJmr42Nv0l4HNVddW065kvVfUocCNw6pRLmdRJwD9u49RXAKck+ZPpljSZqnq4vW8DvsxgCHefZrg/m9Mm7OPaLyAvATZV1SemXc+kkswkObwtvxB4C/Ct6VY1mar6cFUtqaplDP4f+kZV/fqUy5qzJC9qv7wnyYuAtwL7/N1nhvuQqnoS2DFtwibgyv192oQknwf+AnhVki1Jzp92TRM6CXgXg6vBDe112rSLmsAi4IYkdzG4uLi+qvb7Wwc7czRwc5I7gduAa6vqa1OuaVbeCilJHfLKXZI6ZLhLUocMd0nqkOEuSR0y3CWpQ4a7FkSSZePMRJlk0d6YNTDJ47P3mpfjfD7JXUn+9U7tlyU5awGO954kv7+b9b+Q5LL5Pq72fQv2TUzSmH4T+PS0i9idJAe1ZyBm6/ezwC9V1cv2Qlljqaq7kyxJckxVPTjterT3eOWuBZfk77Z5vX9xxOpfBb7W+r0nyVVJvpbkO0n+w9A+Hh9aPmvH1Wi7Ir64zfF+f5J/2Oaw37TzFWuSjye5PcnaJDOt7eXteOuT/O8krx7a7yeS3AD83k77OTTJZ9r83nckeVNb9XXgqPZg1T8Yca5vTPJ/W51nDe3vt5J8s13xf2yo/c9aXRuHJ6tK8t4k307y5wwe6trRfnaSe9rc8DcNHfcrDJ4U1fNJVfnyNe8vYBmDR7RfBdwBrBjR51hg/dDn9wD3Az8DHAp8D1ja1j0+1O8s4LK2fBmD+UvCYHrmHwG/wODCZf2O4zKYAO6dbfkjwO+35bXA8rb8BgaPyu/Y71eBA0fU/SHgM2351cCDrd5lwD27+PO4DPhiq+s4BlNLw+BR9jWt/gPaMd/Y1h3R3l/Y/ixfyuCJ1geBGeAFwP8ZOpe7gcVt+fChY58EfGXafyd87d2XwzJaSDPA1cCv1uhpHBYB23dqW1tVjwEkuRd4Gc+ehnmUr1RVJbkbeKSq7m7bb2QQuBsYTD37hdb/T4Cr2sySvwR8cTBlDQCHDO33i1X11Ijj/TLwXwCq6ltJvge8ksEPlt35s6p6Grg3ydGt7a3tdUf7fBiwHLgJ+FdJ3tHal7b2nwVurKrt7Ry/0I4Ng6C/LMmVwPCEatuAn5ulNnXGcNdCeoxBMJ8EjAr3v2FwxTvsiaHlp3jm7+jwPBm72ubpnbZ/ml3/HS8GV8qP1mC63VF+sov2UVNDj2O4tgy9//uq+uNnHSA5mcEkYn+/qn6a5EaeOe+Rc4ZU1fuSvIHBl2RsSLKiqv66bfc3c6xZ+ynH3LWQ/hY4E3h3kl8bsf7bDK6sx/FIktckOQB4x6y9n+sABsM5AL8G3FyDeeC/m+RsGMw4meR1Y+zrJuCdbZtXAscA982hJhhMUvfP2r8iSLI4yVEMhqZ+2IL91cCJrf+twMlJXtqmPj57x46SvLyqbq2qjwDf55npq1/JfjCLoeaXV+5aUFX1kwy+vOH6JD+pqqt3WveXSV5RVZtn2dVqBuPRDzEIqsP2sJSfAK9Nsp7Bvyj+aWt/J3Bxkt9h8BV3VwB3zrKvPwT+qA0DPQm8p6qeGBraGVtVfT3Ja4C/aNs/Dvw6g18yv6/NFnkfcEvrvzXJRxnM9LkVuB04sO3uPyZZzuBfA2uHzuNNwLV7XJz2a84KqalqY8p/r6p+Z9q19CjJIcCfA79cY9zOqX545a6pqqovJ3nptOvo2DHAaoP9+ccrd0nqkL9QlaQOGe6S1CHDXZI6ZLhLUocMd0nqkOEuSR36//OACFoNdqh+AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "results = []\n",
    "\n",
    "N_binom = 5\n",
    "mu_binom = 0.6\n",
    "Ntrials = 1000\n",
    "\n",
    "for x in range(Ntrials):\n",
    "    nheads = sum(Coin(mu_binom).toss(N_binom))\n",
    "    results.append(nheads)\n",
    "      \n",
    "hist(results, bins=range(N_binom+2), rwidth=.8, align='left')\n",
    "\n",
    "xlabel(\"k (number of heads)\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compare with the actual distribution:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 0, 'k')"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAEGCAYAAABrQF4qAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAR/UlEQVR4nO3df6xfdX3H8edrRXDB6VBuFtMWW7Ua61xgudY/iCxRwDoMJRnGYlxwIWlcaOZCllmjwazGBDVx/lMjzWzi3FhF2eLNqOuIogtx6L1A1bWs41KZ3NWEq2U61IHF9/64x+Xr9Vvuae+3/fZ+eD6Sb+75/Dr3fQJ93ZPzPd/zTVUhSWrXr427AEnS6WXQS1LjDHpJapxBL0mNM+glqXHnjLuAxS688MJat27duMuQpBXlvvvu+35VTQwbO+uCft26dczMzIy7DElaUZL854nGvHQjSY0z6CWpcQa9JDWuV9An2ZzkcJLZJDuGjL8rybeTHEhyT5KNXf+6JD/t+g8k+eSoD0CS9MyWfDM2ySpgF3AFMAdMJ5mqqkMD026rqk92868GPgZs7sYerqqLR1u2JKmvPmf0m4DZqjpSVU8Be4EtgxOq6kcDzfMBn5QmSWeJPkG/Gnh0oD3X9f2SJDcmeRj4CPAnA0PrkzyQ5KtJXj/sFyTZlmQmycz8/PxJlC9JWkqfoM+Qvl85Y6+qXVX1MuA9wPu77u8BF1XVJcBNwG1Jnj9k7e6qmqyqyYmJoff7S5JOUZ+gnwPWDrTXAEefYf5e4BqAqnqyqn7Qbd8HPAy84tRKlSSdij6fjJ0GNiRZD/wXsBV4++CEJBuq6qGueRXwUNc/ARyrqqeTvBTYABwZVfHS2WDdjjvHXUIvj9xy1bhL0JgsGfRVdTzJdmA/sArYU1UHk+wEZqpqCtie5HLgZ8DjwPXd8suAnUmOA08D76qqY6fjQCRJw/V61k1V7QP2Leq7eWD73SdYdwdwx3IKlCQtj5+MlaTGGfSS1DiDXpIaZ9BLUuMMeklqnEEvSY0z6CWpcQa9JDXOoJekxhn0ktQ4g16SGmfQS1LjDHpJapxBL0mNM+glqXEGvSQ1zqCXpMYZ9JLUOINekhpn0EtS4wx6SWpcr6BPsjnJ4SSzSXYMGX9Xkm8nOZDkniQbB8be2607nORNoyxekrS0JYM+ySpgF/BmYCNw3WCQd26rqtdU1cXAR4CPdWs3AluBVwObgU90+5MknSF9zug3AbNVdaSqngL2AlsGJ1TVjwaa5wPVbW8B9lbVk1X1HWC2258k6Qw5p8ec1cCjA+054HWLJyW5EbgJOBd4w8DaexetXT1k7TZgG8BFF13Up25JUk99zugzpK9+paNqV1W9DHgP8P6TXLu7qiaranJiYqJHSZKkvvoE/RywdqC9Bjj6DPP3Atec4lpJ0oj1CfppYEOS9UnOZeHN1anBCUk2DDSvAh7qtqeArUnOS7Ie2AB8Y/llS5L6WvIafVUdT7Id2A+sAvZU1cEkO4GZqpoCtie5HPgZ8Dhwfbf2YJLbgUPAceDGqnr6NB2LJGmIPm/GUlX7gH2L+m4e2H73M6z9EPChUy1QkrQ8fjJWkhpn0EtS4wx6SWqcQS9JjTPoJalxBr0kNc6gl6TGGfSS1DiDXpIaZ9BLUuMMeklqnEEvSY0z6CWpcQa9JDXOoJekxhn0ktQ4g16SGmfQS1Ljen2VoKRnl3U77hx3Cb08cstV4y5hRfCMXpIaZ9BLUuN6BX2SzUkOJ5lNsmPI+E1JDiX5VpIvJXnJwNjTSQ50r6lRFi9JWtqS1+iTrAJ2AVcAc8B0kqmqOjQw7QFgsqp+kuSPgY8Ab+vGflpVF4+4bklST33O6DcBs1V1pKqeAvYCWwYnVNXdVfWTrnkvsGa0ZUqSTlWfoF8NPDrQnuv6TuQG4IsD7ecmmUlyb5Jrhi1Isq2bMzM/P9+jJElSX31ur8yQvho6MXkHMAn83kD3RVV1NMlLgS8n+XZVPfxLO6vaDewGmJycHLpvSdKp6XNGPwesHWivAY4unpTkcuB9wNVV9eQv+qvqaPfzCPAV4JJl1CtJOkl9gn4a2JBkfZJzga3AL909k+QS4FYWQv6xgf4LkpzXbV8IXAoMvokrSTrNlrx0U1XHk2wH9gOrgD1VdTDJTmCmqqaAjwLPAz6XBOC7VXU18Crg1iQ/Z+GPyi2L7taRJJ1mvR6BUFX7gH2L+m4e2L78BOu+BrxmOQVKkpbHT8ZKUuMMeklqnEEvSY0z6CWpcQa9JDXOoJekxvkNUzrj/PYi6czyjF6SGmfQS1LjDHpJapxBL0mNM+glqXEGvSQ1zqCXpMYZ9JLUOINekhpn0EtS4wx6SWqcQS9JjTPoJalxBr0kNa5X0CfZnORwktkkO4aM35TkUJJvJflSkpcMjF2f5KHudf0oi5ckLW3JoE+yCtgFvBnYCFyXZOOiaQ8Ak1X1O8DngY90a18IfAB4HbAJ+ECSC0ZXviRpKX3O6DcBs1V1pKqeAvYCWwYnVNXdVfWTrnkvsKbbfhNwV1Udq6rHgbuAzaMpXZLUR5+gXw08OtCe6/pO5AbgiyezNsm2JDNJZubn53uUJEnqq0/QZ0hfDZ2YvAOYBD56MmurandVTVbV5MTERI+SJEl99Qn6OWDtQHsNcHTxpCSXA+8Drq6qJ09mrSTp9OkT9NPAhiTrk5wLbAWmBickuQS4lYWQf2xgaD9wZZILujdhr+z6JElnyDlLTaiq40m2sxDQq4A9VXUwyU5gpqqmWLhU8zzgc0kAvltVV1fVsSQfZOGPBcDOqjp2Wo5EkjTUkkEPUFX7gH2L+m4e2L78GdbuAfacaoGSpOXxk7GS1DiDXpIaZ9BLUuMMeklqnEEvSY0z6CWpcQa9JDXOoJekxhn0ktQ4g16SGmfQS1LjDHpJapxBL0mNM+glqXEGvSQ1zqCXpMYZ9JLUOINekhpn0EtS4wx6SWqcQS9JjesV9Ek2JzmcZDbJjiHjlyW5P8nxJNcuGns6yYHuNTWqwiVJ/Zyz1IQkq4BdwBXAHDCdZKqqDg1M+y7wTuDPhuzip1V18QhqlSSdgiWDHtgEzFbVEYAke4EtwP8HfVU90o39/DTUKElahj6XblYDjw6057q+vp6bZCbJvUmuGTYhybZuzsz8/PxJ7FqStJQ+QZ8hfXUSv+OiqpoE3g58PMnLfmVnVburarKqJicmJk5i15KkpfQJ+jlg7UB7DXC07y+oqqPdzyPAV4BLTqI+SdIy9Qn6aWBDkvVJzgW2Ar3unklyQZLzuu0LgUsZuLYvSTr9lgz6qjoObAf2Aw8Ct1fVwSQ7k1wNkOS1SeaAtwK3JjnYLX8VMJPkm8DdwC2L7taRJJ1mfe66oar2AfsW9d08sD3NwiWdxeu+BrxmmTVKkpbBT8ZKUuMMeklqnEEvSY0z6CWpcQa9JDXOoJekxhn0ktQ4g16SGmfQS1LjDHpJapxBL0mNM+glqXEGvSQ1zqCXpMYZ9JLUOINekhpn0EtS4wx6SWqcQS9JjTPoJalxBr0kNa5X0CfZnORwktkkO4aMX5bk/iTHk1y7aOz6JA91r+tHVbgkqZ8lgz7JKmAX8GZgI3Bdko2Lpn0XeCdw26K1LwQ+ALwO2AR8IMkFyy9bktRXnzP6TcBsVR2pqqeAvcCWwQlV9UhVfQv4+aK1bwLuqqpjVfU4cBeweQR1S5J66hP0q4FHB9pzXV8fvdYm2ZZkJsnM/Px8z11LkvroE/QZ0lc9999rbVXtrqrJqpqcmJjouWtJUh99gn4OWDvQXgMc7bn/5ayVJI1An6CfBjYkWZ/kXGArMNVz//uBK5Nc0L0Je2XXJ0k6Q5YM+qo6DmxnIaAfBG6vqoNJdia5GiDJa5PMAW8Fbk1ysFt7DPggC38spoGdXZ8k6Qw5p8+kqtoH7FvUd/PA9jQLl2WGrd0D7FlGjZKkZfCTsZLUOINekhpn0EtS4wx6SWqcQS9JjTPoJalxBr0kNa7XffSStJKt23HnuEvo5ZFbrjot+zXoz3LP9v9BJS2fl24kqXEGvSQ1zqCXpMYZ9JLUOINekhpn0EtS4wx6SWqcQS9JjTPoJalxBr0kNc6gl6TGGfSS1LheQZ9kc5LDSWaT7Bgyfl6Sz3bjX0+yrutfl+SnSQ50r0+OtnxJ0lKWfHplklXALuAKYA6YTjJVVYcGpt0APF5VL0+yFfgw8LZu7OGqunjEdUuSeupzRr8JmK2qI1X1FLAX2LJozhbg093254E3JsnoypQknao+Qb8aeHSgPdf1DZ1TVceBHwIv6sbWJ3kgyVeTvH7YL0iyLclMkpn5+fmTOgBJ0jPrE/TDzsyr55zvARdV1SXATcBtSZ7/KxOrdlfVZFVNTkxM9ChJktRXn6CfA9YOtNcAR080J8k5wAuAY1X1ZFX9AKCq7gMeBl6x3KIlSf31CfppYEOS9UnOBbYCU4vmTAHXd9vXAl+uqkoy0b2ZS5KXAhuAI6MpXZLUx5J33VTV8STbgf3AKmBPVR1MshOYqaop4FPAZ5LMAsdY+GMAcBmwM8lx4GngXVV17HQciCRpuF5fDl5V+4B9i/puHtj+X+CtQ9bdAdyxzBolScvgJ2MlqXEGvSQ1zqCXpMYZ9JLUOINekhpn0EtS4wx6SWqcQS9JjTPoJalxBr0kNc6gl6TGGfSS1DiDXpIaZ9BLUuMMeklqnEEvSY0z6CWpcQa9JDWu11cJriTrdtw57hJ6eeSWq8ZdgqRnCc/oJalxBr0kNa5X0CfZnORwktkkO4aMn5fks93415OsGxh7b9d/OMmbRle6JKmPJYM+ySpgF/BmYCNwXZKNi6bdADxeVS8H/hL4cLd2I7AVeDWwGfhEtz9J0hnS54x+EzBbVUeq6ilgL7Bl0ZwtwKe77c8Db0ySrn9vVT1ZVd8BZrv9SZLOkD533awGHh1ozwGvO9Gcqjqe5IfAi7r+exetXb34FyTZBmzrmk8kOdyr+jPnQuD7o9xhPjzKvZ201o4H2jum1o4H2jums+14XnKigT5BnyF91XNOn7VU1W5gd49axiLJTFVNjruOUWnteKC9Y2rteKC9Y1pJx9Pn0s0csHagvQY4eqI5Sc4BXgAc67lWknQa9Qn6aWBDkvVJzmXhzdWpRXOmgOu77WuBL1dVdf1bu7ty1gMbgG+MpnRJUh9LXrrprrlvB/YDq4A9VXUwyU5gpqqmgE8Bn0kyy8KZ/NZu7cEktwOHgOPAjVX19Gk6ltPprL2sdIpaOx5o75haOx5o75hWzPFk4cRbktQqPxkrSY0z6CWpcQb9M1jq0Q8rTZI9SR5L8m/jrmUUkqxNcneSB5McTPLucde0XEmem+QbSb7ZHdNfjLumUUiyKskDSf5x3LWMQpJHknw7yYEkM+OuZyleoz+B7lEN/wFcwcJtotPAdVV1aKyFLUOSy4AngL+uqt8edz3LleTFwIur6v4kvwHcB1yzwv8bBTi/qp5I8hzgHuDdVXXvEkvPakluAiaB51fVW8Zdz3IleQSYrKqRfmDqdPGM/sT6PPphRamqf2HhrqgmVNX3qur+bvt/gAcZ8snrlaQWPNE1n9O9VvTZWJI1wFXAX427lmcrg/7Ehj36YUWHSMu6J6ZeAnx9vJUsX3eZ4wDwGHBXVa30Y/o48OfAz8ddyAgV8M9J7use4XJWM+hPrNfjGzR+SZ4H3AH8aVX9aNz1LFdVPV1VF7PwSfJNSVbsZbYkbwEeq6r7xl3LiF1aVb/LwlN9b+wui561DPoT8/ENK0B3HfsO4G+r6u/HXc8oVdV/A19h4RHfK9WlwNXdNe29wBuS/M14S1q+qjra/XwM+AfO8qfyGvQn1ufRDxqj7o3LTwEPVtXHxl3PKCSZSPKb3favA5cD/z7eqk5dVb23qtZU1ToW/g19uareMeayliXJ+d2b/yQ5H7gSOKvvZDPoT6CqjgO/ePTDg8DtVXVwvFUtT5K/A/4VeGWSuSQ3jLumZboU+EMWzhIPdK/fH3dRy/Ri4O4k32LhZOOuqmrilsSG/BZwT5JvsvDsrjur6p/GXNMz8vZKSWqcZ/SS1DiDXpIaZ9BLUuMMeklqnEEvSY0z6KUekqxr5amfevYx6CWpcQa9dJKSvLR7tvprx12L1IdBL52EJK9k4dk6f1RV0+OuR+rjnHEXIK0gE8AXgD9Y6Y/D0LOLZ/RSfz9k4TsKLh13IdLJ8Ixe6u8p4Bpgf5Inquq2cRck9WHQSyehqn7cfZnGXUl+XFVfGHdN0lJ8eqUkNc5r9JLUOINekhpn0EtS4wx6SWqcQS9JjTPoJalxBr0kNe7/ADgOJlyiKrUIAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "ks = arange(0,N_binom+1,1)\n",
    "\n",
    "def binomialpmf(k,N,mu):\n",
    "    return C(N, k) * mu**k * (1-mu)**(N-k)\n",
    "\n",
    "bar(ks, binomialpmf(ks, N_binom, mu_binom), align='center')\n",
    "xlabel(\"k\")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}