{
 "metadata": {
  "name": "2.1-binary-variables"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%pylab inline\n",
      "from IPython.display import display, Math, Latex\n",
      "maths = lambda s: display(Math(s))\n",
      "latex = lambda s: display(Latex(s))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "Welcome to pylab, a matplotlib-based Python environment [backend: module://IPython.zmq.pylab.backend_inline].\n",
        "For more information, type 'help(pylab)'.\n"
       ]
      }
     ],
     "prompt_number": 61
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "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",
      "    "
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 62
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# 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",
      "\n",
      "print \"Fair Coin   -- \",\"Heads:\", fheads, \"Tails:\", ftails\n",
      "print \"Biased Coin -- \",\"Heads:\", bheads, \"Tails:\", btails\n",
      "\n",
      "#print \"Tosses:\", btosses[:100],\"...\"\n",
      "\n",
      "b = bar([0,1,2,3],[fheads,ftails,bheads,btails], color=['r','b'], width=.5)\n",
      "\n",
      "b[0].set_label(\"Heads\")\n",
      "b[1].set_label(\"Tails\")\n",
      "gca().axes.xaxis.set_ticklabels([])\n",
      "xlabel(\"Fair Coin/Biased Coin\")\n",
      "legend()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Fair Coin   --  Heads: 493 Tails: 507\n",
        "Biased Coin --  Heads: 602 Tails: 398\n"
       ]
      },
      {
       "output_type": "pyout",
       "prompt_number": 63,
       "text": [
        "<matplotlib.legend.Legend at 0x60ee750>"
       ]
      },
      {
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAXMAAAEDCAYAAADHmORTAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAHnpJREFUeJzt3X1QVXXix/H34cE2FRVXuNqFX9cEw+sT+HC1B43SizuV\nhFqM6Na1B22zbWsfcs2tSXdbwbImtaHtgTXWbUMbJ2Ad10FN2mx3FvOhdpc20IWCy4MhouLDEHB+\nf6h3IRFQQeT4ec2cmcs533PO93s4fO7he7/3HMM0TRMREenW/Lq6AiIicukU5iIiFqAwFxGxAIW5\niIgFKMxFRCxAYS4iYgFthvmXX35JTEyMb+rbty+rV6+muroat9vN0KFDiYuLo6amxrdOcnIykZGR\nREVFkZOT06kNEBERMC5knHljYyN2u528vDzWrFnDgAEDWLRoEStWrODw4cOkpKSQn5/PnDlz2LVr\nF16vl6lTp1JQUICfn/4JEBHpLBeUsNu2bSMiIoLw8HCys7PxeDwAeDweMjMzAcjKyiIpKYnAwEAc\nDgcRERHk5eV1fM1FRMTngsI8IyODpKQkACorK7HZbADYbDYqKysBKCsrIywszLdOWFgYXq+3o+or\nIiItCGhvwbq6Ov785z+zYsWKc5YZhoFhGOdd97vLWisrIiLnd76e8XZfmf/lL39h7NixhISEAKev\nxisqKgAoLy8nNDQUALvdTklJiW+90tJS7HZ7ixXqrtPzzz/f5XVQ/bu+Hldb3VX/rp9a0+4wf++9\n93xdLADx8fGkp6cDkJ6eTkJCgm9+RkYGdXV1FBUVUVhYiMvlau9uRETkIrSrm+X48eNs27aNt956\nyzdv8eLFJCYmkpaWhsPhYMOGDQA4nU4SExNxOp0EBASQmpqqbhURkU7WrjDv1asXVVVVzeb179+f\nbdu2tVh+yZIlLFmy5NJrd4WKjY3t6ipcEtW/63TnuoPqfyW7oHHmHbZTw2iz/0dERJprLTvbPZpF\nRORC9O/fn8OHD3d1Nbql4OBgqqurL2gdXZmLSKfQ3/nFO9+xa+2Y6jv2IiIWoDAXEbEAhbmIiAUo\nzEVEOklsbCxpaWmXZV8KcxG5LPr36eO7j1NnTP379GlXPRwOB9u3b28275133mHSpEkd3ua27lvV\nkTQ0UUQui8PHjtGZY1uMY8faV+4yBuzlpCtzEbnqNQ33srIyZs2aRWhoKDfccANr1qzxLcvLy+Om\nm24iODiY6667jieeeIJvv/3Wt3zr1q1ERUXRr18/nnjiiWY3yNq/fz+33XYb/fr1IyQkhNmzZ3do\nGxTmInLV+e5Y7bM/m6bJ9OnTiYmJoaysjO3bt/Pqq6/6Hn8ZEBDAqlWrOHToEH//+9/Zvn07qamp\nAFRVVTFr1iyWL1/OoUOHGDJkCJ988onvjeK5557jBz/4ATU1NXi9Xn7yk590aJsU5iJyVTFNk4SE\nBIKDg33T448/jmEY7Nq1i6qqKp599lkCAgIYPHgwjzzyCBkZGQCMGTMGl8uFn58f119/PQsWLOCj\njz4CYPPmzYwYMYKZM2fi7+/PU089xcCBA3377dGjB8XFxXi9Xnr06MHNN9/coe1SmIvIVcUwDLKy\nsjh8+LBvSk1NxTRNvvrqK8rKypoFfXJyMgcPHgSgoKCAu+++m0GDBtG3b19+9atfcejQIeDcp6wB\nhIeH+16/+OKLmKaJy+VixIgRrF27tkPbpTAXkave2W6W8PBwBg8e3Czojx49yqZNmwB47LHHcDqd\n7N+/nyNHjvDb3/6WxsZGAK677rpmD+YxTbPZzzabjTfffBOv18sbb7zBwoUL+e9//9thbVCYi4ic\n4XK5CAoK4sUXX+TkyZM0NDTwr3/9i08//RSA2tpagoKC6NmzJ//5z394/fXXfeveeeed/Pvf/+aD\nDz6gvr6e1atX+57GBvD+++9TWloKQL9+/TAMAz+/jotghbmIXBbBQUEY0GlTcFDQRdft7HBFPz8/\nNm3axL59+7jhhhsICQlhwYIFHD16FICVK1fypz/9iT59+rBgwQJmz57t+4BzwIABvP/++yxevJgB\nAwawf/9+br31Vt8+Pv30UyZOnEhQUBD33HMPq1evxuFwXHSdz2mD7pooIp1Bf+cXT3dNFBG5SinM\nRUQsQGEuImIBCnMREQtQmIuIWIDCXETEAtoV5jU1Ndx7770MGzYMp9PJP/7xD6qrq3G73QwdOpS4\nuDhqamp85ZOTk4mMjCQqKsp3gxoREek87QrzJ598kjvvvJMvvviCzz//nKioKFJSUnC73RQUFDBl\nyhRSUlIAyM/PZ/369eTn57NlyxYWLlzo+7qriIh0jjbD/MiRI3z88cc89NBDwOlbQPbt25fs7Gw8\nHg8AHo+HzMxMALKyskhKSiIwMBCHw0FERAR5eXmd2AQRkc4XFBREcXExAPPmzeO5557r2gp9R5tP\nGioqKiIkJIQHH3yQzz77jLFjx/Lqq69SWVmJzWYDTt9AprKyEjh957CJEyf61g8LC8Pr9Z6z3aVL\nl/pex8bGEhsbe4lNEbny9e/Th8PtfCLOlSg4KIjqM19tv1B9+vTn2LHDHVyj/wkKCubo0eo2y/Xu\n3dv3Ffzjx4/zve99D39/fwDefPNNkpKSWlzvWJPf2+V6WlFubi65ubntKttmmNfX17Nnzx5ee+01\nxo8fz1NPPeXrUjmrrYa1tKxpmItcLTr70Wmdrb2PZmvJ6SDvvNYfO9a+cK2trfW9Hjx4MGlpadxx\nxx0XvL/LcauC717oLlu27Lxl2+xmCQsLIywsjPHjxwNw7733smfPHgYOHOi7I1h5eTmhoaEA2O32\nZrd9LC0txW63X1RDREQul7YeCefn59fiLWurqqq4++67CQ4O5vvf/z6TJ0/uknvStBnmAwcOJDw8\nnIKCAgC2bdvG8OHDmT59Ounp6QCkp6eTkJAAQHx8PBkZGdTV1VFUVERhYSEul6sTmyAiculaeyRc\nS872OLz88suEh4dTVVXFwYMHSU5O7pIHRrfZzQKwZs0a5s6dS11dHUOGDGHt2rU0NDSQmJhIWloa\nDoeDDRs2AOB0OklMTMTpdBIQEEBqaqoln4QtItYyZswY3+umj4R78sknW12vR48elJeXU1xczJAh\nQ7jllls6u6otaleYjx49ml27dp0zf9u2bS2WX7JkCUuWLLm0momIXEYFBQX87Gc/Y/fu3Zw4cYL6\n+nrGjRt33vJnu1Kefvppli5dSlxcHAALFizgl7/85WWpc1P6BqiICK0/Eq41vXv3ZuXKlRw4cIDs\n7GxeeeUVPvzww8tQ4+YU5iIitP5IuO9q+gHnpk2b2L9/P6Zp0qdPH/z9/X1DHS8nhbmIXBZBQcF0\n3kPjjDPbv3itPRIOOOf12Z/379+P2+0mKCiIm2++mccff5zbbrvtkupyMfTYOJHLyDCM7j3OnPaP\nr9bf+cXTY+NERK5SCnMREQtQmIuIWEC7xpmLiFyo4OBgfWHwIgUHX/iHufoAVOQyupo+AJWOpw9A\nRUQsTmEuImIBCnMREQtQmIuIWIDCXETEAhTmIiIWoDAXEbEAhbmIiAUozEVELEBhLiJiAQpzEREL\nUJiLiFiAwlxExAIU5iIiFtCuMHc4HIwaNYqYmBhcLhcA1dXVuN1uhg4dSlxcHDU1Nb7yycnJREZG\nEhUVRU5OTufUXEREfNp1P/PBgweze/du+vfv75u3aNEiBgwYwKJFi1ixYgWHDx8mJSWF/Px85syZ\nw65du/B6vUydOpWCggL8/P73vqH7mXetPn36c+zY4a6uxkULCgrm6NHqrq7GRdH9zOVSdMj9zL+7\ngezsbDweDwAej4fMzEwAsrKySEpKIjAwEIfDQUREBHl5eRdbd+kEp4Pc7LZTd34jEuks7XpsnGEY\nTJ06FX9/fx599FHmz59PZWUlNpsNAJvNRmVlJQBlZWVMnDjRt25YWBher/ecbS5dutT3OjY2ltjY\n2EtohoiI9eTm5pKbm9uusu0K808++YRBgwbxzTff4Ha7iYqKarbcMIxWn/XX0rKmYS4iIuf67oXu\nsmXLzlu2Xd0sgwYNAiAkJIQZM2aQl5eHzWajoqICgPLyckJDQwGw2+2UlJT41i0tLcVut19wI0RE\npP3aDPMTJ05w7NgxAI4fP05OTg4jR44kPj6e9PR0ANLT00lISAAgPj6ejIwM6urqKCoqorCw0DcC\npqmzV/Pdcerfp09H/g5ERC5Zm90slZWVzJgxA4D6+nrmzp1LXFwc48aNIzExkbS0NBwOBxs2bADA\n6XSSmJiI0+kkICCA1NTUFrtZuvPn4caZNzcRkStFu4YmdvhONTyrS51+c+2+9YfuO7RV575cig4Z\nmigiIlcuhbmIiAUozEVELEBhLiJiAQpzERELUJiLiFiAwlxExAIU5iIiFqAwFxGxAIW5iIgFKMxF\nRCxAYS4iYgEKcxERC1CYi4hYgMJcRMQCFOYiIhagMBcRsQCFuYiIBSjMRUQsQGEuImIBCnMREQtQ\nmIuIWIDCXETEAtoV5g0NDcTExDB9+nQAqqurcbvdDB06lLi4OGpqanxlk5OTiYyMJCoqipycnM6p\ntYiINNOuMF+1ahVOpxPDMABISUnB7XZTUFDAlClTSElJASA/P5/169eTn5/Pli1bWLhwIY2NjZ1X\nexERAdoR5qWlpWzevJlHHnkE0zQByM7OxuPxAODxeMjMzAQgKyuLpKQkAgMDcTgcREREkJeX14nV\nFxERgIC2Cvz0pz/lpZde4ujRo755lZWV2Gw2AGw2G5WVlQCUlZUxceJEX7mwsDC8Xm+L213a5HXs\nmUlERP4nNzeX3NzcdpVtNcw3bdpEaGgoMTEx592gYRi+7pfzLW/J0nZVT0Sk4/Tp059jxw53dTU6\nRath/re//Y3s7Gw2b97MqVOnOHr0KPfffz82m42KigoGDhxIeXk5oaGhANjtdkpKSnzrl5aWYrfb\nO7cFIiLtdDrIza6uxiU4/4Vzq33my5cvp6SkhKKiIjIyMrjjjjtYt24d8fHxpKenA5Cenk5CQgIA\n8fHxZGRkUFdXR1FREYWFhbhcrg5siIiItKTNPvOmznaZLF68mMTERNLS0nA4HGzYsAEAp9NJYmIi\nTqeTgIAAUlNTW+2CERGRjmGYZ4eoXM6dGka3/0enCw5bhzn9Btt96w9Gtz3+Ove7lpXPfX0DVETE\nAhTmIiIWoDAXEbEAhbmIiAUozEVELEBhLiJiAQpzERELUJiLiFiAwlxExAIU5iIiFqAwFxGxAIW5\niIgFKMxFRCxAYS4iYgEKcxERC1CYi4hYgMJcRMQCFOYiIhagMBcRsQCFuYiIBSjMRUQsQGEuImIB\nCnMREQtoNcxPnTrFhAkTiI6Oxul08swzzwBQXV2N2+1m6NChxMXFUVNT41snOTmZyMhIoqKiyMnJ\n6dzai4gIAIZpmmZrBU6cOEHPnj2pr6/n1ltvZeXKlWRnZzNgwAAWLVrEihUrOHz4MCkpKeTn5zNn\nzhx27dqF1+tl6tSpFBQU4OfX/D3DMAxa3ekVzgDaOGxXNMMwoJv/Brrr8de537WsfO632c3Ss2dP\nAOrq6mhoaCA4OJjs7Gw8Hg8AHo+HzMxMALKyskhKSiIwMBCHw0FERAR5eXkd1QoRETmPgLYKNDY2\nMmbMGA4cOMBjjz3G8OHDqaysxGazAWCz2aisrASgrKyMiRMn+tYNCwvD6/W2uN2lTV7HnplERKSp\n3DNT29oMcz8/P/bt28eRI0eYNm0aO3bsaLbcMIwz/7q07HzLlrareiIiV7NYml/qLjtvyXaPZunb\nty933XUXu3fvxmazUVFRAUB5eTmhoaEA2O12SkpKfOuUlpZit9svpOYiInIRWg3zqqoq30iVkydP\nsnXrVmJiYoiPjyc9PR2A9PR0EhISAIiPjycjI4O6ujqKioooLCzE5XJ1chNERKTVbpby8nI8Hg+N\njY00NjZy//33M2XKFGJiYkhMTCQtLQ2Hw8GGDRsAcDqdJCYm4nQ6CQgIIDU1tdUuGBER6RhtDk3s\nlJ1qeFaXsvLwrCudzv2uZeVzX98AFRGxAIW5iIgFKMxFRCxAYS4iYgEKcxERC1CYi4hYgMJcRMQC\nFOYiIhagMBcRsQCFuYiIBSjMRUQsQGEuImIBCnMREQtQmIuIWIDCXETEAhTmIiIWoDAXEbEAhbmI\niAUozEVELEBhLiJiAQpzERELUJiLiFiAwlxExALaDPOSkhJuv/12hg8fzogRI1i9ejUA1dXVuN1u\nhg4dSlxcHDU1Nb51kpOTiYyMJCoqipycnM6rvYiIAGCYpmm2VqCiooKKigqio6Opra1l7NixZGZm\nsnbtWgYMGMCiRYtYsWIFhw8fJiUlhfz8fObMmcOuXbvwer1MnTqVgoIC/Pz+975hGAat7vQKZwBt\nHLYrmmEY0M1/A931+Ovc71pWPvfbvDIfOHAg0dHRAPTu3Zthw4bh9XrJzs7G4/EA4PF4yMzMBCAr\nK4ukpCQCAwNxOBxERESQl5fXUS0REZEWBFxI4eLiYvbu3cuECROorKzEZrMBYLPZqKysBKCsrIyJ\nEyf61gkLC8Pr9Z6zraVNXseemUREpKncM1Pb2h3mtbW1zJo1i1WrVhEUFNRsmWEYZ/59aVlLy5a2\nd8ciIletWJpf6i47b8l2jWb59ttvmTVrFvfffz8JCQnA6avxiooKAMrLywkNDQXAbrdTUlLiW7e0\ntBS73X5B1RcRkQvTZpibpsnDDz+M0+nkqaee8s2Pj48nPT0dgPT0dF/Ix8fHk5GRQV1dHUVFRRQW\nFuJyuTqp+iIiAu0YzbJz504mT57MqFGjfN0lycnJuFwuEhMT+frrr3E4HGzYsIF+/foBsHz5cn7/\n+98TEBDAqlWrmDZtWvOd6hP9LmXlT/SvdDr3u5aVz/02w7wz6ITuWlY+oa90Ove7lpXPfX0DVETE\nAhTmIiIWoDAXEbEAhbmIiAUozEVELEBhLiJiAQpzERELUJiLiFiAwlxExAIU5iIiFqAwFxGxAIW5\niIgFKMxFRCxAYS4iYgEKcxERC1CYi4hYgMJcRMQCFOYiIhagMBcRsQCFuYiIBSjMRUQsQGEuImIB\nCnMREQtoM8wfeughbDYbI0eO9M2rrq7G7XYzdOhQ4uLiqKmp8S1LTk4mMjKSqKgocnJyOqfWIiLS\nTJth/uCDD7Jly5Zm81JSUnC73RQUFDBlyhRSUlIAyM/PZ/369eTn57NlyxYWLlxIY2Nj59RcRER8\n2gzzSZMmERwc3GxednY2Ho8HAI/HQ2ZmJgBZWVkkJSURGBiIw+EgIiKCvLy8Tqi2iIg0FXAxK1VW\nVmKz2QCw2WxUVlYCUFZWxsSJE33lwsLC8Hq9LW5jaZPXsWcmERFpKvfM1LaLCvOmDMPAMIxWl7dk\n6aXuWETE8mJpfqm77LwlL2o0i81mo6KiAoDy8nJCQ0MBsNvtlJSU+MqVlpZit9svZhciInIBLirM\n4+PjSU9PByA9PZ2EhATf/IyMDOrq6igqKqKwsBCXy9VxtRURkRa12c2SlJTERx99RFVVFeHh4fz6\n179m8eLFJCYmkpaWhsPhYMOGDQA4nU4SExNxOp0EBASQmpraaheMiIh0DMM0TfOy79QwuOw77UAG\n0AWHrcOcfoPtvvUHo9sef537XcvK576+ASoiYgEKcxERC1CYi4hYgMJcRMQCFOYiIhagMBcRsQCF\nuYiIBSjMRUQsQGEuImIBCnMREQtQmIuIWIDCXETEAhTmIiIWoDAXEbEAhbmIiAUozEVELEBhLiJi\nAQpzERELUJiLiFiAwlxExAIU5iIiFqAwFxGxgE4J8y1bthAVFUVkZCQrVqzojF2IiEgThmmaZkdu\nsKGhgRtvvJFt27Zht9sZP3487733HsOGDfvfTg2DDt3pZWYAHXzYLivDMKCb/wa66/HXud+1rHzu\nd/iVeV5eHhERETgcDgIDA5k9ezZZWVkdvRsREWmiw8Pc6/USHh7u+zksLAyv19vRuxERkSYCOnqD\np/+NaUe5jt7xZdbedl65unf9u/Px7741P607H/vTunv9W9bhYW632ykpKfH9XFJSQlhYWLMy3bnP\nTUTkStTh3Szjxo2jsLCQ4uJi6urqWL9+PfHx8R29GxERaaLDr8wDAgJ47bXXmDZtGg0NDTz88MPN\nRrKIiEjH6/ChiSIicvnpG6AiIhagMBcRsQCFuYiIBSjMRUQsQGEul8zf35+YmBjf9PXXX5+37C23\n3NKubdbW1vLoo48SERHBuHHjuP3228nLy2t1nfZuOyMjg+XLl/POO+8QEhJCTEwMI0aM4L777uPk\nyZMAvPHGG6xbt65d27sY8+bNY+PGjS0uW7lyJcOGDSMmJgaXy9VmPebPn88XX3zRGdWUbqTDhybK\n1adnz57s3bu3XWU/+eSTc+bV19cTEND8VHzkkUcYMmQI+/fvB6C4uJj8/PwL3nZLtmzZwpNPPsnn\nn39OUlISq1evBmDu3LmsX7+eefPm8eijj7ZrWxfLMIwWv0n5u9/9ju3bt7Nr1y569+7NsWPH+OCD\nD1rd1ltvvdVZ1ZRuRFfm0uGOHz/O1KlTGTt2LKNGjSI7O9u3rHfv3gDk5uYyadIk7rnnHoYPH95s\n/QMHDpCXl8cLL7zgm+dwOLjzzjsBeOWVVxg5ciQjR45k1apVLW47NjaW++67j2HDhvHDH/7QV8Y0\nTfbt20dMTAymafq+jVxfX8/x48fp378/AEuXLuXll18GToely+UiOjqae++913f1/v777zNy5Eii\no6O57bbbgNN3DX366adxuVyMHj2aN99807ffH//4x0RFReF2uzl48GCL34ROTk7m9ddf97UlKCiI\nBx54AIDt27czZswYRo0axcMPP0xdXR0AsbGx7Nmzx3cMnn32WaKjo7nppps4ePBgu35nYgGmyCXy\n9/c3o6OjzejoaHPmzJlmfX29efToUdM0TfObb74xIyIifGV79+5tmqZp7tixw+zVq5dZXFx8zvay\nsrLMGTNmtLivTz/91Bw5cqR54sQJs7a21hw+fLi5b9++c7bdt29f0+v1mo2NjeZNN91k7ty50zRN\n09y9e7f5wAMPmKZpmmvXrjVDQkLM6Oho02azmZMnTzYbGhpM0zTNpUuXmitXrjRN0zQPHTrk2/+z\nzz5rrlmzxjRN0xw5cqRZVlZmmqZpHjlyxDRN03zjjTfMF154wTRN0zx16pQ5btw4s6ioyNy4caPp\ndrvNxsZGs6yszOzXr5+5cePGZm07cuSIGRwc3GK7T548aYaHh5uFhYWmaZrmAw88YL766qumaZpm\nbGysuXv3btM0TdMwDHPTpk2maZrmokWLfHUR69OVuVyya6+9lr1797J37142btxIY2MjzzzzDKNH\nj8btdlNWVtbiFaLL5eL6668/Z35rN3LauXMnM2fO5Nprr6VXr17MnDmTv/71ry1u+7rrrsMwDKKj\noykuLgZOd7GcvcIHmD17Nnv37qWiooIRI0bw0ksvnbOtf/7zn0yaNIlRo0bx7rvv+rp7brnlFjwe\nD2+//Tb19fUA5OTk8Ic//IGYmBgmTpxIdXU1hYWFfPzxx8yZMwfDMBg0aBB33HFH6wf1O7788ksG\nDx5MREQEAB6Pp8V29+jRg7vuuguAsWPH+tot1qcwlw737rvvUlVVxZ49e9i7dy+hoaGcOnXqnHK9\nevVqcX2n08lnn31GY2PjOcsMo/nN+U3TbDH8r7nmGt9rf39/GhoaANi6dStxcXEtbuvuu+9uFpBn\ntztv3jxSU1P5/PPPef75533dLK+//jovvPACJSUljB07lurqagBee+0135vbgQMHcLvdvrq2pk+f\nPvTu3ZuioqIW293U+bYVGBjoe+3n5+d7kxHrU5hLhzt69CihoaH4+/uzY8cOvvrqqwtaf8iQIYwb\nN47nn3/eN6+4uJjNmzczadIkMjMzOXnyJMePHyczM5NJkya1uU3TNDly5Aj19fUEBwf75jW1c+dO\n35Wv2aQ/vba2loEDB/Ltt9/yxz/+0Vf+wIEDuFwuli1bRkhICCUlJUybNo3U1FRfiBYUFHDixAkm\nT57M+vXraWxspLy8nB07drRYz2eeeYbHH3+cY8eO+fa9bt06brzxRoqLizlw4AAA69atIzY2tj2H\nU64SGs0il+y7V41z585l+vTpjBo1inHjxp3zyMDzrdfU22+/zc9//nMiIiK49tprGTBgACtXriQm\nJoZ58+bhcrmA08PyRo8e3ea2DcNg69atvqvks/PWr1/Pzp07aWxsJDw8nHfeece37Ow2fvOb3zBh\nwgRCQkKYMGECtbW1ACxatIjCwkJM02Tq1KmMHj2aUaNGUVxczJgxYzBNk9DQUDIzM5kxYwYffvgh\nTqeT//u//+Pmm29usd2PPfYYtbW1jB8/nsDAQAIDA/nFL37BNddcw9q1a7nvvvuor6/H5XLxox/9\nqNXfxflGzIg16UZbctWYP38+8+fP970RiFiJwlxExALUZy4iYgEKcxERC1CYi4hYgMJcRMQCFOYi\nIhagMBcRsYD/B4s27jhN1bS0AAAAAElFTkSuQmCC\n"
      }
     ],
     "prompt_number": 63
    },
    {
     "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",
     "collapsed": false,
     "input": [
      "# Show first 50\n",
      "print biasedcoin.pmf(btosses)[:50]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "[ 0.4  0.6  0.4  0.6  0.6  0.4  0.4  0.4  0.4  0.6  0.6  0.6  0.6  0.6  0.6\n",
        "  0.6  0.6  0.6  0.4  0.4  0.4  0.4  0.4  0.6  0.4  0.4  0.4  0.6  0.6  0.4\n",
        "  0.6  0.6  0.4  0.4  0.6  0.4  0.4  0.6  0.6  0.6  0.4  0.4  0.6  0.6  0.6\n",
        "  0.6  0.6  0.4  0.4  0.4]\n"
       ]
      }
     ],
     "prompt_number": 64
    },
    {
     "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",
     "collapsed": false,
     "input": [
      "# 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)))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 65
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print biasedcoin.likelihood(btosses)\n",
      "print biasedcoin.loglikelihood(btosses)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "1.16661962924e-292\n",
        "-672.200736793\n"
       ]
      }
     ],
     "prompt_number": 66
    },
    {
     "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",
     "collapsed": false,
     "input": [
      "mus = arange(0,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()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "pyout",
       "prompt_number": 67,
       "text": [
        "<matplotlib.legend.Legend at 0x4cf5a90>"
       ]
      },
      {
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAY8AAAEUCAYAAADTO7pnAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3Xl8TPf6wPHPSGJPIgtBgpBFZBGxKypKLGmlllarC1q0\nF1Wl997uv3Lbor1VrbZ0o3XVXrW0diUoGsRWgoQmZJOQBIkgy5zfH6cZQSKT5MxMluf9ep0XOXPm\nnGcOOc98d52iKApCCCFEKdSwdABCCCEqH0keQgghSk2ShxBCiFKT5CGEEKLUJHkIIYQoNUkeQggh\nSs1kyePmzZt06dKFdu3a4evryxtvvAFAeno6ISEheHt7069fP65cuWJ4z8yZM/Hy8sLHx4etW7ca\n9kdGRhIQEICXlxeTJ082VchCCCGMZLLkUbt2bXbu3MnRo0c5fvw4O3fu5Pfff2fWrFmEhIQQHR1N\nnz59mDVrFgBRUVGsWLGCqKgoNm/ezIQJEygYgjJ+/HgWLFhATEwMMTExbN682VRhCyGEMIJJq63q\n1q0LQE5ODvn5+Tg4OLB+/XpGjRoFwKhRo1i7di0A69atY8SIEdjY2ODu7o6npycREREkJyeTmZlJ\n586dARg5cqThPUIIISzDpMlDr9fTrl07XFxc6N27N35+fqSkpODi4gKAi4sLKSkpACQlJeHm5mZ4\nr5ubG4mJiffsd3V1JTEx0ZRhCyGEKIG1KU9eo0YNjh49ytWrV+nfvz87d+6843WdTodOp9PkWlqd\nRwghqpuyzFJllt5W9vb2PPzww0RGRuLi4sLFixcBSE5OplGjRoBaooiPjze8JyEhATc3N1xdXUlI\nSLhjv6ura5HXURRFNkXh3XfftXgMFWWTeyH3Qu7F/beyMlnyuHz5sqEn1Y0bN9i2bRtBQUGEhYWx\naNEiABYtWsTgwYMBCAsLY/ny5eTk5BAbG0tMTAydO3emcePG2NnZERERgaIoLF682PAeIYQQlmGy\naqvk5GRGjRqFXq9Hr9fz7LPP0qdPH4KCghg+fDgLFizA3d2dlStXAuDr68vw4cPx9fXF2tqaefPm\nGaqi5s2bx+jRo7lx4wahoaEMGDDAVGELIYQwgk4pT7mlAtHpdOUqglUl4eHhBAcHWzqMCkHuxW1y\nL26Te3FbWZ+dkjyEEKIaK+uz06S9rYQQVYOjoyMZGRmWDkOUg4ODA+np6ZqdT0oeQogSye9X5Vfc\nv2FZ/21lYkQhhBClJslDCCFEqUnyEEIIUWqSPIQQ1U5oaCiLFy822/X8/f3ZvXu32a5nDtJgLoQo\nUUX+/XJ3dyc1NRUrKytAjTU6OprGjRtrep3o6GjeeustwsPDyc3NpUWLFowePZrJkydTo0bF/x4u\nDeZCCFGITqfj119/JTMzk8zMTK5du1auxJGfn3/PvnPnztGlSxdatGjBiRMnuHLlCqtWrSIyMpLM\nzMzyhF9pSfIQQlQ5V65c4ZFHHqFRo0Y4OjoyaNCgO5ZyCA4OZsGCBQD88MMPdO/enalTp+Ls7Mz0\n6dPvOd+7775Ljx49+Pjjjw1LSnh7e/Pjjz9ib28PwPr16/Hz88PBwYHevXtz+vRpw/vd3d3ZsWMH\nANOmTWP48OGMGjUKOzs7/P39iYyMNNm9MBVJHkKISu/uahe9Xs+YMWO4cOECFy5coE6dOrz00kuG\n1+9eDuLAgQN4eHiQmprKm2++ec/5f/vtNx577LFirx8dHc1TTz3F3LlzuXz5MqGhoQwaNIi8vDzD\n9Qr75ZdfGDFiBFevXiUsLOyO2CoLSR5CiHLT6bTZykJRFAYPHoyDgwMODg4MHToUR0dHhgwZQu3a\ntalfvz5vvvkmu3btKvYcTZs2ZeLEidSoUYPatWvf83paWhpNmjQp9v0rVqzgkUceoU+fPlhZWfHP\nf/6TGzdusG/fviKP79mzJwMGDECn0/HMM89w7Nix0n9wC5PpSYQQ5WbJtnSdTse6det46KGHDPuy\ns7OZMmUKW7ZsMUyrkpWVhaIoRS4c16xZs/tew8nJiaSkpGJfT05Opnnz5nfE1KxZs2JXPS2o+gJ1\nue6bN2+i1+srRcN7gcoTqRBCGGn27NlER0dz4MABrl69yq5du+67+FFJK5H27duX1atXF/t606ZN\nOX/+vOFnRVGIj48vduG6qqBKJo8Wn7bg6s2rlg5DCGEhWVlZ1KlTB3t7e9LT04tsBC+N6dOns2/f\nPv7973+TkpICwNmzZ3n22We5du0aw4cPZ8OGDezYsYPc3Fxmz55N7dq1eeCBB7T4OBVSlUweNjVs\nSLmeYukwhBAW8sorr3Djxg2cnZ154IEHGDhwYLGli7sbz4vSqlUr9u/fT1xcHH5+fjRo0IDHHnuM\nTp06Ub9+fUPPq0mTJtGwYUM2bNjAL7/8grX1vS0DRV2vpOtXRFVykGD3hd2Z1WcWPVv0tHBUQlQN\nFXmQoDCODBI0QuP6jbmYddHSYQghRJVVJZOHSz0XqbYSQggTkuQhhBCi1Kpm8qjvQkqWJA8hhDCV\nKpk8GtdvLCUPIYQwoSqZPFzquUiDuRBCmFDVTB5SbSWEECZVNZPH3w3m0i9dCCFMo0omj3o162Fd\nw5rMnOq5SIsQ4rbx48fz/vvvm/WahdcL0VpFWdK2SiYPkHYPIaoLd3d36tati62tLY6OjjzyyCMk\nJCQYXp8/fz5vv/22WWMqacqT6OhoHn/8cRo2bEiDBg0IDAxkzpw56PX6Es994sQJHnzwQS3DLZOq\nmzyk3UOIaqHwMrTJycm4uLgwadIkS4dVrKqypG3VTR4yUFCIaqdWrVoMGzaMqKgow77Ro0fzzjvv\nAJCRkXHf5Wl/+OEHPDw8sLOzo1WrVixdutTw2sKFC/H19cXR0ZEBAwZw4cIFw2vbtm3Dx8eHBg0a\nMGnSpPtO/15VlrStuslDSh5CVBsFD+rs7GxWrFhBt27dDK8VrkJSFKXY5WmvX7/O5MmT2bx5M9eu\nXWP//v20a9cOgHXr1jFz5kzWrFnD5cuX6dmzJyNGjADg8uXLDBs2jBkzZpCWloaHhwd79+4tttqq\nqixpW2VXEmxcTwYKCmEuuunaTCmuvFv6HpIFy9BaW1tz/fp1GjVqxObNm+85BjAsT1vgzTffvGMF\nwho1avDnn3/i5uaGi4uLoWTw1Vdf8cYbb9C6dWsA3njjDWbMmMGFCxcIDw/H39+foUOHAup08LNn\nzy423tIsaQvwz3/+k88++4x9+/YV2dZRsKQtwDPPPMOnn35a/M3SUJVNHi71XTicfNjSYQhRLZTl\noa+VwsvQKorC2rVr6dWrF6dOnaJRo0Z3HHu/5Wnr1avHihUr+PjjjxkzZgzdu3dn9uzZtG7dmvPn\nzzN58mReffXVO86XmJhIcnIybm5ud+y/37K2VWVJ26pbbSVtHkJUOzqdjiFDhmBlZcXvv/9+x34o\neXnafv36sXXrVi5evIiPjw/jxo0DoHnz5nzzzTdkZGQYtuvXr9OtWzeaNGlCfHy84VoFS9AWp6os\naVt1k4e0eQhRbRQ8/BVFYd26dWRkZNCmTRvDvoLX77c8bWpqKuvWreP69evY2NhQr149rKysAPjH\nP/7BjBkzDA3xV69eZdWqVQCEhoZy8uRJ1qxZQ15eHnPnzuXixeKHCVSVJW2rbPKQyRGFqD4GDRqE\nra0t9vb2vPPOO/zvf/8zJI/CDeb3W55Wr9czZ84cXF1dcXJyYs+ePcyfPx+AwYMH89prr/Hkk09i\nb29PQEAAW7ZsAcDZ2ZlVq1bx+uuv4+zszNmzZ+nRo0exsVaVJW2r5DK0ANdzruP8X2ey38yulOsD\nC1GRyDK0lV+lWYY2Pj6e3r174+fnh7+/P3PnzgXUfslubm4EBQURFBTEpk2bDO+ZOXMmXl5e+Pj4\nsHXrVsP+yMhIAgIC8PLyYvLkyUZdv17NeljprGSKEiGEMAGT9baysbFhzpw5tGvXjqysLDp06EBI\nSAg6nY6pU6cyderUO46PiopixYoVREVFkZiYSN++fYmJiUGn0zF+/HgWLFhA586dCQ0NZfPmzYau\nafdT0O5hV8vOVB9TCCGqJZOVPBo3bmwYYFO/fn3atGlj6GpWVBFp3bp1jBgxAhsbG9zd3fH09CQi\nIoLk5GQyMzPp3LkzACNHjmTt2rXGxSDtHkIIYRJmGecRFxfHkSNH6Nq1K3v37uXzzz/nf//7Hx07\ndmT27Nk0aNCApKQkunbtaniPm5sbiYmJ2NjY3NGH2tXVtdj+ztOmTTP8PTg4WO2uKz2uhCg3BwcH\naTus5BwcHAAIDw8nPDy83OczefLIysriscce47PPPqN+/fqMHz+e//u//wPgnXfe4dVXX9Vs6uLC\nyQNgxYYVMrOuEBpIT0+3dAhCI8HBwQQHBxt+LtxduTRM2lU3NzeXYcOG8cwzzzB48GAAGjVqZOhe\nNnbsWA4cOACoJYrCA2sSEhJwc3PD1dX1jumVExISjB4sIwMFhRDCNEyWPAomIPP19eWVV14x7E9O\nTjb8fc2aNQQEBAAQFhbG8uXLycnJITY2lpiYGDp37kzjxo2xs7MjIiICRVFYvHixIRGVRJKHEEKY\nhsmqrfbu3cuPP/5I27ZtCQoKAmDGjBksW7aMo0ePotPpaNmyJV9//TUAvr6+DB8+HF9fX6ytrZk3\nb56hjnXevHmMHj2aGzduEBoaalRPK/i7wfycJA8hhNBalR0kCLAvfh9Tt0zlj7F/WCgqIYSo2Crc\nIMGKQKqthBDCNKp28vh7kGAVKVwJIUSFUaWTR/2a9amhq0FWTpalQxFCiCqlSicP+Lv0IVVXQgih\nqaqfPOq5yEBBIYTQWNVPHrIolBBCaK7YcR4lTUfg6OioeTCm0KR+E5Iyi18vWAghROkVmzzat29v\n6P974cIFw6RaGRkZtGjRgtjYWLMFWR7eTt5Ep0dbOgwhhKgQFAW0mOOy2GqruLg4YmNjCQkJ4ddf\nfyUtLY20tDQ2bNhASEhI+a9sJj7OPpy+fNrSYQghRIXQvj2c1uCRWGKbx/79+wkNDTX8PHDgQPbt\n21f+K5uJJA8hhLjtr7/AxaX85ykxeTRt2pT333/fUBL54IMPjJ7VtiJobt+c9BvpZN6S5WiFENXb\ntWuQlwcNGpT/XCUmj2XLlpGamsqQIUMYOnQoqampLFu2rPxXNpMauhpqu0eatHsIIaq3xERwddWm\nzaPEWXWdnJyYO3cumZnqN3dbW9vyX9XMCqquOjTtYOlQhBDCYhISoNDCrOVSYsnjzz//JCgoCD8/\nP/z8/OjQoQMnTpzQ5uomcvAg6PW3f/Zx8uF0mrR7CCGqt8REMyaPF154gU8++YQLFy5w4cIFZs+e\nzQsvvKDN1U0kNBQuXbr9szSaCyGEmUse2dnZ9O7d2/BzcHAw169f1+bqJuLmpt6kApI8hBBCfS5q\n1d+pxOTRsmVL3nvvPUNvq/fff59WrVppc3UTadYMCi2HjpeTF2fTz5Kvz7dcUEIIYWFmrbZauHAh\nqampDB06lGHDhnHp0iUWLlyozdVN5O6SR12bujSu35i4K3EWi0kIISxNy2qrEntbOTo68vnnn1eq\n3lZ3lzzgdtWVh6OHZYISQggLM2u1VWXsbXV3yQOk3UMIUb3dvAlXr0KjRtqcr0r2tiqy5CHddYUQ\n1VhSEjRtCjU0WoijSva2atZMSh5CCFGYllVWYESbR0Fvq2effRZFUViyZEmF723l6qr2KtDrb2dZ\nSR5CiOpMy8ZyqKK9rWrXBnv7OwcKNqrXiDx9HpezL1suMCGEsBAtu+lCKXpbVTZubmq7R8HUwzqd\nDh9nH85cPoNzc2fLBieEEGaWkADu7tqdr8TkcebMGT7++GPi4uLIy8sD1Afxjh07tIvCBAraPTp2\nvL2voOqqe/PulgtMCCEsICEBumv46CsxeTz++OOMHz+esWPHYmVlBajJo6IrKHkUJj2uhBDVldmr\nrWxsbBg/frx2VzST4gYKfnv4W8sEJIQQFmS2BvP09HTS0tIYNGgQX375JcnJyaSnpxu2iq6ogYLt\nGrfjUNIhFEWxTFBCCGEBeXmQkgJNmmh3zmJLHu3bt7+jeurjjz++4/XY2FjtojCBokoeze2bU0NX\ng/NXz+PewN0icQkhhLmlpICzM9jYaHfOYpNHXFycdlexgKJKHjqdjm7NurE/fr8kDyFEtaF1lRXc\nJ3ns2LGDhx56iNWrVxfZQD506FBtI9GYq6s6HL/wQEGAbm7d2J+wnxEBIywXnBBCmJHWo8vhPslj\n165dPPTQQ/zyyy+VMnkUDBRMTYXGjW/v7+rWlRUnV1guMCGEMDOzljymT58OwA8//KDtFc2ooN2j\ncPLo0KQDUZeiuJF7gzo2dSwXnBBCmInW3XThPslj9uzZ9+zT6XQoioJOp2Pq1KnaRmICBe0enTrd\n3lfHpg5+Df04lHSIni16Wi44IYQwk4QEaNtW23MW21U3MzOTrKysO7aCfQULQ91PfHw8vXv3xs/P\nD39/f+bOnQuoXYBDQkLw9vamX79+XLlyxfCemTNn4uXlhY+PD1u3bjXsj4yMJCAgAC8vLyZPnmz0\nhyuqxxWoVVd/JPxh9HmEEKIyM0WbB4qJJCcnK0eOHFEURVEyMzMVb29vJSoqSvnXv/6lfPjhh4qi\nKMqsWbOU1157TVEURTl58qQSGBio5OTkKLGxsYqHh4ei1+sVRVGUTp06KREREYqiKMrAgQOVTZs2\n3XO9oj7KzJmK8q9/3Rvb0uNLlSHLh2jyOYUQoqJr1UpRoqOLfq2saaDEWXXPnDlDnz598PPzA+D4\n8eO8//77JSalxo0b065dOwDq169PmzZtSExMZP369YwaNQqAUaNGsXbtWgDWrVvHiBEjsLGxwd3d\nHU9PTyIiIkhOTiYzM5POnTsDMHLkSMN7SlJcyaNbM7XHlSKDBYUQVZyiqG0eWpc8Skwe48aNY8aM\nGdSsWROAgIAAli1bVqqLxMXFceTIEbp06UJKSgouf0916+LiQkpKCgBJSUm4FWrRcXNzIzEx8Z79\nrq6uJCYmGnXdohaFAmhh3wKA81fPl+pzCCFEZXP5MtStq25aKnFuq+zsbLp06WL4WafTYVOKYYpZ\nWVkMGzaMzz77DFtb2zte0+l0mk6yOG3aNMPfg4ODad48uMiSh06no6tbVxksKISo8s6fV79IFwgP\nDyc8PLzc5y0xeTRs2JCzZ88afv7pp59oYuQEKbm5uQwbNoxnn32WwYMHA2pp4+LFizRu3Jjk5GQa\n/b0au6urK/GFnvQJCQm4ubnh6upKQqHiQ0JCAq7FlL8KJw+AW7cgOfnegYIggwWFENXD8eN39rQK\nDg4mODjY8HPBsIzSKrHa6osvvuDFF1/kzJkzNG3alDlz5jB//vwST6woCmPGjMHX15dXXnnFsD8s\nLIxFixYBsGjRIkNSCQsLY/ny5eTk5BAbG0tMTAydO3emcePG2NnZERERgaIoLF682PCektSqBQ0a\nqPO63K2bWzfpcSWEqPKOHoW/m5+1VVKL+qFDhxRFUXtMXb16VVEURfnll19KbInfs2ePotPplMDA\nQKVdu3ZKu3btlE2bNilpaWlKnz59FC8vLyUkJETJyMgwvOeDDz5QPDw8lNatWyubN2++IwZ/f3/F\nw8NDmTRpUpHXK+6jdOigKH931LrD9ZzrSp336yjZOdklfhYhhKisevZUlO3bi3/diDRQJN3fby5W\n+/btWbRoEQEBAQAsW7aMOXPmcODAAROksrIrGMB4t8GDYeRIKGo2le4Lu/N/D/4f/T37myFCIYQw\nL70eHBzg3Dl1Vt2iFPfsLEmJ1VY//fQTo0aN4vTp03z77bfMmzePbdu2lfpCluLurt64ooR6hrLp\n7CazxiOEEOYSFwd2dsUnjvIoMXm0atWKZcuWMWTIEFavXs2WLVuwt7fXPhITaddOrfMrSqhXKBtj\nNpo3ICGEMBOTtXdwn95WBdVUBdLT09Hr9XTp0gWdTsfx48dNE5HGgoLgo4+Kfq1d43Zk5WQRkxaD\nl5OXeQMTQggTO3YMAgNNc+5ik8cvv/ximiuama+vWnS7fh3q1bvzNZ1OxwDPAWw6u0mShxCiyjl6\nFJ591jTnLrbaytHREXd3d+zs7IrcKgsbGzWBFFdQCvWSdg8hRNVkkWqrESNGsGHDhnvWMi9Q0dcw\nLywoCI4cgW7d7n0tpFUIz617juzcbOraaDx+XwghLCQ9HTIyoFUr05y/2OSxYcMGoPKvZQ63k0dR\n7Gvb06FJB3bG7uRh74fNG5gQQpjIsWPqyPK7Z9fQSrHJ4/Dhw/d9Y/v27TUPxlTat4cFC4p/PdQr\nlI1nN0ryEEJUGaassoL7JI+pU6fed9LCnTt3miQgU2jbFk6dgtxctQ3kbgM9BxK2PMywSqIQQlR2\nR49CzyIWS41MisTH2Yd6Nevd+2IpFJs8tJh1saKoW1cdLBgVVXS3Nf9G/uTp8ziTdgYfZx+zxyeE\nEFo7dgwmTbp3/2OrHmP7s9vxcPQo1/lNVBtW8bRvD8XVxOl0OkK9QvnlTNXoniyEqN5yciA6Gv5e\nw+/2/vwckjKTaG7fvNzXqDbJ436N5gDDfYez7ETpFrkSQoiKKCoKWraEOnXu3B93JY5mds2wsTJ+\nTabiSPL4W7B7MCnXUzh16ZT5ghJCCBMorrH8bPpZPB09NblGiYtBRUZG3tOIbG9vT4sWLbC2LvHt\nFUZQkFoHWNTCUABWNax40v9Jlvy5hPcfKnmNdiGEqKgOHFCr6u92Nv1suds6CpRY8pg4cSJdunRh\n3LhxjBs3jq5du/LYY4/h7e3Nli1bNAnCHBwcwMkJCi2KeI+nA55m6Z9LyzQ9sRBCVBS//QZ9+ty7\n/1zGOTwdtCl5lJg8mjZtytGjR4mMjCQyMpKjR4/SqlUrtm3bxr///W9NgjCXkqqughoHUcu6FvsT\n9psvKCGE0NCFC+ro8sJLzxYwa8njzJkz+BVqsvf19eX06dN4eHhUujERJSUPnU7H0wFPs+TPJeYL\nSgghNFRQ6iiqev5c+jnN2jxKTB5+fn6MHz+eXbt2ER4ezoQJE/D19eXWrVvYFDXirgJr3x4OHbr/\nMU8FPMWqk6vIzc81T1BCCKGhbdugb9979+fr84m7EkfLBi01uU6JyeOHH37Aw8ODTz/9lM8++4xW\nrVqxaNEibGxs2LFjhyZBmEv37mpD0s2bxR/TyqEVno6ebPur8qyWKIQQoHYI+u03CAm597X4a/E0\nrNeQOjZ17n2xDErsLlW3bl0mTZpE//7qOt+tW7emZs2aANja2moShLk0aAABAfD770Vn5gJPBTzF\nj8d/JNQr1HzBCSFEOZ04oS4726LFva+dSz+Hh4M27R1gRMkjPDwcb29vJk6cyMSJE/H29mbXrl2a\nBWBu/frB1q33P+ZJ/yfZdHYTl7MvmycoIYTQwPbtxX8x1nKMBxiRPKZOncrWrVvZvXs3u3fvZuvW\nrUyZMkWzAMytXz8oqYexc11nBvsM5rvD35knKCGE0EBx7R3wdzddcyaPvLw8WrdubfjZ29ubvLw8\nzQIwt06d1K5sycn3P+6lTi8x/9B88vSV97MKIaqPW7dg71546KGiXz+bfta81VYdOnRg7NixhIeH\ns3PnTsaOHUvHjh01C8DcrK3Vm7t9+/2P69C0A662rjJZohCiUvjjD/DxUQdEF8Xs1Vbz58+nTZs2\nzJ07l88//xw/Pz/mz5+vWQCW0L9/yVVXAC91fokvDn5h+oCEEKKc7tfeoSgK5zLOaTZAEECnVJG5\nOHQ6ndHTisTGQteuatXV/ZZozMnPocWnLfht5G/4NvTVKFIhhNBe164wY0bR1VbJmckEfhVI6r9S\n73mtNM/OwortqhsQEFDsm3Q6HcePHy/1xSqKli3B3h6OH7//Mo01rWryQocX+OLAF8x7eJ75AhRC\niFKIj4eYGHUsW1G0LnXAfZLHL79U7br+gl5XJa3x+2KHF/Gb58eMPjNoULuBeYITQohSWLEChgyB\nWrWKfl3r9g64T5uHu7v7fbfKrn//ksd7ADS1bUpY6zC+OCBtH0KIimnZMhgxovjXtZxNt0C1WQzq\nbsHBEBEBWVklH/tmjzf5LOIzMm9lmjwuIYQojehoSEpSn2nF0XI23QLVNnnY2kLPnmBM7Vxr59b0\n8+jHlwe/NH1gQghRCsuXw+OPg5VV8ceYtdqqsOzsbM6cOaPphSuCp5+GH3807ti3er7FnD/mkJVj\nRFFFCCHMQFFKrrIC7ee1AiOSx/r16wkKCjJMjHjkyBHCwsI0DcJSBg9WR2Sm3tt77R6+DX0Jdg9m\n/sHKPcZFCFF1HDumzhLetWvxx6TfSCdfyce5rrOm1y4xeUybNo2IiAgc/h62GBQUxF9//aVpEJZS\nvz488gisXGnc8e88+A6z988mOzfbtIEJIYQRli2DJ5+E+63LdzL1JK2dWmu+eF+JycPGxoYGDe7s\nolrjfiPrKplnnjG+6sq/kT89mvdg3kEZ8yGEsCy9Xm3vKKnKas+FPfRo3kPz6xu1kuCSJUvIy8sj\nJiaGSZMm8cADD2geiKX07auOOI+JMe7493q/x4d7PyQtO820gQkhxH3s2aPWntxnPLd63IU9PNji\nQc2vX2Ly+Pzzzzl58iS1atVixIgR2NnZ8emnnxp18ueffx4XF5c7RqtPmzYNNzc3goKCCAoKYtOm\nTYbXZs6ciZeXFz4+PmwtNAgjMjKSgIAAvLy8mDx5cmk+X4msrdVi3xIjly1v07ANT/g9wfRd0zWN\nQwghSuOLL2D8+PtXWeXp89gXv88kJQ+UEkRGRpZ0SLF2796tHD58WPH39zfsmzZtmjJ79ux7jj15\n8qQSGBio5OTkKLGxsYqHh4ei1+sVRVGUTp06KREREYqiKMrAgQOVTZs23fN+Iz5KsQ4cUBRPT0X5\n+3IlunT9kuL8kbNy6tKpMl9TCCHK6vx5RXFwUJRr1+5/3KHEQ4rfl373Paasz06jFoPy8fHhnXfe\n4cSJE6VKTD179jQ0tN+VsO7Zt27dOkaMGIGNjQ3u7u54enoSERFBcnIymZmZdO7cGYCRI0eydu3a\nUsVRko6J++oCAAAgAElEQVQd1QkSDxww7njnus683v11/r3t35rGIYQQxpg/H0aOVMer3c/u87tN\nUmUFRi5Du3PnTpydnXnxxRcJCAjgvffeK9dFP//8cwIDAxkzZgxXrlwBICkpCTc3N8Mxbm5uJCYm\n3rPf1dWVxMTEcl3/bjodjBoF33xj/Hte6vwSJy+d5Le/ftM0FiGEuJ8bN+C77+Cll0o+dvcFCyYP\ngCZNmjB58mS++uorAgMD+c9//lPmC44fP57Y2FiOHj1KkyZNePXVV8t8rrtNmzbNsIWHh5fqvePG\nwc8/GzfmA6CWdS0+6vsRU7ZMITc/t/TBCiFEGSxdCl26gGcJA8b1ip495/fQs3nPO/aHh4ff8aws\nq2Jn1S0QFRXFypUr+emnn3BycuKJJ57gk08+KfMFGzVqZPj72LFjGTRoEKCWKOLj4w2vJSQk4Obm\nhqurKwkJCXfsd3V1LfLc5bkRDRvC8OFqcfDdd417z9A2Q/n28Ld8sv8TXuvxWpmvLYQQxlAUmDsX\n/vvfko89dekU9rXtcbW783kZHBxMcKGJsKZPL1vnnxJLHs8//zwNGjRgy5Yt7Nq1iwkTJtyRAEor\nudDi4WvWrDH0xAoLC2P58uXk5OQQGxtLTEwMnTt3pnHjxtjZ2REREYGiKCxevJjBgweX+fr388or\navK4edO443U6HfMfns9/9/2Xs+lnTRKTEEIU2L0bcnIgJKTkY03VRbdAiSWPP/74o8wnHzFiBLt2\n7eLy5cs0a9aM6dOnEx4eztGjR9HpdLRs2ZKvv/4aAF9fX4YPH46vry/W1tbMmzfPMCJy3rx5jB49\nmhs3bhAaGsqAAQPKHNP9tGkD7durxcLnnzfuPS0dWvJGjzf4x6//YNuz2zQfxSmEEAU++ggmT75/\n99wCu8/vJqSVEVmmjIpdhvbxxx9n1apVRa4oWBFXEizrUop327YNpkyBP/807h8I1L7UXb7rwsud\nX2ZUu1HljkEIIe62f786Ji06uvhFnwooikKzOc3YNXpXiVOxa74M7WeffQbAr7/+es+Jq/K36759\n1aSxfbtxRUMA6xrWfDvoWwYuGcgAzwG41HcxbZBCiGrnnXfUraTEARB7JRYFhVYOrUwWT7FtHk2b\nNgXUKqO7VxGcN6/qzu2k06klj48/Lt372jdpz9j2YxmzfowmJSAhhCgQHg5xceqQAmMUjO8w5Rf9\nEhvMtxaxVuvGjRtNEkxF8fTTatFw9+7Sve/dXu+Scj2F+Ydk2nYhhDYURS1xTJsGNjbGvWdjzEYe\ncn/IpHEVmzzmz59PQEAAZ86cISAgwLC5u7vTtm1bkwZlabVqwX/+A2+8of7DGaumVU2WDF3Cu+Hv\ncurSKdMFKISoNrZuhbS0kmfPLZB5K5Mt57YwtM1Qk8ZVbIP51atXycjI4PXXX+fDDz80VMXY2tri\n5ORk0qDKQqsG8wL5+dCuHcyYAX8PRTHaN5HfMP/QfP4Y8we1rI2ooBRCiCLo9epCT//6l7rUrDGW\nHF/C0hNL2fDUBqOOL+uzs9jkcbfU1FRuFhoA0bx581JfzJS0Th6grm/+5ptw9Oj91we+m6IoDFkx\nhFYOrfikf9kHVAohqrfvv4evv4Z9+9T594zxyNJHeNL/SZ5p+4xRx5f12WnUMrReXl60bNmSXr16\n4e7uzsCBA0t9ocrokUfAzk4d91EaOp2OhY8uZO3ptaw8aeQyhUIIUUhGhlp1Pm+e8Ykj/UY6ey7s\n4dHWj5o2OIxIHm+//Tb79+/H29ub2NhYfvvtN7p06WLywCoCnQ5mzYL/+z+4dat073Ws48jq4auZ\nuHEiUZeiTBOgEKLKevttGDZMHbhsrNVRq+nn0Q/bWiVMt6sBo5ahdXZ2Rq/Xk5+fT+/evTl06JDJ\nA6soevaEwEDj5pK5W1CTIP4b8l+GrhjKtVvXtA9OCFElRUbC6tXw/vule9+yE8sY4W9ky3o5lZg8\nHBwcyMzMpGfPnjz99NO8/PLL1K9f3xyxVRiffQaffgpnyzB91eh2o+ndsjej145Gr+i1D04IUaXo\n9TBxIsycCUUsh1Ss5Mxkjlw8QqhXqOmCK6TEBvOsrCzq1KmDXq9nyZIlXLt2jaeffrrC9bgyRYN5\nYR9/rHaZ27LF+GlLCtzKu0XfxX3p0bwHM/vMNE2AQogq4csv1WWxf//d+LYOgM/++IzDFw+zaPCi\nUl3P5L2tKjpTJ4/cXHXFwddfN76/dWGXsy/TbUE3Xuv+GmPbj9U+QCFEpXfmDPToAXv3grd36d7b\n5bsuTA+ezgDP0k0ca7LeVra2tvdsbm5uDBkyhL/++qvUF6ysbGzULnOvvgp/L35YKs51ndn41Ebe\n3vE2285t0z5AIUSllpenLi07bVrpE0dEQgSp11Pp26qvSWIrSoklj7fffptmzZox4u+v28uXL+fc\nuXMEBQXx1VdflXrFPlMxdcmjwIQJ6nofCxeW7f17zu9h2MphbB+5nbYuVXukvhDCeO+9B3v2lK1q\n/PFVj9OzeU9e7vJyqa9rsmqrtm3b3jP9ert27Th69CiBgYEcO3as1Bc1BXMlj6wsCApSG7Mee6xs\n51h5ciVTtkwhfFQ4Xk5e2gYohKh0IiNh4EA4cgSKWSi1WOfSz9F1QVdiJ8dSv2bpOzOZrNqqbt26\nrFixAr1ej16vZ+XKldSuXdtw0eqmfn110ODEiVBo1dxSGe43nGm9phGyOISEawklv0EIUWVduaKu\n0/H556VPHABz/pjDCx1eKFPiKI8SSx7nzp1j8uTJhhUFu3btyqeffoqrqyuRkZH06NHDLIGWxFwl\njwIzZ6rFy99+K93UJYV9vO9jFhxZwO7Ru2lYr6G2AQohKjy9HgYPBnd3dW3y0rqcfRmvz72ImhBF\nE9smZYpBeluZOXnk50OfPtC/vzqFQFm9s/Md1p9Zz/Znt0sCEaKa+eAD2LgRdu6EmjVL//73dr1H\n3NU4FoQtKHMMJqu2OnPmDH369MHPzw+A48eP835phz1WQVZWsHixOoDwt9/Kfp7/BP+HQd6D6L2o\nNylZKdoFKISo0LZuVcd0rFpVtsRxI/cGXx78kle7vap9cEYoMXmMGzeOGTNmUPPvTxcQEMCyZctM\nHlhl0KwZLF+uLh5V1l7LOp2O93q/x2O+j9F7UW8uZl3UNkghRIUTHQ3PPgvLlsHfi7aW2pcHv6Sr\nW1d8G/pqG5yRSkwe2dnZd0yEqNPpsDF2OatqIDgY3npLrbfMyirbOXQ6HdOCpzHCfwS9fujFhasX\nNI1RCFFxpKZCaKg6b1WvXmU7x6Xrl/hw74d82PdDbYMrhRKTR8OGDTlbaFKnn376iSZNytYwU1W9\n9JI6+nz06NKtPHi3d3q9w/iO4+mxsAcnU09qFp8QomLIzoawMLV31bhxZT/Pu+Hv8lTAU7R2bq1d\ncKVkVG+rF154gX379uHg4EDLli1ZsmQJ7u7uZgrROOZuML/brVvQuzc8+KA6jXt5LDm+hKlbp7Lm\niTU80OwBbQIUQlhUfr46NszWFhYtKv1AwAInU0/Se1FvTr90Gsc6juWOy+S9ra5fv45er8fW1vTz\nxJeFpZMHqOsM9+gBL7wAU6aU71ybz27m2TXP8t2g73jUx/QLuwghTEevV58LsbGwaVPZGsgLDPhx\nAAM9BzK562RNYivrs9O6pANu3rzJ6tWriYuLIz8/H0VR0Ol0/N///V+ZAq3KnJzUsR89ekCjRmpD\nelkN8BzAxqc2MmTFEGLSY3i126vVclCmEJWdoqhV26dPw+bN5Uscm2I2EXsllgmdJmgXYBmVmDwe\nffRRGjRoQIcOHQwjy0XxmjdXv1k89BA4OqpTDpRVJ9dO7B+zn0HLBnH68mnmPTyPmlbl+J8nhDAr\nRVFrISIjYds2dYaKsrp26xoTNk5g/sPzsbGyfKelEqut/P39OXHihLniKbOKUG1V2P798Oij8L//\nwYDSzZB8j6ycLJ5a/RRXb11l5WMrcanvok2QQgiTURT4179gxw51LFhpFnYqyrhf1Bb2bwd9q0F0\nt5lskOADDzxwz8SIomTdusHateoUyxs3lu9c9WvWZ80Ta+jVohedvu1EREKENkEKIUwiPx/+8Q/Y\nvVstcZQ3cWyI3sD2v7bzSb9PtAlQAyWWPNq0acPZs2dp2bIltWrVUt+k01W4hFLRSh4F/vhD7Zq3\ncCE88kj5z7f+zHrGrh/L+w+9z7j246QdRIgKJidHHQB46RKsW6f2riqPtOw02n7VlqVDl9LLvYwD\nQ+7DZL2t4uLiitwvXXWNFxEBgwbBnDnla0QvEJ0WzbCVw2jr0pavHv4K21oVswecENVNVhY8/jjU\nqqXOPlHeZmJFUXhy9ZM0tW3KnP5ztAnyLjIxYgVOHgAnTqijSidPVlcjLK/s3GymbJnCjtgdLB+2\nnA5NO5T/pEKIMktIUL8kdugAX30F1iV2RyrZ3Ii5LDiygD/G/EEdmzrlP2ERJHlU8OQB6vofAwao\nM/F+/HHpFrcvzooTK3hp00u81v01pnSdglWNMs4PL4Qos8hItYPM5Mnwz3+WfQBgYTtjdzJi9Qj2\nj9lPS4eW5T9hMSR5VILkAZCers6D5eSk9sTSYsxlbEYsI9eOpIauBosGL8K9gXv5TyqEMMqqVery\n1N98A0OGaHPOuCtxdP2uK0uGLqFPqz7anLQYJuttJbTl6Ajbt4OzMzzwQNln4y2spUNLwkeF87DX\nw3T6thMLjyysFIlUiMosNxemToV//1sdHKxV4sjOzWbIiiG81v01kyeO8pCSh4UoijqX//vvw5Il\n6sJSWjiecpzn1j2HUx0nvhn0jZRChDCB5GR44gl10N+PP6pfCrWQm5/L4BWDaVi3Id8/+r1ZelNK\nyaOS0enUKQuWLlW79f3nP2rf8PJq69KWiLER9GnZh07fdmJuxFzy9RqcWAgBqFOMdOigfuH79Vft\nEke+Pp9Ra0dhpbPi20HfVvhu+CZNHs8//zwuLi4EBAQY9qWnpxMSEoK3tzf9+vXjypUrhtdmzpyJ\nl5cXPj4+bN261bA/MjKSgIAAvLy8mDxZm8nAKoqHHoJDh9RlKPv3hxQNFhO0rmHNaz1e4/fnfmf1\nqdV0+a4LBxMPlv/EQlRjN2+qU4288IL6pe/dd7Xp9AJql9xJmyaRlJnEisdWVIjpR0pi0uTx3HPP\nsXnz5jv2zZo1i5CQEKKjo+nTpw+z/p6/PCoqihUrVhAVFcXmzZuZMGGCoSg1fvx4FixYQExMDDEx\nMfecs7Jr2lRtB+neHYKC1G8zWmjt3JrwUeG83OVlwpaHMWHDBDJuZGhzciGqkWPHoEsXtcfk0aPq\nInBaURSF1397nYjECNaPWG+yLrmaU0wsNjZW8ff3N/zcunVr5eLFi4qiKEpycrLSunVrRVEUZcaM\nGcqsWbMMx/Xv31/Zv3+/kpSUpPj4+Bj2L1u2THnxxRfvuY4ZPopZ7NqlKO7uijJ2rKJcu6bdedOz\n05Xxv45XGv23kTLvwDwlNz9Xu5MLUUXduqUo776rKA0bKsrChYqi12t7/nx9vjJhwwSl/dftlUvX\nL2l7ciOV9dmpwTCW0klJScHFRZ3Yz8XFhZS/62mSkpLo2rWr4Tg3NzcSExOxsbHBzc3NsN/V1ZXE\nxMQizz1t2jTD34ODgwnW8uuBmTz4oPotZ+pUCAyE775Tq7bKy6GOA/MenseLHV7klS2vMP/QfOb0\nn1Ohe3MIYUmHDsGYMdCsGRw5Aq6u2p4/T5/Hc+ue4/yV8+wYuQP72vbaXqAY4eHhhIeHl/s8Zk8e\nhel0Ok0bhQonj8rMzk5NGr/+qi5t27evOqhQi4a5wMaB7Bi5g59P/cwLv75Aa6fWfNj3QwJcAkp+\nsxDVwJUr8Pbb8NNP8N//wjPPaDPor7Ds3GyeWv0Ut/JvsfmZzdS1qavtBe7j7i/W06dPL9N5zN7b\nysXFhYsXLwKQnJxMo0aNALVEER8fbzguISEBNzc3XF1dSUhIuGO/q9ZfASqoRx6BkyfV7oB+fmqX\nXi16I+t0Oob5DuPUxFMM9BxI38V9Gb12NHFX4sp/ciEqKUVRf8d8fdUxHFFRak9IrRNH4rVEHvz+\nQWxr2bLuyXVmTRxaMnvyCAsLY9GiRQAsWrSIwYMHG/YvX76cnJwcYmNjiYmJoXPnzjRu3Bg7Ozsi\nIiJQFIXFixcb3lMd2NrC3LmwZg3Mnq1Wax09qs25a1rVZFKXSUS/FE0z+2Z0+KYDEzdOJPFa0dWC\nQlRVERHqoN1PPoHVq+Hrr7XrglvYwcSDdPmuC4/5Psb/Bv+vci/upm3Ty52efPJJpUmTJoqNjY3i\n5uamLFy4UElLS1P69OmjeHl5KSEhIUpGRobh+A8++EDx8PBQWrdurWzevNmw/9ChQ4q/v7/i4eGh\nTJo0qchrmfijVAh5eYryzTeK4uKiKP/4h6KkpGh7/tSsVOWfW/+pOMxyUF7e9LKScDVB2wsIUcHE\nxirK008rStOmivLDD4qSn2+6a/1w5Ael4UcNlbWn1pruImVQ1menjDCvhDIy1EGFixerE7FNnQr1\n6ml3/uTMZD7e/zHfH/meJ/yf4LXur8lIdVGlpKTABx+o1VQvvaSu+FeeJWLvJysni4kbJ3Iw8SAr\nH1+JfyN/01yojGSEeTXi4KCuDXLgAJw6Bd7e6lQnt25pc/4mtk2Y3W82Z146g0NtBzp804Fnfn6G\nYxePaXMBISzk8mV46y21XaNGDfX3Z/p00yWOYxeP0fGbjljprDg47mCFSxzlIcmjEmvVSh3pun49\nbNoEXl5qXW1Ojjbnb1ivITP6zOCvl/+irUtbQpeGMuDHAWw9t7XalPJE1XDpErz+OrRurSaQyEj4\n9FP4u7+O5nLzc3l/9/v0XdyXt3q+xcJHF1KvpobVAxWAVFtVIRER6pQJUVHqglNjx2pbnXUr7xZL\n/1zKpxGfkqfPY3KXyTzT9plK21tEVH2xsWoj+JIl6kSGb7wBzZub9poFk5M2rNuQbwd9SzP7Zqa9\nYDnJeh6SPAwOHoRZs+D332HiRHWtAWdn7c6vKArhceHM+WMO++L3MTJwJP/o+A+8nby1u4gQ5XDw\noFq1u3UrjBsHL78MTZqY9ppZOVm8t/s9Fh5ZyKw+s3g+6PkKP7khSJuHKKRTJ7W74a5dcOGCWp31\nwgtqiUQLOp2O3i17s37Eeg69cIiaVjXpsbAHIYtDWHVyFTn5GtWbCVEKubmwbBl066aWMjp2VNfL\nmTnTtIlDURRWnVxFmy/bkJSZxPF/HGdM+zGVInGUh5Q8qoHUVHVN5fnz1cGG48dDWBjYaDhx5828\nm/x86me+PfwtJ1NPMjJwJM8HPY9vQ1/tLiJEEWJj4dtv4fvvoU0btQfiI4+AlRlWZD6QeIB/bfsX\n6TfS+TL0Sx5s8aDpL6oxqbaS5FGiW7fg55/VJHL2rDpvz3PPqQ3vWopJi2HBkQUsPr4YV1tXRrcb\nzZP+T+JYxwSjrkS1dPOm2lHk++/VKqpnn1VL123amOf6Z9PP8taOt9h7YS/Tg6czqt0orGtYdLan\nMpPkIcmjVE6cUL+tLV0KAQFqEhk6VNsG9nx9Ptv/2s4Px35gY8xGerXoxdMBTzOo9SBpZBelpihq\n9/TFi2H5cmjX7vb/2zpmmsU8Ji2GD/Z8wK/RvzKl6xSmdJtS6f8vS/KQ5FEmt27BL7+o3+D27lWL\n+08/DSEhYK3hF6lrt66x5tQalp5YSkRCBKFeoTzh9wT9PftT27q2dhcSVU5UlNqWsXSp+n/ymWdg\n5Eho0cJ8MfyZ8icf7fuIzWc3M6nzJF7u8jINajcwXwAmJMlDkke5pabCypXqmsznzsGQITB8uLrw\njZaJJCUrhZ9P/czKqJUcvXiUgZ4DGdpmKAM8B1C/polGa4lKQ1HUCUFXrVK3zEx4/HH1S0379tpP\nVFh8HArb/9rO7P2zOZ5ynEmdJzGx80TsatmZJwAzkeQhyUNTcXG3f3ljY9USyeDBaomkroal9ItZ\nF1l3eh0/n/6Z/fH7CXYPJqx1GI94P0Lj+o21u5Co0PLzYf9+WLtW3XJzYdgwNWl06aLdcq/GuHrz\nKv879j++ivwKHTr++cA/GeE/glrWtcwXhBlJ8pDkYTLnz8O6dep28CD06gWhofDww9oOuMq4kcGm\ns5tYf2Y9W85twdvJm1DPUEK9QunQtAM1dNKzvCq5dEkdh7FhA2zZoi66NHgwPPqo2p5hzp6uiqKw\nP2E/3x/9np+ifqKfRz/GdxxPrxa9qn6XW0kekjzMISND/UXfsAE2b1andwgJgX791KSiVYN7Tn4O\nv1/4nU1nN7ExZiOXsy8T0iqEfh79CGkVQhNbE4/4Epq7cQP27YNt29Tt7Fno3Vv9EjJwIBRaMNRs\nzl85z9I/l/LDsR/QoWN0u9GMbje6WpV6JXlI8jC7/Hw4fFh9EGzdqi7b2a6d+kDo3Ru6dtWuiuv8\nlfNsPbeVrX9tZUfsDprUb8JDLR/ioZYP0atFLxzqOGhzIaGZGzfUKXN274YdO9T/HwEBt79sdOmi\n7VgjYyVlJvFT1E8sP7Gc6LRohvkOY3TgaLq6da3ypYyiSPKQ5GFx16+rPbZ27lS3P/+Etm2hRw/o\n3l1NJo01+EKXr8/nyMUj7IjdwW+xv7Evfh+ejp70atGLB1s8SI/mPWhUz0Qz3oliJSer7Rb796v/\nD44fV5NFz57ql4kePdTFzcxNURSi06JZe3ota06v4UzaGR7xfoQR/iPo26pv5V6QSQOSPCR5VDjZ\n2eo3zz171AfKH39AgwbqN85OndTpI9q3L/8DJSc/h8ikSHaf382u87vYF7+PRvUa0b15dx5we4Cu\nbl3xbeiLVQ0zDDmuJjIy1FLnoUNqO9jBg+qXh65d1a17d/XfWcvOFaWRnZvN7vO72RizkY0xG7mR\nd4Ow1mEM8RlCsHtwtU8YhUnykORR4en1EB2tJpRDh9Tt+HG1rjsoSK3yCgxUv626upa9wVSv6DmZ\nepK98XvZF7+PiMQIkjOT6di0I51cO9Gpqbo1t29eLaspSkOvV3vbHT+uliSPHoUjR9RpzQMD1S8B\nBZunp3kbuQsr+AKxI3YH22O3czDxIEFNggwdLtq6tJV/62JI8pDkUSnl5sKZM+oD6ciR2w+pnBzw\n91enm/D1Vf9s3VrtkVOWOYvSstM4kHiAg0kHOZR0iINJB8nNzyWoSRBBjdWtrUtbWju3rrTTTJTH\nrVvq2J7oaDh9Wh2YFxWl/t3JSa1+DAhQE0ZQkJoozNl99m4ZNzI4kHiA/Qn72XNhDwcSD+Dp6Emw\nezB9W/blwRYPYlvLAnVklZAkD0keVUpqqjqFyqlTt7foaPUbb6tW6kzBrVqBh4f6Z4sW6laa3l7J\nmckcuXiEw8mHOXrxKMdTjpNwLQEfZx/8G/nj19APv0Z++Db0pYV9i0pd7aUokJ6ujt+Ji1NLE+fO\nqT2ezp2DpCS123Xr1urm63s7advbWzb26znXOXrxKJHJkUQmR3Iw8SDx1+Lp2LQjXd260rN5Tx5o\n9kCVGfFtbpI8JHlUC9nZ6gMvJkadbvvcOfXP8+fV6efr1VMfgq6uanWYq6s6HXfB5uKirm1Ss5gq\n7+s51zmReoKTl05y8tJJTqSe4NSlU1zOvoy3kzc+zj54OXnh7eiNl5MXHg4eONd1tliViKJAVpa6\nJndqKly8qDZcJyVBYiIkJEB8vPpnzZrg7n578/C4vbVsaZmeT4Xl5udyNv0sUZeiOHnpJMdTjnMs\n5RiJ1xLxbehLx6YdDZt/I/9qWUI0BUkekjyqPUVRH6AJCbe3xET1YVqwpaaqpZf69dUk4uSkbo6O\namN+wWZnpzbk29mpxyo1M0nNjybh5ikSbsRwPjOGv65GE3f1L/L0ebRyaEUrh1a4N3DHvYE7Lexb\n0Ny+Oc3sm+FUx+mO5KLXq9V1OTlqddHNm2q31hs31Ebn69fVhJCZCdeuqdvVq2ojdXq6+mdamvo5\nLl9W2xlcXG5vrq7QtKmaLN3c1Kq+Zs0s09Ppbjn5OVy4eoG/Mv7ir4y/OJt+lui0aGLSY4i7Eoeb\nnRu+DX3xbehL20ZtCWwciLeTtyQKE5LkIclDGEmvv/MBnJamPpSvXoUrV9TXMjNvbwUP9Oxs9c9b\nt9Ttxg01CehrZmDV8BxWTnEoDWLBPg7FPg69bTx623iwuoUuqylcc4VMV5RrTbC66YL1zcbY3HKh\ntr4h9WhEXRpiV7cW9eqpCat+fbXKyN5eTWIODrc3Z+fbm6V6NBWmKAoZNzNIyUoh9XoqyVnJJF5L\nJCkriYRrCcRfjefC1QukXk/F1c6VVg6taNmgJV6OXmpJzskbDwcP6tiYaXpcYSDJQ5KHsJCCkkR+\nvvr3gk2nU7esnCxSbySRkp3IxexEUrMvkpqdwsWsi1zMusil65e4lH2JS9cvUdOqJo51HHGq60SD\n2g0Mm30te2xr2WJbU93q1axHXZu61LOpR23r2tS2rk0t61rUsqqFjZUNNjVssLGywUpnRQ1dDWro\nahh+RxQUFEUhT59HvpJPnj6PnPwcbuXd4lb+LW7l3SI7N5sbeTfIzs0m81YmmTnqdvXmVa7cvMKV\nm1dIv5FO2o000rLTSL+RTl2burjUd6FRvUY0qd8EVztXmtZvSlPbprRooJbEmto2lVJEBSPJQ5KH\nqOQURSErJ8vwUC54SF+5eYWrN68aHuDXbl0jOzf7jq3gwX8z7yZ5+jxy83PJ1eeiV/R3bDp06HQ6\ndOiwrmGNdQ1rrGpYUdOqJrWsahkSUL2a9ahjXYe6NnUNSat+zfrY17K/I6k513XGqa4TTnWcquzE\ngVWdJA9JHkIIUWplfXbKNKVCCCFKTZKHEEKIUpPkIYQQotQkeQghhCg1SR5CCCFKTZKHEEKIUpPk\nIYQQotQkeQghhCg1SR5CCCFKTZKHEEKIUrNY8nB3d6dt27YEBQXRuXNnANLT0wkJCcHb25t+/fpx\n5aTgFAAAAAguSURBVMoVw/EzZ87Ey8sLHx8ftm7daqmwK4Xw8HBLh1BhyL24Te7FbXIvys9iyUOn\n0xEeHs6RI0c4cOAAALNmzSIkJITo6Gj69OnDrFmzAIiKimLFihVERUWxefNmJkyYgF6vt1ToFZ78\nYtwm9+I2uRe3yb0oP4tWW909Gdf69esZNWoUAKNGjWLt2rUArFu3jhEjRmBjY4O7uzuenp6GhCOE\nEML8LFry6Nu3Lx07duTbb78FICUlBRcXFwBcXFxISUkBICkpCTc3N8N73dzcSExMNH/QQgghALDY\nqix79+6lSZMmXLp0iZCQEHx8fO54XafT3Xdd6KJes9Q60hXR9OnTLR1ChSH34ja5F7fJvSgfiyWP\nJk2aANCwYUOGDBnCgQMHcHFx4eLFizRu3Jjk5GQaNWoEgKurK/Hx8Yb3JiQk4Orqesf5ZC0PIYQw\nH4tUW2VnZ5OZmQnA9evX2bp1KwEBAYSFhbFo0SIAFi1axODBgwEICwtj+fLl5OTkEBsbS0xMjKGH\nlhBCCPOzSMkjJSWFIUOGAJCXl8fTTz9Nv3796NixI8OHD2fBggW4u7uzcuVKAHx9fRk+fDi+vr5Y\nW1szb948qaISQghLUiqZTZs2Ka1bt1Y8PT2VWbNmFXnMpEmTFE9PT6Vt27bK4cOHzRyh+ZR0L378\n8Uelbdu2SkBAgPLAAw8ox44ds0CU5mHM/wtFUZQDBw4oVlZWyurVq80YnXkZcy927typtGvXTvHz\n81N69epl3gDNqKR7cenSJaV///5KYGCg4ufnp3z//ffmD9IMnnvuOaVRo0aKv79/sceU9rlZqZJH\nXl6e4uHhocTGxio5OTlKYGCgEhUVdccxGzZsUAYOHKgoiqL88ccfSpcuXSwRqskZcy/27dunXLly\nRVEU9ZeoOt+LguN69+6tPPzww8pPP/1kgUhNz5h7kZGRofj6+irx8fGKoqgP0KrImHvx7rvvKq+/\n/rqiKOp9cHR0VHJzcy0Rrknt3r1bOXz4cLHJoyzPzUo1PcmBAwfw9PTE3d0dGxsbnnzySdatW3fH\nMYXHinTp0oUrV64YuvxWJcbci27dumFvbw+o9yIhIcESoZqcMfcC4PPPP+exxx6jYcOGFojSPIy5\nF0uXLmXYsGGG7u/Ozs6WCNXkjLkXTZo04dq1awBcu3YNJycnrK0t1o/IZHr27ImDg0Oxr5fluVmp\nkkdiYiLNmjUz/FzUeI+ijqmKD01j7kVhCxYsIDQ01ByhmZ2x/y/WrVvH+PHjgarbrduYexETE0N6\nejq9e/emY8eOLF682NxhmoUx92LcuHGcPHmSpk2bEhgYyGeffWbuMCuEsjw3K1WKNfYXXrmr225V\nfFCU5jPt3LmThQsXsnfvXhNGZDnG3ItXXnmFWbNmodPpUNTqWjNEZn7G3Ivc3FwOHz7Mb7/9RnZ2\nNt26daNr1654eXmZIULzMeZezJgxg3bt2hEeHs65c+cICQnh2LFj2NramiHCiqW0z81KlTzuHu8R\nHx9/x8jzoo4pakxIVWDMvQA4fvw448aNY/PmzfcttlZmxtyLyMhInnzySQAuX77Mpk2bsLGxISws\nzKyxmpox96JZs2Y4OztTp04d6tSpw4MPPsixY8eqXPIw5l7s27ePt956CwAPDw9atmzJmTNn6Nix\no1ljtbQyPTc1a5Exg9zcXKVVq1ZKbGyscuvWrRIbzPfv319lG4mNuRfnz59XPDw8lP3791soSvMw\n5l4UNnr06Crb28qYe3Hq1CmlT58+Sl5ennL9+nXF399fOXnypIUiNh1j7sWUKVOUadOmKYqiKBcv\nXlRcXV2VtLQ0S4RrcrGxsUY1mBv73KxUJQ9ra2u++OIL+vfvT35+PmPGjKFNmzZ8/fXXALz44ouE\nhoayceNGPD09qVevHt9//72FozYNY+7Ff/7zHzIyMgz1/DY2NlVyQklj7kV1Ycy98PHxYcCAAbRt\n25YaNWowbtw4fH19LRy59oy5F2+++SbPPfccgYGB6PV6PvroIxwdHS0cufZGjBjBrl27uHz5Ms2a\nNWP69Onk5uYCZX9u6hSlilb+CiGEMJlK1dtKCCFExSDJQwghRKlJ8hBCCFFqkjyEEEKUmiQPIYQQ\npSbJQwghRKlJ8hDCBFJSUmjRogUbN268Y39mZiY+Pj6sWbPGQpEJoQ1JHkKYwK+//kp8fLxhKeUC\n27dvJzo6uspOFSOqD0keQpjAjh07sLOzo0OHDnfs37VrFzVr1qRbt24WikwIbUjyEMIEdu7cSc+e\nPe+ZmTQ8PJxOnTpRq1YtC0UmhDYkeQihsVOnTnHx4kWCg4Pv2J+RkcGJEyd48MEHLROYEBqS5CGE\nxnbs2AFwT/LYs2cPer1ekoeoEiR5CKGxHTt2YGVlRWBg4B37w8PDsbKyonv37ly9epUffvjBMgEK\noQFJHkJoSFEUdu3aRcP/b+/uUSOEwigMn1axndbKzmWIVroDCxF7XZvuQDutXIXDMExjI4iQIiSQ\n8maUm8D7rOB0L/7Ad7v9uIW9bZvatlUYhvI8T33fy3Vdi0uB9xAP4ETzPOv1eun5fGpZFkmfZ1/r\nupbned+/7nZdpyzLbE4F3vKvjkEBf93X946qqpTnuYIg0L7vappGdV2rLEsVRaEoiuQ4juW1wO9x\nDAo4UZqmGsdRj8fD9hTgUry2Ak5yHIeGYVAcx7anAJcjHsBJpmnSuq5KksT2FOByxAM4yf1+l+/7\nStPU9hTgcnzzAAAY48kDAGCMeAAAjBEPAIAx4gEAMEY8AADGiAcAwBjxAAAYIx4AAGMfx+sBYVdn\nda0AAAAASUVORK5CYII=\n"
      }
     ],
     "prompt_number": 67
    },
    {
     "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",
     "collapsed": false,
     "input": [
      "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)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "latex": [
        "Fair Coin: $\\mu_{ML} = 0.493$"
       ],
       "output_type": "display_data",
       "text": [
        "<IPython.core.display.Latex at 0x5df4b10>"
       ]
      },
      {
       "latex": [
        "Biased Coin: $\\mu_{ML} = 0.602$"
       ],
       "output_type": "display_data",
       "text": [
        "<IPython.core.display.Latex at 0x5ddbdd0>"
       ]
      }
     ],
     "prompt_number": 68
    },
    {
     "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",
     "collapsed": false,
     "input": [
      "@vectorize\n",
      "def C(N,k):\n",
      "    return math.factorial(N)/(math.factorial(N-k)*math.factorial(k))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 69
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "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)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "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"
       ]
      }
     ],
     "prompt_number": 70
    },
    {
     "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",
     "collapsed": false,
     "input": [
      "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)\")"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "pyout",
       "prompt_number": 80,
       "text": [
        "<matplotlib.text.Text at 0x748b910>"
       ]
      },
      {
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAXUAAAEMCAYAAAA70CbBAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAHeJJREFUeJzt3X9QlHUeB/D3g9A5JCpesji75noCwvJrVzmgKbo1QfNU\nBrUorNwU606vbjybqXPm6qBphK76w+y8a26scPpB6iVYlx5ed5tm460pnHdDBhkkLAvGIaIiLbDf\n+8NrE/m1C7s+7Jf3a4aZ9dnn+3w/C/jeZ798n++jCCEEiIhICkFqF0BERL7DUCcikghDnYhIIgx1\nIiKJMNSJiCTCUCciksiQod7V1YW0tDQYjUYYDAZs2bIFAFBQUACdTgeTyQSTyYQDBw642xQVFSE6\nOhqxsbGoqKjwb/VERNSHMtw89c7OToSGhqKnpwd33HEHXnzxRXz00UcICwvD5s2b++xbXV2N1atX\n4/jx47Db7cjMzERNTQ2CgviBgIjoRhg2bUNDQwEATqcTvb29CA8PBwAM9F5QXl6OvLw8hISEQK/X\nIyoqCjabzcclExHRYIYNdZfLBaPRCI1GgwULFiA+Ph4AsH37diQnJyM/Px/t7e0AgKamJuh0Ondb\nnU4Hu93up9KJiOh6wcPtEBQUhKqqKly4cAGLFy+G1WrFhg0b8MwzzwAAnn76aTzxxBPYuXPngO0V\nRfFoGxERDW+4lV08HuyeMmUKli5dis8++wwRERFQFAWKomD9+vXuIRatVouGhgZ3m8bGRmi12kEL\nC9Sv3/72t6rXMB5rZ/3qf7F+db88MWSot7a2uodWrly5gkOHDsFkMqG5udm9z759+5CYmAgAyM7O\nRmlpKZxOJ+rq6lBbW4vU1FSPCiEiotEbcvjF4XDAYrHA5XLB5XLhoYcewsKFC7FmzRpUVVVBURTM\nnj0br776KgDAYDAgNzcXBoMBwcHB2LFjB4daiIhuoGGnNPqlU0Xx+KPEWGS1WmE2m9UuY0QCuXaA\n9auN9avLk+xkqBMRBQhPspNXBRERSWTYKY1E493kydNw8eJ5tcvoJywsHB0dbWqXQWMMh1+IhnH1\nj/1j8feV/4/GGw6/EBGNMwx1IiKJMNSJiCTCUCcikghDnYhIIgx1IiKJMNSJiCTCUCcikghDnYhI\nIgx1IiKJMNSJiCTCUCcikghDnYhIIgx1IiKJMNSJiCTCUCcikghDnYhIIgx1IiKJMNSJiCQyZKh3\ndXUhLS0NRqMRBoMBW7ZsAQC0tbUhKysLMTExWLRoEdrb291tioqKEB0djdjYWFRUVPi3eiIi6mPY\nG093dnYiNDQUPT09uOOOO/Diiy9i//79uOWWW/Dkk0/i+eefx/nz51FcXIzq6mqsXr0ax48fh91u\nR2ZmJmpqahAU1Pe9gzeepkDCG0/TWOGTG0+HhoYCAJxOJ3p7exEeHo79+/fDYrEAACwWC8rKygAA\n5eXlyMvLQ0hICPR6PaKiomCz2Ub7OoiIyEPBw+3gcrkwb948nDlzBhs2bEB8fDxaWlqg0WgAABqN\nBi0tLQCApqYmpKenu9vqdDrY7fYBj1tQUOB+bDabYTabR/EyiIjkY7VaYbVavWozbKgHBQWhqqoK\nFy5cwOLFi/GPf/yjz/OKovz/4+nABnvu2lAnIqL+rj/hLSwsHLaNx7NfpkyZgqVLl+LEiRPQaDRo\nbm4GADgcDkRERAAAtFotGhoa3G0aGxuh1Wo97YKIiEZpyFBvbW11z2y5cuUKDh06BJPJhOzsbJSU\nlAAASkpKkJOTAwDIzs5GaWkpnE4n6urqUFtbi9TUVD+/BCIi+s6Qwy8OhwMWiwUulwsulwsPPfQQ\nFi5cCJPJhNzcXOzcuRN6vR67d+8GABgMBuTm5sJgMCA4OBg7duwYcmiGiIh8a9gpjX7plFMaKYBw\nSiONFT6Z0khERIGDoU5EJBGGOhGRRBjqREQSYagTEUmEoU5EJBGGOhGRRIZd+4WIAtvkydNw8eJ5\ntcvoJywsHB0dbWqXIR1efEQ0jEC/+CjQ66fv8eIjIqJxhqFORCQRhjoRkUQY6kREEmGoExFJhKFO\nRCQRhjoRkUQY6kREEmGoExFJhKFORCQRhjoRkUQY6kREEmGoExFJhKFORCSRIUO9oaEBCxYsQHx8\nPBISEvDyyy8DAAoKCqDT6WAymWAymXDgwAF3m6KiIkRHRyM2NhYVFRX+rZ6IiPoYcj315uZmNDc3\nw2g04tKlS5g/fz7Kysqwe/duhIWFYfPmzX32r66uxurVq3H8+HHY7XZkZmaipqYGQUF93zu4njoF\nkkBfjzzQ66fvjXo99cjISBiNRgDApEmTEBcXB7vdDgADHri8vBx5eXkICQmBXq9HVFQUbDbbSOsn\nIiIveXw7u/r6elRWViI9PR1Hjx7F9u3bsWvXLqSkpOCll17C1KlT0dTUhPT0dHcbnU7nfhO4XkFB\ngfux2WyG2Wwe8YsgIpKR1WqF1Wr1qo1Ht7O7dOkSzGYzfvOb3yAnJwfnzp3D9OnTAQBPP/00HA4H\ndu7ciccffxzp6el44IEHAADr16/HT3/6U6xcubJvpxx+oQAS6MMXgV4/fc8nt7Pr7u7GqlWr8OCD\nDyInJwcAEBERAUVRoCgK1q9f7x5i0Wq1aGhocLdtbGyEVqsdzWsgIiIvDBnqQgjk5+fDYDBg06ZN\n7u0Oh8P9eN++fUhMTAQAZGdno7S0FE6nE3V1daitrUVqaqqfSiciousNOaZ+9OhRvPnmm0hKSoLJ\nZAIAbN26Fe+88w6qqqqgKApmz56NV199FQBgMBiQm5sLg8GA4OBg7Nix4/8f/YiI6EbwaEzd551y\nTJ0CSKCPSQd6/fQ9n4ypExFR4GCoExFJhKFORCQRhjoRkUQY6kREEmGoExFJhKFORCQRhjoRkUQY\n6kREEmGoExFJhKFORCQRhjoRkUQY6kREEmGoExFJhKFORCQRj288TTRSkydPw8WL59Uuo5+wsHB0\ndLSpXQaRT/EmGeR3gX6TBtbvL8wBb/EmGURE4wxDnYhIIgx1IiKJMNSJiCTCUCcikghDnYhIIkOG\nekNDAxYsWID4+HgkJCTg5ZdfBgC0tbUhKysLMTExWLRoEdrb291tioqKEB0djdjYWFRUVPi3eiIi\n6mPIeerNzc1obm6G0WjEpUuXMH/+fJSVleH111/HLbfcgieffBLPP/88zp8/j+LiYlRXV2P16tU4\nfvw47HY7MjMzUVNTg6Cgvu8dnKc+vgT6PGnW7y/MAW+Nep56ZGQkjEYjAGDSpEmIi4uD3W7H/v37\nYbFYAAAWiwVlZWUAgPLycuTl5SEkJAR6vR5RUVGw2Wy+eC1EROQBj5cJqK+vR2VlJdLS0tDS0gKN\nRgMA0Gg0aGlpAQA0NTUhPT3d3Uan08Futw94vIKCAvdjs9kMs9k8gvKJiORltVphtVq9auNRqF+6\ndAmrVq3Ctm3bEBYW1uc5RVH+//FuYIM9d22oExFRf9ef8BYWFg7bZtjZL93d3Vi1ahUeeugh5OTk\nALh6dt7c3AwAcDgciIiIAABotVo0NDS42zY2NkKr1Xr1IoiIaOSGDHUhBPLz82EwGLBp0yb39uzs\nbJSUlAAASkpK3GGfnZ2N0tJSOJ1O1NXVoba2FqmpqX4sn4iIrjXk7JdPPvkEd955J5KSktzDKEVF\nRUhNTUVubi7Onj0LvV6P3bt3Y+rUqQCArVu34rXXXkNwcDC2bduGxYsX9++Us1/GlUCffcH6/YU5\n4C1PspNL75LfBXqosH5/YQ54i0vvEhGNMwx1IiKJMNSJiCTCUCcikghDnYhIIgx1IiKJMNSJiCTC\nUCcikghDnYhIIgx1IiKJMNSJiCTCUCcikghDnYhIIgx1IiKJMNSJiCTCUCcikghDnYhIIgx1IiKJ\nMNSJiCTCUCcikghDnYhIIgx1IiKJMNSJiCQybKivW7cOGo0GiYmJ7m0FBQXQ6XQwmUwwmUw4cOCA\n+7mioiJER0cjNjYWFRUV/qmaiIgGpAghxFA7HDlyBJMmTcKaNWvw73//GwBQWFiIsLAwbN68uc++\n1dXVWL16NY4fPw673Y7MzEzU1NQgKKjve4eiKBimW5KIoigAxuLP27PfQ9bvL8wBb3mSncOeqWdk\nZCA8PLzf9oEOXF5ejry8PISEhECv1yMqKgo2m82LkomIaDSCR9pw+/bt2LVrF1JSUvDSSy9h6tSp\naGpqQnp6unsfnU4Hu90+YPuCggL3Y7PZDLPZPNJSiIikZLVaYbVavWozolDfsGEDnnnmGQDA008/\njSeeeAI7d+4ccN+rH/36uzbUiYiov+tPeAsLC4dtM6LZLxEREVAUBYqiYP369e4hFq1Wi4aGBvd+\njY2N0Gq1I+mCiIhGYESh7nA43I/37dvnnhmTnZ2N0tJSOJ1O1NXVoba2Fqmpqb6plIiIhjXs8Ete\nXh4+/vhjtLa2YubMmSgsLITVakVVVRUURcHs2bPx6quvAgAMBgNyc3NhMBgQHByMHTt2DDr8QkRE\nvjfslEa/dMopjeNKoE+pY/3+whzwlk+mNBIRUeBgqBMRSYShTkQkEYY6EZFEGOpERBJhqBMRSYSh\nTkQkEYY6EZFEGOpERBJhqBMRSYShTkQkEYY6EZFEGOpERBJhqBMRSYShTkQkEYY6EZFEGOpERBJh\nqBMRSYShTkQkEYY6EZFEGOpERBJhqBMRSYShTkQkkWFDfd26ddBoNEhMTHRva2trQ1ZWFmJiYrBo\n0SK0t7e7nysqKkJ0dDRiY2NRUVHhn6qJiGhAw4b62rVrcfDgwT7biouLkZWVhZqaGixcuBDFxcUA\ngOrqarz77ruorq7GwYMHsXHjRrhcLv9UTkRE/Qwb6hkZGQgPD++zbf/+/bBYLAAAi8WCsrIyAEB5\neTny8vIQEhICvV6PqKgo2Gw2P5RNREQDCR5Jo5aWFmg0GgCARqNBS0sLAKCpqQnp6enu/XQ6Hex2\n+4DHKCgocD82m80wm80jKYWISFpWqxVWq9WrNiMK9WspigJFUYZ8fiDXhjoREfV3/QlvYWHhsG1G\nNPtFo9GgubkZAOBwOBAREQEA0Gq1aGhocO/X2NgIrVY7ki6IiGgERhTq2dnZKCkpAQCUlJQgJyfH\nvb20tBROpxN1dXWora1Famqq76olIqIhDTv8kpeXh48//hitra2YOXMmnn32Wfz6179Gbm4udu7c\nCb1ej927dwMADAYDcnNzYTAYEBwcjB07dgw5NENERL6lCCHEDe9UUaBCt6SSq2/sY/Hn7dnvIev3\nF+aAtzzJTl5RSkQkEYY6EZFEGOpERBJhqBMRSYShTkQkkVFfUUr+N3nyNFy8eF7tMgYUFhaOjo42\ntcsgov/jlMYAMHanpAGeTEsbu/WPjymBgV4/fY9TGomIxhmGOhGRRBjqREQSYagTEUmEs1+IaEwb\nq7O/xurML85+CQBjd/YCwNkvamL96rrxOcbZL0RE4wxDnYhIIgx1IiKJMNSJiCTCUCcikghDnYhI\nIgx1IiKJMNSJiCTCUCcikghDnYhIIqNa+0Wv12Py5MmYMGECQkJCYLPZ0NbWhvvuuw9ff/019Ho9\ndu/ejalTp/qqXiIiGsKoztQVRYHVakVlZSVsNhsAoLi4GFlZWaipqcHChQtRXFzsk0KJiGh4ox5+\nuX5xmf3798NisQAALBYLysrKRtsFERF5aFTDL4qiIDMzExMmTMDPfvYzPPLII2hpaYFGowEAaDQa\ntLS0DNi2oKDA/dhsNsNsNo+mFCIi6VitVlitVq/ajGrpXYfDgRkzZuCbb75BVlYWtm/fjuzsbJw/\n//3ax9OmTUNbW981h7n0rnfG7tKjAJfeVRPrV5eES+/OmDEDADB9+nSsWLECNpsNGo0Gzc3NAK6G\nfkRExGi6ICIiL4w41Ds7O3Hx4kUAwOXLl1FRUYHExERkZ2ejpKQEAFBSUoKcnBzfVEpERMMa8fBL\nXV0dVqxYAQDo6enBAw88gC1btqCtrQ25ubk4e/bsoFMaOfzinbH78RPg8IuaWL+6xubwC29nFwDG\n7i81wFBXE+tX19gMdV5RSkQkEYY6EZFEGOpERBJhqBMRSYShTkQkEYY6EZFEGOpERBJhqBMRSYSh\nTkQkEYY6EZFEGOpERBJhqBMRSYShTkQkEYY6EZFEGOpERBJhqBMRSYShTkQkEYY6EZFEGOpERBJh\nqBMRSSRY7QJuhMbGRnR0dKhdRj9BQUGIiYlBUBDfW4nINxRxo2+HDc/uiO1LN988BRMmzMBY+2DS\n1XUW77//ZyxevHjI/cbu3dQBT+6oPnbrHx93s2f9/nJjcwzwLDv9cqZ+8OBBbNq0Cb29vVi/fj2e\neuopf3Tjse7ubnR2ngQQ6qMjWgGYR32UKVOWw+l0jvo43rHCF7WrxwrWryYrWP/Y5vNT197eXjz2\n2GM4ePAgqqur8c477+Dzzz/3dTcqs6pdwChY1S5glKxqFzBKVrULGCWr2gWMklXtAvzO56Fus9kQ\nFRUFvV6PkJAQ3H///SgvL/d1N0RENACfD7/Y7XbMnDnT/W+dTod//vOfvu7GK0FBQZg8eRUUxTcv\nt6urBhMnnhj1cb791oagoJ/7oCIioqt8HupX/6jhu/185dtvD/r4eDU+Oc6yZcs83NOX369CHx7L\n05/l2Kzf899D1v891v+dG51jnvB5qGu1WjQ0NLj/3dDQAJ1O12cfFSbcEBGNCz4fU09JSUFtbS3q\n6+vhdDrx7rvvIjs729fdEBHRAHx+ph4cHIxXXnkFixcvRm9vL/Lz8xEXF+frboiIaAB+uRpnyZIl\n+OKLL/Dll19iy5Yt/Z7fs2cP4uPjMWHCBJw8edIfJfjFwYMHERsbi+joaDz//PNql+OVdevWQaPR\nIDExUe1SRqShoQELFixAfHw8EhIS8PLLL6tdkle6urqQlpYGo9EIg8Ew4P+Lsa63txcmkwnLly9X\nuxSv6fV6JCUlwWQyITU1Ve1yvNbe3o577rkHcXFxMBgMOHbs2OA7CxV8/vnn4osvvhBms1mcOHFC\njRK81tPTI+bMmSPq6uqE0+kUycnJorq6Wu2yPHb48GFx8uRJkZCQoHYpI+JwOERlZaUQQoiLFy+K\nmJiYgPr+CyHE5cuXhRBCdHd3i7S0NHHkyBGVK/LOSy+9JFavXi2WL1+udile0+v14r///a/aZYzY\nmjVrxM6dO4UQV39/2tvbB91XlevmY2NjERMTo0bXIxbo8+8zMjIQHh6udhkjFhkZCaPRCACYNGkS\n4uLi0NTUpHJV3gkNvXpFs9PpRG9vL6ZNm6ZyRZ5rbGzEhx9+iPXr1wfsRIdArfvChQs4cuQI1q1b\nB+DqEPeUKVMG3X9sLYYyhg00/95ut6tY0fhVX1+PyspKpKWlqV2KV1wuF4xGIzQaDRYsWACDwaB2\nSR771a9+hRdeeCFgF59TFAWZmZlISUnBn/70J7XL8UpdXR2mT5+OtWvXYt68eXjkkUfQ2dk56P5+\n+wllZWUhMTGx39f777/vry79aizORx2PLl26hHvuuQfbtm3DpEmT1C7HK0FBQaiqqkJjYyMOHz4M\nq9Wqdkke+eCDDxAREQGTyRSwZ7tHjx5FZWUlDhw4gN///vc4cuSI2iV5rKenBydPnsTGjRtx8uRJ\n3HzzzSguLh50f78tvXvo0CF/HVoVnsy/J//q7u7GqlWr8OCDDyInJ0ftckZsypQpWLp0KT777DOY\nzWa1yxnWp59+iv379+PDDz9EV1cXOjo6sGbNGuzatUvt0jw2Y8YMAMD06dOxYsUK2Gw2ZGRkqFyV\nZ3Q6HXQ6HX784x8DAO65554hQ131z1KB8s7P+ffqEkIgPz8fBoMBmzZtUrscr7W2tqK9vR0AcOXK\nFRw6dAgmk0nlqjyzdetWNDQ0oK6uDqWlpbjrrrsCKtA7Oztx8eJFAMDly5dRUVERULPAIiMjMXPm\nTNTUXL2K/W9/+xvi4+MHb3AD/nDbz3vvvSd0Op2YOHGi0Gg04u6771ajDK99+OGHIiYmRsyZM0ds\n3bpV7XK8cv/994sZM2aIm266Seh0OvHaa6+pXZJXjhw5IhRFEcnJycJoNAqj0SgOHDigdlkeO3Xq\nlDCZTCI5OVkkJiaK3/3ud2qXNCJWqzXgZr989dVXIjk5WSQnJ4v4+PiA+78rhBBVVVUiJSVFJCUl\niRUrVgw5+0WVm2QQEZF/qD78QkREvsNQJyKSCEOdiEgiDHUiIokw1Mkn6uvrPZomdu7cOSxdutTv\n9dyoC5Py8vKQnJyMbdu29dn+8MMP489//rPP+3vjjTfw+OOPD/r8qVOnkJ+f7/N+KXD47eIjooG8\n8sorePjhh/3ez2iuAO7t7cWECROG3a+5uRmfffYZamtrfdr/aCQlJeHMmTM4d+4cIiIiVKmB1MUz\ndfK5r776CvPmzcOJE/3v47p37173mfobb7yBlStXYsmSJYiJicFTTz3l3u/aM+29e/di7dq1AK6e\nAW/cuBG33XYb5syZA6vVCovFAoPB4N7nO5s3b0ZCQgIyMzPR2toKADhz5gyWLFmClJQU3Hnnnfji\niy/cx/35z3+O9PT0PnUAV5fNXbt2LZKSkjBv3jz35f2LFi2C3W6HyWTCJ5980u+1Hj58GLfffjvm\nzJnT56z9hRdeQGpqKpKTk1FQUODevmLFCqSkpCAhIaHP+iSvv/465s6di7S0NHz66afu7Xv27EFi\nYiKMRiN+8pOfuLcvWbIEe/bs6VcPjRM3bPY8Sa2urk4kJCSI06dPC5PJJE6dOtVvH4fD0Wfp39df\nf1386Ec/Eh0dHaKrq0vMmjVLNDY2CiGEmDRpknu/vXv3iocfflgIIYTFYhF5eXlCCCHKy8tFWFiY\n+M9//iNcLpeYP3+++Ne//iWEEEJRFPH2228LIYR49tlnxWOPPSaEEOKuu+4StbW1Qgghjh07Ju66\n6y73cZcvXy5cLle/ul988UWRn58vhBDi9OnT4tZbbxXffvutqK+vH3QpY4vFInJzc4UQQlRXV4uo\nqCghhBB//etfxaOPPiqEEKK3t1csW7ZMHD58WAghRFtbmxBCiM7OTpGQkCDa2tpEU1OTuPXWW0Vr\na6twOp3i9ttvF48//rgQQojExETR1NQkhBDiwoUL7r7//ve/u/um8YfDL+Qz586dQ05ODvbt24fY\n2Nh+z3/99dfuNTiAq0MUCxcuRFhYGADAYDDg66+/hlarHbQPRVHcN2lISEhAZGSk+5Lp+Ph41NfX\nIykpCUFBQbjvvvsAAA8++CBWrlyJy5cv49NPP8W9997rPp7T6XQf99577x1w2OTo0aP45S9/CQCY\nO3cuZs2ahZqamiHH7RVFca9PExcXh5aWFgBARUUFKioq3EsEXL58GV9++SUyMjKwbds2lJWVAbi6\n1G1NTQ0cDgfMZjN++MMfAgDuu+8+9+Xit99+OywWC3Jzc7Fy5Up33zNmzEB9ff2gtZHcGOrkM1On\nTsWsWbNw5MiRAUMd6L/Wzw9+8AP34wkTJqCnpwdA3zHpK1eu9Glz0003Abi66uG17YOCgtztr+9T\nURS4XC6Eh4ejsrJywNq+W+/ck7o98V2d17ffsmULHn300T77Wq1WfPTRRzh27BgmTpyIBQsWoKur\nq9+bzLXH+cMf/gCbzYa//OUvmD9/Pk6cOIFp06a5Xy+NTxxTJ5+56aab8N5772HXrl145513+j0/\na9YsNDc3u/89VFBqNBqcPn0aLpcL+/bt8zqkXC6Xe1z57bffRkZGBsLCwjB79mzs3bvX3f+pU6eG\nPVZGRgbeeustAEBNTQ3Onj2LuXPnelXPdxYvXozXXnsNly9fBnB1nf5vvvkGHR0dCA8Px8SJE3H6\n9GkcO3YMiqIgLS0NH3/8Mdra2tDd3Y09e/a4vxdnzpxBamoqCgsLMX36dDQ2NgIAHA4HZs2aNaL6\nKPDxTJ18RlEUhIaG4oMPPkBWVhbCwsKwbNky9/ORkZHo6elBZ2cnQkNDoSjKoGFdXFyMZcuWYfr0\n6UhJSXGH4Hf9DPT4WjfffDNsNhuee+45aDQavPvuuwCAt956Cxs2bMBzzz2H7u5u5OXlISkpachj\nbdy4ERs2bEBSUhKCg4NRUlKCkJCQIdsMVmdWVhY+//xz3HbbbQCAsLAwvPnmm7j77rvxxz/+EQaD\nAXPnznU/HxkZiYKCAtx2222YOnVqn5Udn3zySdTW1kIIgczMTPfrsNlsuPPOOweti+TGBb3ohioo\nKEBcXJx7vJt8z2w2Y/fu3ZzSOE5x+IVuqF/84hcoKSlRuwxpnTp1ClFRUQz0cYxn6kREEuGZOhGR\nRBjqREQSYagTEUmEoU5EJBGGOhGRRBjqREQS+R/FJKkarMVXywAAAABJRU5ErkJggg==\n"
      }
     ],
     "prompt_number": 80
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Compare with the actual distribution:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "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\")"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "pyout",
       "prompt_number": 81,
       "text": [
        "<matplotlib.text.Text at 0x76eed10>"
       ]
      },
      {
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAXgAAAEKCAYAAAAYd05sAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAF6lJREFUeJzt3X9M1Pfhx/HXKbc0FX9hqaV3l97Wo97ZAdIeEteQ3DqR\naO0FaJPhXNt0jBEX02zZP9v+GfSPpWTxD1f2B2lsM2NjWdJVXKeEknrR0cFl6mITp8VO5kGlda1O\nhC7I7fP9w/a+RfRzd3pw3NvnIyHhw+f9vs/r8MPLDx8+nzuHZVmWAADGWZDtAACA2UHBA4ChKHgA\nMBQFDwCGouABwFAUPAAYKmnBd3d3y+/3q7i4WG1tbTPWd3V1qaysTOXl5Xr00Uf17rvvJtZ5vV6V\nlpaqvLxca9euzWxyAIAth9118PF4XKtWrVJvb69cLpcqKiq0d+9eBQKBxJjx8XEtWrRIkvT++++r\nrq5OZ86ckSR9/etf19GjR1VQUDDLTwMAcD3bI/hoNCqfzyev1yun06mGhgZ1dXVNG/NluUvSlStX\ndM8990xbz31UAJAdtgU/MjIij8eTWHa73RoZGZkxbt++fQoEAtq4caN++9vfJr7ucDi0fv16BYNB\nvfLKKxmMDQBIJs9upcPhSOlBamtrVVtbqyNHjuiZZ57R6dOnJUl9fX0qKirShQsXVF1dLb/fr6qq\nqlvaBgBgumRnSGyP4F0ul2KxWGI5FovJ7XbfdHxVVZWmpqb06aefSpKKiookSYWFhaqrq1M0Gr1p\nyFz9+NWvfpX1DOTPfo47MX8uZzchfypsCz4YDGpwcFBDQ0OanJxUZ2enwuHwtDEffvhhYmPHjh2T\nJK1YsUITExMaGxuTdO0PsT09PSopKUkpFDBfLVlSIIfDkbGP1tbWjD3WkiVczIDpbE/R5OXlqb29\nXTU1NYrH42psbFQgEFBHR4ckqbm5WW+++aZ2794tp9Op/Px8vfHGG5Kk0dFR1dfXS5Kmpqa0detW\nbdiwYZafDjC7xsYuSsrkhQMtX3zcvrExTndiOtvLJOckgMOR8q8b81EkElEoFMp2jFtG/vRc+5tR\nJvfXiKRQhh5rbn+W2HeyK5XupOCBNGS+4DOJn6U7SSrdyUsVAIChKHgAMBQFDwCGouABwFAUPAAY\nioIHAENR8ABgKAoeAAxFwQOAoSh4ADAUBQ8AhqLgAcBQFDwAGIqCBwBDUfAAYCgKHgAMRcEDgKEo\neAAwFAUPAIai4AHAUBQ8ABiKggcAQyUt+O7ubvn9fhUXF6utrW3G+q6uLpWVlam8vFyPPvqo3n33\n3ZTnAgBmj8OyLOtmK+PxuFatWqXe3l65XC5VVFRo7969CgQCiTHj4+NatGiRJOn9999XXV2dzpw5\nk9JcSXI4HLKJAMwrDodD0nzdX/lZupOk0p22R/DRaFQ+n09er1dOp1MNDQ3q6uqaNubLcpekK1eu\n6J577kl5LgBg9uTZrRwZGZHH40ksu91uDQwMzBi3b98+/eIXv9D58+fV09OT1lxJamlpSXweCoUU\nCoXSeQ4AYLxIJKJIJJLWHNuCv/braHK1tbWqra3VkSNH9Mwzz+jUqVNphfhqwQMAZrr+4Le1tTXp\nHNtTNC6XS7FYLLEci8XkdrtvOr6qqkpTU1P67LPP5Ha705oLAMgs24IPBoMaHBzU0NCQJicn1dnZ\nqXA4PG3Mhx9+mDjRf+zYMUnSihUrUpoLAJg9tqdo8vLy1N7erpqaGsXjcTU2NioQCKijo0OS1Nzc\nrDfffFO7d++W0+lUfn6+3njjDdu5AIC5YXuZ5JwE4DJJ5BAuk8R8cduXSQIAchcFDwCGouABwFAU\nPAAYioIHAENR8ABgKAoeAAxFwQOAoSh4ADAUBQ8AhqLgAcBQFDwAGIqCBwBDUfAAYCgKHgAMRcED\ngKEoeAAwFAUPAIai4AHAUBQ8ABiKggcAQ+VlOwCAubNkSYHGxi5mO8YMixcv1+XLn2U7hnEclmVZ\nWQ3gcCjLEYCUORwOSfN1f03+szR/89MD6UqlO5Oeounu7pbf71dxcbHa2tpmrH/99ddVVlam0tJS\nPfbYYzpx4kRindfrVWlpqcrLy7V27dpbeAoAgFtle4omHo9r+/bt6u3tlcvlUkVFhcLhsAKBQGLM\nN77xDR0+fFhLly5Vd3e3fvSjH6m/v1/Stf9hIpGICgoKZvdZAABmsD2Cj0aj8vl88nq9cjqdamho\nUFdX17Qx69at09KlSyVJlZWVGh4enraeX7sAIDtsj+BHRkbk8XgSy263WwMDAzcdv2vXLm3atCmx\n7HA4tH79ei1cuFDNzc1qamq64byWlpbE56FQSKFQKMX4AHBniEQiikQiac2xLfhrf5BJzaFDh/Tq\nq6+qr68v8bW+vj4VFRXpwoULqq6ult/vV1VV1Yy5Xy14AMBM1x/8tra2Jp1je4rG5XIpFosllmOx\nmNxu94xxJ06cUFNTk/bv36/ly5cnvl5UVCRJKiwsVF1dnaLRaNJAAIDMsC34YDCowcFBDQ0NaXJy\nUp2dnQqHw9PGnDt3TvX19dqzZ498Pl/i6xMTExobG5MkjY+Pq6enRyUlJbPwFAAAN2J7iiYvL0/t\n7e2qqalRPB5XY2OjAoGAOjo6JEnNzc168cUXdfHiRW3btk2S5HQ6FY1GNTo6qvr6eknS1NSUtm7d\nqg0bNszy0wEAfIkbnYA0zN8bhSRudLqzZORGJwBAbqLgAcBQFDwAGIqCBwBDUfAAYCgKHgAMRcED\ngKF4RyfMOd5VCJgb3OiEOZfLN9vM3+xSbuenB9LFjU4AcAej4AHAUBQ8ABiKggcAQ1HwAGAoCh4A\nDEXBA4ChKHgAMBQFDwCGouABwFAUPAAYioIHAENR8ABgKAoeAAyVtOC7u7vl9/tVXFystra2Getf\nf/11lZWVqbS0VI899phOnDiR8lwAwCyybExNTVkPPvigdfbsWWtyctIqKyuzTp48OW3Me++9Z126\ndMmyLMs6ePCgVVlZmfLcL16L3i4CDCTJkqx5+JF8X5y/2XM9Pz2QrlS+Z7ZH8NFoVD6fT16vV06n\nUw0NDerq6po2Zt26dVq6dKkkqbKyUsPDwynPBQDMHtu37BsZGZHH40ksu91uDQwM3HT8rl27tGnT\nprTntrS0JD4PhUIKhUKpZAeAO0YkElEkEklrjm3BX3t7r9QcOnRIr776qvr6+tKe+9WCBwDMdP3B\nb2tra9I5tgXvcrkUi8USy7FYTG63e8a4EydOqKmpSd3d3Vq+fHlacwEAs8P2HHwwGNTg4KCGhoY0\nOTmpzs5OhcPhaWPOnTun+vp67dmzRz6fL625AIDZY3sEn5eXp/b2dtXU1Cgej6uxsVGBQEAdHR2S\npObmZr344ou6ePGitm3bJklyOp2KRqM3nQsAmBuOLy63yV4Ah0NZjoA5du3vM/Px3zz5vjh/s0u5\nnZ8eSFcq3cmdrABgKAoeAAxFwQOAoSh4ADAUBQ8AhqLgAcBQFDwAGIqCBwBDUfAAYCgKHgAMRcED\ngKEoeAAwFAUPAIai4AHAUBQ8ABiKggcAQ1HwAGAoCh4ADEXBA4ChKHgAMBQFDwCGouABwFAUPAAY\nKmnBd3d3y+/3q7i4WG1tbTPWnzp1SuvWrdNdd92lHTt2TFvn9XpVWlqq8vJyrV27NnOpAQBJ5dmt\njMfj2r59u3p7e+VyuVRRUaFwOKxAIJAYs2LFCr388svat2/fjPkOh0ORSEQFBQWZTw4AsGV7BB+N\nRuXz+eT1euV0OtXQ0KCurq5pYwoLCxUMBuV0Om/4GJZlZS4tACBltkfwIyMj8ng8iWW3262BgYGU\nH9zhcGj9+vVauHChmpub1dTUdMNxLS0tic9DoZBCoVDK2wCAO0EkElEkEklrjm3BOxyO28mjvr4+\nFRUV6cKFC6qurpbf71dVVdWMcV8teADATNcf/La2tiadY3uKxuVyKRaLJZZjsZjcbnfKgYqKiiRd\nO41TV1enaDSa8lwAwO2xLfhgMKjBwUENDQ1pcnJSnZ2dCofDNxx7/bn2iYkJjY2NSZLGx8fV09Oj\nkpKSDMUGACRje4omLy9P7e3tqqmpUTweV2NjowKBgDo6OiRJzc3NGh0dVUVFhS5fvqwFCxZo586d\nOnnypD755BPV19dLkqamprR161Zt2LBh9p8RAECS5LCyfJmLw+HgSps7zLW/7czHf/Pk++L8zS7l\ndn56IF2pdCd3sgKAoSh4ADAUBQ8AhqLgAcBQFDwAGIqCBwBDUfAAYCgKHgAMRcEDgKEoeAAwFAUP\nAIai4AHAUBQ8ABiKggcAQ1HwAGAoCh4ADEXBA4ChKHgAMBQFDwCGouABwFAUPAAYioIHAENR8ABg\nqKQF393dLb/fr+LiYrW1tc1Yf+rUKa1bt0533XWXduzYkdZcAMDscViWZd1sZTwe16pVq9Tb2yuX\ny6WKigrt3btXgUAgMebChQv617/+pX379mn58uX62c9+lvJcSXI4HLKJAAM5HA5J8/HfPPm+OH+z\nS7mdnx5IVyrdaXsEH41G5fP55PV65XQ61dDQoK6urmljCgsLFQwG5XQ6054LAJg9eXYrR0ZG5PF4\nEstut1sDAwMpPXA6c1taWhKfh0IhhUKhlLYBAHeKSCSiSCSS1hzbgr/269ytSWfuVwseADDT9Qe/\nra2tSefYnqJxuVyKxWKJ5VgsJrfbnVKY25kLALh9tgUfDAY1ODiooaEhTU5OqrOzU+Fw+IZjrz/Z\nn85cAEDm2Z6iycvLU3t7u2pqahSPx9XY2KhAIKCOjg5JUnNzs0ZHR1VRUaHLly9rwYIF2rlzp06e\nPKn8/PwbzgUAzA3byyTnJACXSd5xcvlSvfmbXcrt/PRAum77MkkAQO6i4AHAUBQ8ABiKggcAQ1Hw\nAGAoCh4ADEXBA4ChKHgAMJTtnawAMJ8sWVKgsbGL2Y4xw+LFy3X58mfZjjEDd7LmoFzfyXP5bsr5\nm13K7fyp9UCu58/oFlPoTgo+B+X6Tp7L+edvdim385u/72R8i7xUAQDcuSh4ADAUBQ8AhqLgAcBQ\nFDwAGIqCBwBDUfAAYCgKHgAMRcEDgKEoeAAwFAUPAIai4AHAUEkLvru7W36/X8XFxWpra7vhmBde\neEHFxcUqKyvT8ePHE1/3er0qLS1VeXm51q5dm7nUAICkbF8PPh6Pa/v27ert7ZXL5VJFRYXC4bAC\ngUBizIEDB3TmzBkNDg5qYGBA27ZtU39/v6Rrr3YWiURUUFAwu88CADCD7RF8NBqVz+eT1+uV0+lU\nQ0ODurq6po3Zv3+/nnvuOUlSZWWlLl26pI8//jixnpcCBoDssD2CHxkZkcfjSSy73W4NDAwkHTMy\nMqKVK1fK4XBo/fr1WrhwoZqbm9XU1HTD7bS0tCQ+D4VCCoVCt/BUAMBckUhEkUgkrTm2BX/txfWT\nu9lR+l/+8hfdf//9unDhgqqrq+X3+1VVVTVj3FcLHgAw0/UHv62trUnn2J6icblcisViieVYLCa3\n2207Znh4WC6XS5J0//33S5IKCwtVV1enaDSa/FkAADLCtuCDwaAGBwc1NDSkyclJdXZ2KhwOTxsT\nDoe1e/duSVJ/f7+WLVumlStXamJiQmNjY5Kk8fFx9fT0qKSkZJaeBgDgeranaPLy8tTe3q6amhrF\n43E1NjYqEAioo6NDktTc3KxNmzbpwIED8vl8WrRokV577TVJ0ujoqOrr6yVJU1NT2rp1qzZs2DDL\nTwcA8CXedDsH5fobD+dy/vmbXcrt/ObvOxnfIm+6DQB3LgoeAAxFwQOAoSh4ADAUBQ8AhqLgAcBQ\nFDwAGIqCBwBDUfAAYCgKHgAMRcEDgKEoeAAwFAUPAIai4AHAUBQ8ABiKggcAQ1HwAGAoCh4ADEXB\nA4ChbN9021RLlhRobOxitmPMsHjxcl2+/Fm2YwAwxB35ptu5/sa95J8tufym1VJu5zd/38n4FnnT\nbQC4cyUt+O7ubvn9fhUXF6utre2GY1544QUVFxerrKxMx48fT2tu7otkO8BtimQ7wG2KZDvAbYpk\nO8BtiGQ7wG2KZDvA7LNsTE1NWQ8++KB19uxZa3Jy0iorK7NOnjw5bcyf//xna+PGjZZlWVZ/f79V\nWVmZ8twvTg/ZRZgVkizJytDHrzL4WKl9L8ifvfyZzZ7r+dl30s2fSals0/YIPhqNyufzyev1yul0\nqqGhQV1dXdPG7N+/X88995wkqbKyUpcuXdLo6GhKcwEAs8e24EdGRuTxeBLLbrdbIyMjKY356KOP\nks4FAMwe28skr/3FOrlrvy3culS3k1mZ3GZrxh4p9e8F+f/fXOfP9P6ay/nZd76UnR6zZ1vwLpdL\nsVgssRyLxeR2u23HDA8Py+126+rVq0nnSrf/nwMA4MZsT9EEg0ENDg5qaGhIk5OT6uzsVDgcnjYm\nHA5r9+7dkqT+/n4tW7ZMK1euTGkuAGD22B7B5+Xlqb29XTU1NYrH42psbFQgEFBHR4ckqbm5WZs2\nbdKBAwfk8/m0aNEivfbaa7ZzAQBzZNav5UnBH/7wB2v16tXWggULrKNHj2Y7TsoOHjxorVq1yvL5\nfNZLL72U7Thpef755617773X+uY3v5ntKGk7d+6cFQqFrNWrV1sPP/ywtXPnzmxHSsvnn39urV27\n1iorK7MCgYD185//PNuRbsnU1JS1Zs0aa/PmzdmOkrYHHnjAKikpsdasWWNVVFRkO05aLl68aD31\n1FOW3++3AoGA9de//vWmY+dFwf/jH/+wTp8+bYVCoZwp+FSv85+vDh8+bB07diwnC/78+fPW8ePH\nLcuyrLGxMeuhhx7Kqe+9ZVnW+Pi4ZVmWdfXqVauystI6cuRIlhOlb8eOHdb3vvc968knn8x2lLR5\nvV7r008/zXaMW/Lss89au3btsizr2v5z6dKlm46dFy9V4Pf79dBDD2U7Rlpy/Tr/qqoqLV++PNsx\nbsl9992nNWvWSJLy8/MVCAT00UcfZTlVeu6++25J0uTkpOLxuAoKCrKcKD3Dw8M6cOCAfvjDH+bs\nhRK5mPs///mPjhw5oh/84AeSrp0KX7p06U3Hz4uCz0Wp3COA2Tc0NKTjx4+rsrIy21HS8r///U9r\n1qzRypUr9e1vf1urV6/OdqS0/PSnP9VvfvMbLViQmxXicDi0fv16BYNBvfLKK9mOk7KzZ8+qsLBQ\nzz//vB555BE1NTVpYmLipuPn7F+nurpaJSUlMz7+9Kc/zVWEjJqP17zeaa5cuaKnn35aO3fuVH5+\nfrbjpGXBggX6+9//ruHhYR0+fFiRSCTbkVL29ttv695771V5eXlOHgVLUl9fn44fP66DBw/qd7/7\nnY4cOZLtSCmZmprSsWPH9OMf/1jHjh3TokWL9NJLL910/Jy9Hvw777wzV5uaE6ncI4DZc/XqVT31\n1FP6/ve/r9ra2mzHuWVLly7VE088ob/97W8KhULZjpOS9957T/v379eBAwf03//+V5cvX9azzz6b\nuFw6FxQVFUmSCgsLVVdXp2g0qqqqqiynSs7tdsvtdquiokKS9PTTT9sW/Lz7/SpXjgi4zj97LMtS\nY2OjVq9erZ/85CfZjpO2f//737p06ZIk6fPPP9c777yj8vLyLKdK3a9//WvFYjGdPXtWb7zxhh5/\n/PGcKveJiQmNjY1JksbHx9XT06OSkpIsp0rNfffdJ4/How8++ECS1Nvbq4cffvim4+dFwb/11lvy\neDzq7+/XE088oY0bN2Y7UlJfvc5/9erV+u53v5tT1/lv2bJF3/rWt/TBBx/I4/Ek7l/IBX19fdqz\nZ48OHTqk8vJylZeXq7u7O9uxUnb+/Hk9/vjjWrNmjSorK/Xkk0/qO9/5TrZj3bJcO1358ccfq6qq\nKvH937x5szZs2JDtWCl7+eWXtXXrVpWVlenEiRP65S9/edOxWX9HJwDA7JgXR/AAgMyj4AHAUBQ8\nABiKggcAQ1HwwHWGhoZy5rI5wA4FDwCGouABG//85z/1yCOP6OjRo9mOAqRtzl6qAMg1p0+f1pYt\nW/T73/+eUzbISRQ8cAOffPKJamtr9dZbb8nv92c7DnBLOEUD3MCyZcv0wAMP5MyrDAI3whE8cANf\n+9rX9Mc//lE1NTXKz8/Xli1bsh0JSBsFD9yAw+HQ3XffrbffflvV1dVavHixNm/enO1YQFp4sTEA\nMBTn4AHAUBQ8ABiKggcAQ1HwAGAoCh4ADEXBA4Ch/g8UXEPsb8R9QAAAAABJRU5ErkJggg==\n"
      }
     ],
     "prompt_number": 81
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 72
    }
   ],
   "metadata": {}
  }
 ]
}