distributions.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "# Distributions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Some notes on notation\n",
    "\n",
    "\"$|$\" and \";\" are used interchangably for conditional distribution. The main difference is that *variables* behind \";\" are considered fixed parameters (*hyper-parameters*) whereas those behind \"$|$\" are considered model parameters (true unknown). Of course, we can essentially model the hyper-parameters as model parameters and have \"$hyper^2$\"-parameter to model them. So the use of \"$|$\" and \";\" is really for the sake of increasing clarity (whenever possible) but the notation is by no mean to be completely rigorous. \n",
    "\n",
    "E.g., $Norm(x|\\mu;\\sigma^2)$ denotes a normal pdf with unknown mean and known variance\n",
    "\n",
    "$Norm(x|\\mu, \\sigma^2)$ indicates that both mean and variance are unknown"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Some notes on notation\n",
    "\n",
    "We denote $Norm(x; \\mu, \\sigma^2)$ or $Norm(x|\\mu, \\sigma^2)$ as the normal pdf. On the other hand, we use\n",
    "$$ X \\sim Norm(\\mu, \\sigma^2) $$ \n",
    "to indicate that $X$ is a normal random variable with mean $\\mu$ and variance $\\sigma^2$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Bernoulli Distribution (single trial)\n",
    "\n",
    "$$ Bern(x|p) = p^x (1-p)^{1-x},$$\n",
    "\n",
    "where $X$ is 0 or 1.\n",
    "\n",
    "$$ E[X]=p $$\n",
    "$$ Var[X] = p(1-p) $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Binomial Distribution ($N$ trials)\n",
    "\n",
    "$$ Bin(x|p)={N \\choose x} p^x (1-p)^{N-x}, $$\n",
    "\n",
    "where $X \\in \\{0,1,\\cdots,N\\}$.\n",
    "\n",
    "$$ E[X] = Np$$\n",
    "$$ Var[X] = Np(1-p)$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.text.Text at 0x7f7c85989cd0>"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYkAAAEMCAYAAAAxoErWAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAHfxJREFUeJzt3XucHGWd7/HPN0CQW7iImjUQAoscBa9BES/AKApBkHjc\ndQ1eYJWzykEOeNsFPeshqHsUV1ddXVfdk2UB0SAIAhEleGTwskgiBLklJCKExJCgGCIXkRC++0fV\nhJ5O10wPmZrunnzfr9e8pqvqqXp+XTNdv3qep6patomIiGhlQqcDiIiI7pUkERERlZIkIiKiUpJE\nRERUSpKIiIhKSRIREVEpSSIiIiolSURERKUkiS2MpFslHdqt9Uq6S9Jrx6KukW5rNGKr2vZok/R/\nJZ1ax7Z7kaTrJT2v03H0oiSJcUbS3ZIekfQHSfdLukLSlIHltp9v+8djHddo1dvw/tZJ+r2kn0p6\nrySNpK52D/ijub+a66zrbyFpd+CdwNdGebt7Sfpeud9XSfqSpNqOIZJ2lXSppIfKfXfcEGX7Jf2x\n/L9/UNLipiL/CHyirljHsySJ8cfA0bYnAX8G3Ad8qbMhjaqB97czsBfwaeB0YM5oViJpq9Hc3hj7\na+BK238a5e1+BVgDPAt4MXAYcPIo19Fc36PAM4B3AP86RGvAwMm2J9neyXZzuSuA10h6Vn3hjk9J\nEuOTAGw/BlwM7L9xQcPZbPn6Q5J+KWmtpG9JmthQ9rmSrimX3SLpjU3b+XC57oOS/k3SMyVdWZ7N\nzZe0c0W9p0v6VVnuVklveorv70Hb84C3AidI2r+irpVlXYslvUbSecBU4Ipy/ocb1vs7Sb8EHpK0\nVYsWx0GSbitbaXOa9tcTkvZpmD5H0sfL1wN1zivr/Nvmbbexvyv/Vk2OAq4d4T5tx97At22vt30f\n8APggBrqQdL2wJuBv7f9R9s/Ay6naCFVrla1oEyYNwBHjGqgW4AkiXGs/KC9FbhuiGJvofjg7A28\niOIsFElbU5x9/YDiTO5U4AJJz2lY983A4cB+wLHAlcAZwNOBrcp1WvkV8KqytXMW8I3NOcOzvRBY\nCRzSOF/SfsD7gAPLuo4E7rZ9PHAPcEx55vnZhtVmURxkd7G9oUV1bwNeD/w58N+Av28MZYgYB+o8\nuqzzH5tibWd/t/xbtfAC4I6qWJqVXZJry26k5t+XNxT9AnCcpO3KLsyjgO+3W88I69sPeNz2nQ2r\n/pKhk9KnJN0n6SeSDmuxfDHFfosR2LrTAUQtvivpcWAniu6BI4co+0Xba6D48FJ0IwC8AtjB9tnl\n9DWS5gHHAR8v533J9u/KdX8CrLF9czl9KdCyz9/2dxpeXyTpo8BBFAfJp2oVsFvTvA3AROD5ku63\nfU/T8lZnnl+0vWqIer40sFzSP1B05f2fIbbXrKrMwQy/v6v+Vs12AR7cWKE0vdz+s4GFFJ/7o22/\nG8D2G1ttpIUfA+8B/kBxgnmu7ctbFZR0APByilbsT4BnAo/ZPrfN+nYE1jXNW0fxP93K3wG3A49R\n7LMrJL3I9l0NZR4EJrdRdzRIS2J8mml7N4oD5P8CfizpmRVl1zS8foTiwwnFeMaKprLLgSkN043r\n/rHF9I60IOl4SYvKs8e1FGeHuw/xftoxBbi/cUZ5Fvp+YDawRtI3JQ13kFg5guXLKfbTaHg2I9vf\njX+rZmsZfDB9BrAE2N/2ZWWSbnWmXUmSgKsoui+3p/h77Sbp7IpV9qA4859m+zLgAuB/j6DKh4BJ\nTfMm0ZD8GtleaPvhsivsPOBnwBuaiu0EPDCCGIIkifFqoM/eti+lOKN+9Qi3sQrYs2neVOA3mxWY\nNBX4OsUg4662dwVuo72z8KptvoziIPvT5mW259o+hGKQG2DgoFbVNTTcF6w07pO9KPbTgEcoDqAD\nmhPSUNsezf19M0V3TVGpfRVFF9k3ACS9guIATjl9ZTmu9IcWP98ri+1GceD/l/JAvBY4h6LLaRNl\nnUcA88pZ04GBVmc79S0Ftpb05w2bfRHF/0o7zKb/U89rfN/RniSJcU7STIruh9tHuOr1wMPlQO7W\nkvqAY4BvbWZIOwBPAL+TNEHSu4DnP5UNSdpJ0kBM59u+vWn5fuVA9USKbog/Ao+Xi9cA+zBy75M0\nRdJuwEeAuQ3LFgFvK9/XDDY9W189RJ2jub+vBPqa5r0W+P/l6xOA8wYGxm2/obwiaFKLn6PLMvcD\ndwH/sxzQ36XcTmOyOUfSvzfUeQRPDqC/E/jsCOp7BLgE+Lik7SW9imLc6/zmNytpZ0lHSNq2jO3t\nFONTVzWUmQgcCFw9sl0ZSRLj08BVO+sorg0/3vaSclnj2exQA63rKT6Ub6A4A/wy8E7byyrWHe4M\n3OV2FwOfA35OcdA8gMEtgHa+KvGK8r3dQ3Gg/izw7hbb2JbiEtnfUpypP4Mnuzw+BXysHCz94BB1\nN++vbwLzKQbffwX8Q8Py91Pss7UU/eKXNm3r0w11fqhx209hfw/lPOAoSdsCSNoOWGt7oI//IYoT\nh9Uj2CYUFyocRbE/lwLrgQ80LN+TopsHSTtQXCp7iKS/ARbavmSE9b2PomV2H0V31Unl/89Aa+SM\nstw2wCfLcr8t15vZsO8AZgLX2B7pe97iqe6vLy3PqL5AkZDmNAzMDSx/L8UfdQNFf+N7Bg5okj5C\n8eF/HDjN9vxag40YJyR9ErjP9j+PUX3bADcBL7S9oWyl9Nn+0FjUPxxJ1wEnNrc2Y3i1JgkVd2Mu\npbhMchXFlRWzGs5qkbSj7YfK12+k6Ks+SsU17xcAL6PoC/0h8BznS7kjulp52e6/UQy8n2Y7g8U9\nrO5LYA8CltleDiBpLkWzb2OSGEgQpR0p+quhaHrPtf04cLekZeX2rq855ojYDGU3T1+n44jRUXeS\nmMLgy/pWUhzoB5F0MvBBir7FgWvrpzD4JrDfMPhywIiIqFndA9etLmvcpLvI9lds70vxDJ6PjWTd\niIioT90tiZUU13oP2IPB15U3uxD4asO6jdeNt1xXUhJHRMRTYHvY+5PqbkksBPZV8YjhiRTPxRl0\nG7+kfRsmj6EY6KYsN0vSREl7A/sCC1pVYrvrfs4888yOx5CYEtOWGFdiau+nXbW2JFxcCncKxXXl\nA5fALpZ0FsV10/OAUyS9juJmp7UUN+hg+3ZJ36a4CWw9xVVPaTVERIyh2h/wZ/sHFE/LbJx3ZsPr\n9w+x7qcobnqKiIgOyB3XNenr6+t0CJtITO1JTO3rxrgS0+iq/Y7ruklKL1RExAhJwl0wcB0RET0s\nSSIiIiolSURERKUkiYiIqJQkERERlZIkIiKiUpJERERUSpKIiIhKSRIREVEpSSIiIiolSURERKUk\niYiIqJQkERERlZIkIppMnjwNSWP+M3nytE6/9YhN5FHhEU0kAZ34n9KIvlYyYnPkUeEREbHZkiQi\nIqJSkkRERFRKkoiIiEpJEhERUSlJIiIiKiVJREREpSSJiIiolCQRERGVkiQiIqJSkkRERFSqPUlI\nmiFpiaSlkk5vsfwDkm6TdJOkqyXt2bBsg6QbJS2S9N26Y42IiMFqfcCfpAnAUuBwYBWwEJhle0lD\nmcOA620/KukkoM/2rHLZH2xPGqaOPOAvRlUe8Bdbgm55wN9BwDLby22vB+YCMxsL2L7W9qPl5M+B\nKQ2Lh30DERFRn7qTxBRgRcP0SgYngWYnAt9vmN5W0gJJ/ylpZtVKERFRj61r3n6rlkDL9rSkdwAH\nAoc1zJ5qe7WkvYEfSbrZ9l01xBkRES3UnSRWAlMbpvegGJsYRNLrgI8Ah5bdUgDYXl3+vktSP/AS\nYJMkMXv27I2v+/r66OvrG5XgIyLGi/7+fvr7+0e8Xt0D11sBd1AMXN8LLACOs724ocxLgIuAI23f\n2TB/F+AR249J2h34GTCzcdC7LJeB6xhVGbiOLUG7A9e1tiRsb5B0CjCfYvxjju3Fks4CFtqeB3wG\n2AG4SMWnc7ntNwHPA74maUO57qeaE0RERNQr33Ed0SQtidgSdMslsBER0cOSJCIiolKSREREVEqS\niIiISkkSERFRKUkiIiIqJUlERESlJImIiKiUJBEREZWSJCIiolKSREREVEqSiIiISkkSERFRKUki\nIiIqJUlERESlJImIiKiUJBEREZWSJCIiolKSREREVEqSiI6aPHkaksb8Z/LkaZ1+6yOS/RSdol7/\n4nVJ7vX3sCWTBHTi7yeq/m8S06CaK2OK3iYJ2xquXFoSERFRKUkiIiIqJUlERESlJImIiKiUJBER\nEZWSJCIiolKSREREVEqSiIiISrUnCUkzJC2RtFTS6S2Wf0DSbZJuknS1pD0blp1QrneHpOPrjjUi\nIgar9Y5rSROApcDhwCpgITDL9pKGMocB19t+VNJJQJ/tWZJ2BX4BTAcE3ABMt72uqY7ccd3DuvFO\n4sQ0qObccT1Odcsd1wcBy2wvt70emAvMbCxg+1rbj5aTPwemlK+PBObbXmf7AWA+MKPmeCMiokHd\nSWIKsKJheiVPJoFWTgS+X7Hub4ZZNyIiRtnWNW+/VVOmZdtV0juAA4HDRrru7NmzN77u6+ujr69v\nJDFGRIx7/f399Pf3j3i9usckDgZm255RTp8B2PbZTeVeB3wRONT2/eW8WRTjEyeV018FrrF9YdO6\nGZPoYd3Y156YBtWcMYlxqt0xibqTxFbAHRQD1/cCC4DjbC9uKPMS4CLgSNt3NsxvHLieUL4+sByf\naKwjSaKHdePBLzENqjlJYpxqN0nU2t1ke4OkUygGnScAc2wvlnQWsND2POAzwA7ARSo+Ccttv8n2\nWkmfoEgOBs5qThAREVGvfOlQdFQ3niEnpkE1pyUxTnXLJbAREdHDkiQiIqJSkkRERFRKkoiIiEpJ\nEhERUSlJIiIiKiVJREREpSSJiIio1FaSkPQdSUeX3w8RERFbiHYP+v8KvA1YJunTkp5bY0wREdEl\n2koStn9o++0UD9u7G7ha0n9KepekbeoMMCIiOqft7iNJTwf+GvgfwCKKR3tPB66uJbKIiOi4tp4C\nK+kS4LnA+cAbbd9bLrpQ0i/qCi4iIjqrrafASnqD7Sub5m1r+0+1RdamPAW2t3Xj000T06Ca8xTY\ncWq0nwL7yRbzrhtZSBER0WuG7G6SNBmYAmxXfoPcQNaZBGxfc2wREdFhw41JHEkxWL0H8E8N8x8E\nPlpTTBER0SXaHZP4C9vfGYN4RixjEr2tG/vaE9OgmjMmMU6NyndcS3qH7W8A0yR9sHm57X9qsVpE\nRIwTw3U37VD+3rHuQCIiovu01d3UzdLd1Nu6sRslMQ2qOd1N49RodTf981DLbZ860sAiIqJ3DNfd\ndMOYRBEREV0p3U3RUd3YjZKYBtWc7qZxarS6m75g+/2SrqDFf6jtYzcjxoiI6HLDdTedX/7+bN2B\nRERE92m7u0nSRIonwRq4w/ZjdQbWrnQ39bZu7EZJTINqTnfTODUq3U0NGzsa+CpwJ8Xzm/aW9F7b\n39+8MCMiopu1+xTYzwGvsd1n+zDgNcDn21lR0gxJSyQtlXR6i+WHSLpB0npJb25atkHSjZIWSfpu\nm7FGRMQoaaslATxo+1cN07+meMjfkCRNAL4MHA6sAhZKusz2koZiy4ETgA+32MTDtqe3GWNERIyy\n4a5uGjiz/4WkK4FvU3SMvgVY2Mb2DwKW2V5ebm8uMBPYmCRs31Mua9XxOWx/WURE1Ge4lsQbG16v\nAQ4rX/8W2K6N7U8BVjRMr6RIHO3aVtIC4HHgbNuXjWDdiIjYTEMmCdvv2sztt2oJjORSiam2V0va\nG/iRpJtt39VcaPbs2Rtf9/X10dfXN9I4IyLGtf7+fvr7+0e8XrvfJ/E04ETgAOBpA/Ntv3uY9Q4G\nZtueUU6fUazms1uUPQe4wvYlFdtquTyXwPa2bry0MzENqjmXwI5To/0d1+cDkym+qe5aim+qG3bg\nmmLcYl9Je5X3WcwCLh+i/MaAJe1SroOk3YFXAre3GW9ERIyCdlsSi2y/pOzueaGkbYCf2D64jXVn\nAF+kSEhzbH9a0lnAQtvzJL0UuBTYBXgUWG37BZJeAXwN2FCu+3nb/9Fi+2lJ9LBuPENOTINqTkti\nnGq3JdFuklhg+yBJPwZOBlYDC2zvs/mhbp4kid7WjQe/xDSo5iSJcWpU77gGvi5pV+BjFN1FO5av\nIyJiHMujwqOjuvEMOTENqjktiXFqVAeuJT1d0pfKR2TcIOkLkp6++WFGREQ3a/fqprnAfcBfAH8J\n/A64sK6gIiKiO7Q7cH2r7ec3zbvF9gtqi6xN6W7qbd3YjZKYBtWc7qZxarTvk5gvaZakCeXPXwFX\nbV6IERHR7YZsSUh6kOL0RcAOwBPlognAQ7Yn1R7hMNKS6G3deIacmAbVnJbEODUql8Da3mn0QoqI\niF7T7n0SSDoWOLSc7Lc9r56QIiKiW7R7CeyngdMonp10O3BaOS8iIsaxdq9uuhl4se0nyumtgEW2\nX1hzfMPKmERv68a+9sQ0qOaMSYxTo311ExQP4Buw88hDioiIXtPumMSngEWSrqG40ulQ4CO1RRUR\nEV1h2O4mFe3cPSi+QvRlFEnietur6w9veOlu6m3d2I2SmAbVnO6mcWq0HxXeFXdXt5Ik0du68eCX\nmAbVnCQxTo32mMSNkl62mTFFRESPabclsQR4DnA38DBFl5NzdVNsrm48Q05Mg2pOS2KcGu0vHTpy\nM+OJiIgeNGSSkPQ04CRgX+AWiu+ofnwsAouIiM4bbkziXOClFAniKOBztUcUERFdY7jupv0HrmqS\nNAdYUH9IERHRLYZrSawfeJFupoiILc9w3yexgeJqJiiuaNoOeIQnr27K90nEZunGq3YS06Cac3XT\nODVa3yex1eiFFBERvWYkD/iLiIgtTJLEFmTy5GlIGvOfyZOndfqtR0068T+V/6ex1dYd190sYxLt\n68Z+7cQ0qOaeigk6FVfGSUZDHd8nERERW5jak4SkGZKWSFoq6fQWyw+RdIOk9ZLe3LTshHK9OyQd\nX3esERExWK3dTZImAEuBw4FVwEJglu0lDWWmApOADwOX276knL8r8AtgOsUltzcA022va6oj3U1t\n6sYui8Q0qOaeignS3dTLuqW76SBgme3lttcDc4GZjQVs32P7Vjb9TzsSmG97ne0HgPnAjJrjjYiI\nBnUniSnAiobpleW8p7Lub0awbkREjIJ2HxX+VLVqyrTbTmx73dmzZ2983dfXR19fX5tVRERsGfr7\n++nv7x/xenWPSRwMzLY9o5w+g+JxHme3KHsOcEXDmMQsoM/2SeX0V4FrbF/YtF7GJNrUjf3aiWlQ\nzT0VE2RMopd1y5jEQmBfSXtJmgjMAi4fonxjwFcBr5e0czmI/fpyXkREjJFak4TtDcApFIPOtwFz\nbS+WdJakYwAkvVTSCuAvga9KuqVcdy3wCYornK4HzioHsCMiYozkjustSDd2WSSmQTX3VEyQ7qZe\n1i3dTRER0cOSJCIiolKSREREVEqSiIiISkkSERFRKUkiIiIqJUlERESlJImIiKiUJBEREZWSJCIi\nolKSREREVEqSiIiISkkSERFRKUkiIiIqJUlERESlJImIiKiUJBEREZWSJCIiolKSREREVEqSiIiI\nSkkSERFRKUkiIiIqJUlERESlJImIiKiUJBEREZWSJCIiolKSREREVEqSiIiISrUnCUkzJC2RtFTS\n6S2WT5Q0V9IySddJmlrO30vSI5JuLH++UnesEREx2NZ1blzSBODLwOHAKmChpMtsL2kodiLwe9vP\nkfRW4DPArHLZr2xPrzPGiIioVndL4iBgme3lttcDc4GZTWVmAueWry+mSCgDVHN8ERExhLqTxBRg\nRcP0ynJeyzK2NwAPSNqtXDZN0g2SrpH06ppjjYiIJrV2N9G6JeBhyqgscy8w1fZaSdOB70ra3/ZD\nNcQZEREt1J0kVgJTG6b3oBibaLQC2BNYJWkrYJLtteWyxwBs3yjpTmA/4MbmSmbPnr3xdV9fH319\nfaMUfkTE+NDf309/f/+I15PdfGI/esqD/h0U4wz3AguA42wvbihzMvB82ydLmgW8yfYsSbtTDGg/\nIWkf4FrgBbYfaKrDdb6H8UQaaKSNec1U/Y0S06Caeyom6FRcQ8cU7ZGE7WHHfWttSdjeIOkUYD7F\n+Mcc24slnQUstD0PmAOcL2kZcD9PXtl0KPBxSeuBDcB7mxNERETUq9aWxFhIS6J93Xg2mpgG1dxT\nMUFaEr2s3ZZE7riOiIhKSRIREVEpSSIiIiolSURERKUkiYiIqJQkERERlZIkIiKiUpJERERUSpKI\niIhKSRIREVEpSaImkydPQ9KY/kyePK3Tbzui4zrx2RvPn788u6km3fhMm258/k9iGlRzT8UE+T9v\nqrmnnimVZzdFRMRmS5KIiIhKSRIREVEpSSIiIiolSURERKUkiYiIqJQkERERlZIkIiKiUpJERERU\nSpKIiIhKSRIREVEpSSIiIiolSURERKUkiYiIqJQkERERlZIkIiKiUu1JQtIMSUskLZV0eovlEyXN\nlbRM0nWSpjYs+0g5f7GkI+qONSIiBqs1SUiaAHwZOBI4ADhO0nObip0I/N72c4AvAJ8p190f+Cvg\necBRwFdUfOVURERP6e/v73QIT1ndLYmDgGW2l9teD8wFZjaVmQmcW76+GHht+fpYYK7tx23fDSwr\nt7eJfJ9tRHSzY455U88ep7be/Lc/pCnAiobplWx6oN9YxvYGSesk7VbOv66h3G/KeS2M/ffKrlmT\nRk1EtOfhh9fRq8epulsSrSJs3lNVZdpZNyIialR3S2IlMLVheg9gVVOZFcCewCpJWwE7214raWU5\nf6h1S505qx9+iGTs4+rGmGC4uBLTxlp7LibI/3lDrT359xta3UliIbCvpL2Ae4FZwHFNZa4ATgCu\nB94C/KicfzlwgaTPU3Qz7QssaK7Advp9IiJqUmuSKMcYTgHmU3RtzbG9WNJZwELb84A5wPmSlgH3\nUyQSbN8u6dvA7cB64GTb6W6KiBhDynE3IiKq9PQd18PdqNcJkuZIWiPp5k7HAiBpD0k/knS7pFsk\nndrpmAAkbSvpekmLyrjO7HRMAyRNkHSjpMs7HQuApLsl/bLcV5t0uXaCpJ0lXVTe6HqbpJd3OJ79\nyv1zY/l7XRf9r39A0q2SbpZ0gaSJXRDTaeXnbthjQs+2JMob9ZYCh1MMaC8EZtle0uG4Xg08BJxn\n+4WdjKWMZzIw2fZNknYEbgBmdno/AUja3vYj5QULPwNOtd3xg6CkDwAHApNsH9sF8fwaOND22k7H\nMkDSfwDX2j5H0tbA9rb/0OGwgI3HhpXAy22vGK58zbE8G/gp8Fzbj0m6EPie7fM6GNMBwLeAlwGP\nAz8ATrJ9Z6vyvdySaOdGvTFn+6dA13yYba+2fVP5+iFgMZX3m4wt24+UL7elGB/r+BmLpD2ANwD/\nr9OxNBBd9FmVtBNwiO1zAMobXrsiQZReB9zZ6QTRYCtgh4FkSuVVmmPmecDPbf/J9gbgWuC/VxXu\nmn+8p6DVjXpdcfDrVpKmAS+muJKs48punUXAauBq2ws7HRPweeBv6YKE1cDAVZIWSvqbTgcD7AP8\nTtI5ZffO1yVt1+mgGryV4ky542yvAj4H3ENxQ/ADtn/Y2ai4FThU0q6Stqc4KdqzqnAvJ4ncbDcC\nZVfTxcBpZYui42w/YfslFPfAvFzF87o6RtLRwJqy5SU6dWH7pl5p+6UUH+b3lV2anbQ1MB34F9vT\ngUeAMzobUkHSNhSP9Lmo07EASNqFoodjL+DZwI6S3tbJmMqu5rOBHwJXAjdRdDu11MtJop0b9QIo\nm7kXA+fbvqzT8TQruyr6gRkdDuVVwLHlGMC3gNdI6ljf8QDbq8vfvwUupeIZZmNoJbDC9i/K6Ysp\nkkY3OAq4odxX3eB1wK9t/77s2rkEeGWHY8L2ObYPtN1H0T2+rKpsLyeJjTfqlVcLzKK4Aa8bdNNZ\nKMC/A7fb/mKnAxkgaXdJO5evt6P4MHV0MN32R21Ptb0Pxf/Tj2wf38mYJG1ftgKRtANwBEV3QcfY\nXgOskLRfOetwivuZusFxdElXU+ke4GBJT1Nx6/PhFOOCHSXpGeXvqRTjEZX7rO47rmtTdaNeh8NC\n0jeBPuDpku4BzhwY4OtQPK8C3g7cUvb/G/io7R90KqbSnwHnlleiTAAutH1lh2PqRs8CLpVkis/r\nBbbndzgmgFMpnoiwDfBr4F0djqfxZOM9nY5lgO0Fki4GFlHcFLwI+HpnowLgO+WDVAduVF5XVbBn\nL4GNiIj69XJ3U0RE1CxJIiIiKiVJREREpSSJiIiolCQRERGVkiQiIqJSkkRERFRKkoiIiEr/BexL\njQk9NpwRAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f7c85a12410>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import scipy.stats as stats\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "\n",
    "N=8\n",
    "p=0.5\n",
    "ps=[stats.binom.pmf(i,N,p) for i in xrange(N+1)]\n",
    "fig, ax = plt.subplots()\n",
    "rects1 = ax.bar(xrange(N+1),ps)\n",
    "plt.title('Binomial Distribution ($n$='+str(N)+',$p$='+str(p)+')')\n",
    "plt.ylabel('Probability')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Poisson Distribution\n",
    "\n",
    "$$ Poisson(x|\\lambda T) = \\frac{e^{-\\lambda T} {\\lambda T}^x}{x!}, $$\n",
    "\n",
    "where $X$ is a non-negative integer, $\\lambda$ is rate of arrival and $T$ is the length of the observed period. \n",
    "\n",
    "$$ E[X] = \\lambda T$$\n",
    "$$ Var[X] = \\lambda T$$\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.text.Text at 0x7f7cb85a32d0>"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYwAAAEMCAYAAADXiYGSAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAHBlJREFUeJzt3XuYXWVh7/HvL8TghZv3HINJsNEKXqpBOWgLDqASVIiP\nrW28YqVPsYi32gpofUi0rUBthVOPVc+JHOSo8UYVKGjggfFyVBgBBSEhUSQkBqJyk0upIfzOH2sN\nLjZ7Zt6Z7DV7Lr/P8+xn1vVd75pJ9m+v933X2rJNRETEWOb0uwIRETE9JDAiIqJIAiMiIookMCIi\nokgCIyIiiiQwIiKiSAIjIiKKJDAiIqJIAiN2iqSfSDq43/Uo1cv6NsuS9HNJh/ai3M6y2yDpHyW9\ns63y2yLpMkn79rses1UCIwCQdKOkeyX9RtLNkj4j6dFj7Wf72ba/PRl1HEvjHO6UdJuk70o6VpKG\ntymtb0kA9PLcO4/X5u9V0hOANwGfaqP8Lsd7u6QhSfdJ+kzB9oOS/rP+t3iXpHWN1f8EfLi92sZo\nEhgxzMArbe8BLAVeCPxdf6s0bsPnsCewCDgFOAFY3cuDSNqll+X1wVuAC2z/1/ACSW+S9JLmRpL2\nlHSlpHdLOqZ+814p6S8lfV3SGwqP9wuqN/nSv4OB42zvYXt3280rivOAQyQ9ubCs6KEERjQJwPbN\nwIXAswEk7SvpUkm3S7pG0pEP7tD4ZCzpBElb6k+G6yQd0tiu67qCst8r6cf1+i9Imld4DnfZPh/4\nM+BoSfuV1lfSZ4GFwPn1ur+t93ufpB8Dd0vapctVyAGSrpV0q6TVzbpKekDS0xrzZ0r60BjHO7Sx\n/TN7+Hs6AvhWx7JzgbM79jsCONz26cA3gd/aXmn708A/ANePcowH2f6a7XOB20q2r6nbwjrkrgBe\nPo6yokcSGPEwkp4KvAK4UtJcqjeTbwBPBN4JfE7S0zv2eQbwdmD/+irlcODG0dYVlv1aqjeHfYA/\noPp0XMz2ELAFOKi0vrbfDNxEfcVl+5/q3VZQvYnuZXtHl8O9HngZ8HvA7/PQK7QRn/I5yvGG6zqX\n6pN1r35Pz6Hjzd72ncDFwGsai2+y/at6+hAeGjL/CZxWB9RtXX6eO8rxS3xE0i8lfafzygdYR3WO\nMckSGNH0NUm3Ad8GLgU+AhwIPMb2qbbvt30pcD7wuo59dwDzgGdLmmv7Jts/H2NdSdln2N5m+w6q\nN83nTeC8tgKPG0d9h3V+yj3D9tZmU06Hf63X30H1Cfz1o5TVzUjb9Pr3tBdwV5flq4G/GJ6x/b3G\nukOASxrrrrF9qO3H2n5cl59HjXqmo3sf8DRgAfC/gPMk7dNYf1d9DjHJEhjRtLz+z76P7XfUb4xP\nATZ3bLeJ6j/zg2z/DHg3sBLYJunzkv5bl3W/bKwrKXtbY/peYLcJnNcC4NbS+o5iyzjWbwLGKq9U\nr39PtwO7d1n+M2CppMVd1h1C9SGidbaHbN9je7vtzwL/j+qKd9juwB2TUZd4qARGNHX7hLsVeGrH\nsoVUHZkPYXuN7YOoOpyh6nTuXLewsW5rY37UsidK0gup3nC/O5760r0Jaawvj2n+nhZRnd+we4Hm\nqLP54yi7+G9Q6GrgGc0FkpZQNdGdAhzTsW4x8Cjb13Ysv6DuCP9Nl9d/TLBu3ZiH/tvcF/hxD8uP\nQgmMGMtlwD11h+9cSQPAq4AvNDeS9AxJh9Sdpr+lauPeMcq6++uy7x6r7ImQtLuk4bLOtn1daX1r\nt1A1i4zH2yUtkPQ44CRgTWPdVcDrJc2RtAzobJcf7XhFf4NxuAAYGJ6RdADwfqqrrbOAN0pqvkEf\nSperC9uvqEcx7dHl9cpG+btIeiSwCzBX0q4aYaSZqpFZLx/eRtVIrIOoOt2p/177AxdN8NxjJyQw\nYljXT7i2twNHUTUJ/Br4OPAm2xs79tuV6tPpr6g+ET+R6k1opHUfGEfZ43GepDupOpFPAj4KvLXL\neY5WX+p1H6w7cN87Ql3cMf15YC3w0/r1D43176Y619up+h7+vaOsEY/Xwu/ps8AR9Zvy84CTgb+y\nvcP2NuCHwJGS9pP0NuCvgL1UDa0t6Yvp9HdUV1gnAG+opz8wvLK+Ujmxnn0E8PfAL6n+Nm+naiod\nPtflwKW2b5lAPWInqe2vaK0/TZ1OFU6rbZ/asf5Yqn8UO6g6s/7S9npJi6hGQ6yvN/2B7eNarWzE\nLCFp+E3534Adth9orNuVqpN9PMNgJ4Wk7wPHdF4xxuRoNTAkzQE2AIdRfYobAlbYXt/YZjfbd9fT\nR1LdsHNEHRjn2X5uaxWMiIhibTdJHQBstL2pvqxeQ3VJ+aDhsKjtBjzQmJ/I5W9ERLSg7cBYwEOH\nA26hYzgmgKTjJP2Uqh23+UC0xZKuqO9w/aN2qxoREaNpOzC6XSE8rA3M9idsL6HqFPtgvfhmYKHt\n/YH3Ap+XNJEx+BER0QNzWy5/Cw8dZ783Dx2b3umLwCcBbP+Wargjtq+U9DOqseNXNneQ1G6vfUTE\nDGV7XM3+bV9hDAFLJC2qx0+voHp20IPqG4aGvYqqkxxJT6g7zVH10LYlwA3dDmJ7xr5OPvnkvtch\n55fzm43nN5PPzZ7Y5+xWrzBs75B0PNXY9OFhteskrQKGXD1N9HhJL6W6mrgdOLre/WDgQ5K2Uw25\nPdbVc3IiIqIP2m6SwvY3qJ7c2Vx2cmP63SPsdw5wTru1i4iIUrnTe4obGBjodxValfOb3mby+c3k\nc5uo1u/0bpskT/dziIiYbJLwFOv0joiIGSKBERERRRIYERFRJIExTvPnL0ZSK6/58xf3+/QiIkaU\nTu/xH4+JfU1DUekTvqEmImI80ukdERGtSWBERESRBEZERBRJYERERJEERkREFElgREREkQRGREQU\nSWBERESRBEZERBRJYERERJEERkREFElgREREkQRGREQUSWBERESRBEZERBRJYERERJHWA0PSMknr\nJW2QdEKX9cdKulrSVZK+LemZjXUnSdooaZ2kl7dd14iIGFmr37gnaQ6wATgM2AoMAStsr29ss5vt\nu+vpI4HjbB8haT/gc8ALgb2Bi4Gnd369Xr5xLyJi/KbiN+4dAGy0vcn2dmANsLy5wXBY1HYDHqin\njwLW2L7f9o3Axrq8iIjog7ktl78A2NyY30KXN31JxwF/DTwCOLSx7/cbm/2iXhYREX3QdmB0u9x5\nWJuL7U8An5C0Avgg8JbSfQFWrlz54PTAwAADAwPjr2lExAw2ODjI4ODgTpXRdh/GgcBK28vq+RMB\n2z51hO0F3G57r85tJX0DONn2ZR37pA8jImKcpmIfxhCwRNIiSfOAFcC5zQ0kLWnMvoqqk5x6uxWS\n5knaB1gCXN5yfSMiYgStNknZ3iHpeGAtVTittr1O0ipgyPb5wPGSXgr8FrgdOLre9zpJXwKuA7ZT\njZ7Kx++IiD5ptUlqMqRJKiJi/KZik1RERMwQCYyIiCiSwIiIiCIJjIiIKJLAiIiIIgmMiIgoksCI\niIgiCYyIiCiSwIiIiCIJjIiIKJLAiIiIIgmMiIgoksCIiIgiCYyIiCiSwIiIiCIJjIiIKJLAiIiI\nIgmMiIgoksCIiIgiCYyIiCiSwIiIiCIJjIiIKJLAiIiIIq0HhqRlktZL2iDphC7r3yPpWkk/knSR\npKc21u2QdKWkqyR9re26RkTEyGS7vcKlOcAG4DBgKzAErLC9vrHNS4DLbN8n6W3AgO0V9brf2N5j\njGO4zXPocjygreOJyTyXiJi9JGFb49mn7SuMA4CNtjfZ3g6sAZY3N7D9Ldv31bM/ABY0Vo/rZCIi\noj1tB8YCYHNjfgsPDYROxwAXNuZ3lXS5pO9JWj7SThER0b65LZff7Qqha5uLpDcC+wMvaSxeaPsW\nSfsAl0i62vbPO/dduXLlg9MDAwMMDAzsTJ2nlPnzF7Nt26ZWyn7ykxdxyy03tlJ2REwtg4ODDA4O\n7lQZbfdhHAistL2snj8RsO1TO7Z7KXAGcLDtW0co60zgPNvndCyf0X0Y6TOJiDZMxT6MIWCJpEWS\n5gErgHObG0h6PvBJ4KhmWEjaq94HSU8AXgxc13J9IyJiBK02SdneIel4YC1VOK22vU7SKmDI9vnA\nacBjgC+r+ji9yfargX2BT0naUe/7keboqoiImFytNklNhjRJ9fZ4ETE7TMUmqYiImCESGBERUSSB\nERERRRIYERFRJIERERFFEhgREVEkgREREUUSGBERUSSBERERRRIYERFRJIERERFFEhgREVEkgRER\nEUVmRGBIauU1f/7ifp9aRMSUMSMebz6THzeex5tHRBvyePOIiGhNAiMiIookMCIiokgCIyIiihQF\nhqSvSnqlpARMRMQsVRoA/wa8Htgo6RRJz2yxThERMQUVBYbti22/AVgK3AhcJOl7kv5c0iParGBE\nREwNxU1Mkh4PvAX4C+Aq4AyqALmolZpFRMSUUtqHcQ7wHeDRwJG2j7L9RdvvAHYbY99lktZL2iDp\nhC7r3yPpWkk/knSRpKc21h1d73e9pDeP79QiIqKXiu70lvQK2xd0LNvV9n+Nsd8cYANwGLAVGAJW\n2F7f2OYlwGW275P0NmDA9gpJjwV+SHUVI+AKYKntOzuOkTu9e3i8iJgd2rzT+++7LPt+wX4HABtt\nb7K9HVgDLG9uYPtbtu+rZ38ALKinDwfW2r7T9h3AWmBZYX0jIqLH5o62UtJ8qjfwR0l6PtUnfYA9\nqJqnxrIA2NyY30IVIiM5BrhwhH1/we/CJCIiJtmogUH1Kf8twN7AvzSW3wW8v6D8bpc7XdtAJL0R\n2B94yXj3hZWN6YH6FRERwwYHBxkcHNypMkr7MP7Y9lfHXbh0ILDS9rJ6/kTAtk/t2O6lVKOuDrZ9\na71sBVV/xtvq+U8Cl9r+Yse+6cPo4fEiYnaYSB/GqIEh6Y22/6+k99LlXcv2v3TZrbn/LsD1VJ3e\nNwOXA6+zva6xzfOBLwOH2/5ZY3mz03tOPb1/3Z/RPEYCo4fHi4jZYSKBMVaT1GPqn6MOnR2J7R2S\njqfqsJ4DrLa9TtIqYMj2+cBp9XG+rOrdcZPtV9u+XdKHqYLCwKrOsIiIiMmTL1AavfS+f+LPFUZE\ntKHnVxiS/sdo622/czwHi4iI6WusJqkrJqUWEREx5aVJavTS+95ElCapiGhDG01Sp9t+t6Tz6D5K\n6qhx1jEiIqapsZqkzq5/frTtikRExNRW3CQlaR7wTKorjett/7bNipVKk1RvjxcRs0Mb92EMF/xK\n4JPAz6ge2bGPpGNtXzj6nhERMVOUPhpkPfAq2z+t538P+A/bff+q1lxh9PZ4ETE7tPl487uGw6J2\nA9UDCCMiYpYYa5TUa+rJH0q6APgS1cfd11J9GVJERMwSY/VhHNmY3sbvHj3+K+BRrdQoIiKmpNy4\nN3rpfe9TSB9GRLShzVFSj6T6NrxnAY8cXm77reOqYURETFulnd5nA/OpvoHvW1TfwJdO74iIWaR0\nWO1Vtp8v6Wrbz5X0COA7tg9sv4pj1i1NUj08XkTMDm0Oq91e/7xD0rOBPYEnjedAERExvRX1YQCf\nrr8y9YPAuVTfwPfB1moVERFTTkZJjV5635uI0iQVEW1orUlK0uMl/aukKyVdIel0SY+fWDUjImI6\nKu3DWAP8Evhj4E+AXwNfbKtSEREx9ZSOkvqJ7Wd3LLvG9nNaq1mhNEn19ngRMTu0OUpqraQVkubU\nrz8Fvjn+KkZExHQ16hWGpLuoPt4KeAzwQL1qDnC37T1ar+EYcoXR2+NFxOzQ8ysM27vb3qP+Ocf2\n3Po1pzQsJC2TtF7SBkkndFl/UN2Rvr3xdNzhdTvqjvarJH1tPCcWERG9VXofBpKOAg6uZwdtn1+w\nzxzg48BhwFZgSNLXba9vbLYJOBr4my5F3GN7aWkdIyKiPaXDak8B3gVcV7/eVS8bywHARtubbG+n\nGm21vLmB7Zts/4Tu7S7julyKiIj2lF5hvAJ4nu0HACSdBVwFnDjGfguAzY35LVQhUmpXSZcD9wOn\n2v76OPaNiIgeKm6SAvYCbqun9yzcp9sVwnh6WRfavkXSPsAl9cMPf/7wzVY2pgfqV0REDBscHGRw\ncHCnyii9D+N1wCnApVQhcDBwku01Y+x3ILDS9rJ6/kTAtk/tsu2ZwHm2zxmhrK7rM0qqt8eLiNmh\nlfswVL1jfRc4EDinfr1orLCoDQFLJC2SNA9YQfXwwhEP1zjuXvU+SHoC8GKq/pOIiOiD0iuMCd/V\nLWkZcAZVOK22fYqkVcCQ7fMlvQD4d6omr/uAW2w/R9KLgE8BO+p9P2b7/3QpP1cYPTxeRMwOE7nC\nKA2Ms4CP2x6aaOXaksDo7fEiYnZoMzDWA08HbgTuoWo6su3nTqCePZXA6O3xImJ2mEhglI6SOnwC\n9YmIiBlk1MCQ9EjgbcAS4BqqPoj7J6NiERExtYw1Suos4AVUYXEE8M+t1ygiIqaksZ5W++DoKElz\ngcun2rOd0ofR2+NFxOzQxn0Y24cn0hQVETG7jXWFsYNqVBRUI6MeBdzL70ZJ5fswZtjxImJ26Pko\nKdu77FyVIiJipij9itaIiJjlEhgREVEkgREREUUSGPEQ8+cvRlLPX/PnL+73qUXETip6ltRUllFS\n0+V4GZEVMZW08n0YERERkMCIiIhCCYyIiCiSwIiIiCIJjIiIKJLAiIiIIgmMiIgoksCIiIgiCYyI\niCiSwIiIiCKtB4akZZLWS9og6YQu6w+SdIWk7ZJe07Hu6Hq/6yW9ue26RkTEyFp9lpSkOcAG4DBg\nKzAErLC9vrHNQmAP4G+Ac22fUy9/LPBDYCnVN/xdASy1fWfHMfIsqWlxvDxLKmIqmYrPkjoA2Gh7\nk+3twBpgeXMD2zfZ/gkPf5c6HFhr+07bdwBrgWUt1zciIkbQdmAsADY35rfUyyay7y/GsW9ERPTY\nqN/p3QPdLndK2yXGse/KxvRA/YqIiGGDg4MMDg7uVBltB8YWYGFjfm+qvozSfQc69r20+6Yrx12x\niIjZZGBggIGBgQfnV61aNe4y2m6SGgKWSFokaR6wAjh3lO2bVxXfBF4mac+6A/xl9bKIiOiDVgPD\n9g7geKoO62uBNbbXSVol6VUAkl4gaTPwJ8AnJV1T73s78GGqkVKXAavqzu+IiOiDfEXr6KXP4GGu\nk328DKuNmEqm4rDaiIiYIRIYERFRJIERERFFEhgREVEkgREREUUSGBERUSSBERERRRIYERFRJIER\nERFFEhgREVEkgREREUUSGBERUSSBERERRRIYERFRJIERERFFEhgREVEkgREREUUSGBERUSSBERER\nRRIYERFRJIERERFFEhgREVEkgREREUVaDwxJyyStl7RB0gld1s+TtEbSRknfl7SwXr5I0r2Srqxf\nn2i7rhERMbK5bRYuaQ7wceAwYCswJOnrttc3NjsGuM320yX9GXAasKJe91PbS9usY0RElGn7CuMA\nYKPtTba3A2uA5R3bLAfOqqe/QhUuw9Ry/SIiolDbgbEA2NyY31Iv67qN7R3AHZIeV69bLOkKSZdK\n+qOW6xoREaNotUmK7lcIHmMb1dvcDCy0fbukpcDXJO1n++6HF7myMT1QvyIiYtjg4CCDg4M7VYbs\nzvfv3pF0ILDS9rJ6/kTAtk9tbHNhvc1lknYBbrb9pC5lXQq81/aVHcv98Azq2RnQ+fuRhvMsx9vZ\nY0VE/0jC9ria/dtukhoCltQjnuZRdWaf27HNecDR9fRrgUsAJD2h7jRH0tOAJcANLdc3Jtn8+YuR\n1PPX/PmL+31qETNOq01StndIOh5YSxVOq22vk7QKGLJ9PrAaOFvSRuBWfjdC6mDgQ5K2AzuAY23f\n0WZ9Y/Jt27aJNq5otm3LeImIXmu1SWoypElquhyve5NUmsAi+mMqNklFRMQMkcCIiIgiCYyIiCiS\nwIiIiCIJjIiIKJLAiIiIIgmMiIgoksCIiIgiCYyIiCiSwIiIiCIJjIiIKJLAiIiIIgmMiIgoksCI\niIgiCYyIiCiSwIiIiCIJjIiIKJLAiIiIIgmMmFXmz1+MpJ6/5s9f3O9Ti2hdvtN79NJn8HdsT/bx\npsZ3euc7xCMq+U7viIhoTQIjIiKKtB4YkpZJWi9pg6QTuqyfJ2mNpI2Svi9pYWPdSfXydZJe3nZd\nIyJiZK0GhqQ5wMeBw4FnAa+T9MyOzY4BbrP9dOB04LR63/2APwX2BY4APqGqATpi2hgcHOx3FVo1\nk89vJp/bRLV9hXEAsNH2JtvbgTXA8o5tlgNn1dNfAQ6tp48C1ti+3/aNwMa6vIhpY6a/6czk85vJ\n5zZRbQfGAmBzY35LvazrNrZ3AHdKelyXfX/RZd+IKe2jHz09w3hjxmg7MLo1IXWOPRxpm5J9I6a0\ne+65k+qfbW9f27Zt6nq83GcSbZrbcvlbgIWN+b2BrR3bbAaeCmyVtAuwp+3bJW2pl4+2b629ro3u\n3SY5Xu+OleP1/ni9t23bphGPt2rVqkmrx2Sbyec2EW0HxhCwRNIi4GZgBfC6jm3OA44GLgNeC1xS\nLz8X+Jykj1E1RS0BLu88wHhvPImIiIlpNTBs75B0PLCWqvlrte11klYBQ7bPB1YDZ0vaCNxKFSrY\nvk7Sl4DrgO3Acc6ttBERfTPtHw0SERGTY1rf6T3WTYHTmaS9JV0i6TpJ10h6Z7/r1GuS5ki6UtK5\n/a5Lr0naU9KX65tOr5X03/tdp16S9B5JP5F0taTPSZrX7zrtDEmrJW2TdHVj2WMlrZV0vaRvStqz\nn3XcGSOc32n1v88fSfqqpD3GKmfaBkbhTYHT2f3AX9veD3gR8PYZdn4A76JqcpyJzgAusL0v8AfA\nuj7Xp2ckPQV4B7DU9nOpmrZX9LdWO+1MqveSphOBi23/PlXf6kmTXqve6XZ+a4Fn2X4e1X1uY57f\ntA0Mym4KnLZs32L7R/X03VRvODPmPhRJewOvAP53v+vSa5J2Bw6yfSZAffPpb/pcrV7bBXiMpLnA\noxlxBOP0YPu7wO0di5s3FZ8FvHpSK9VD3c7P9sW2H6hnf0A1EnVU0zkwSm4KnBEkLQaeRzWSbKb4\nGPC3zMx7a54G/FrSmXWT26clParfleoV21uBfwZuorqh9g7bF/e3Vq14ku1tUH2AA57Y5/q06a3A\nhWNtNJ0DY1bc2CdpN6pHpryrvtKY9iS9EthWX0GJNm806Y+5wFLgf9peCtxL1bwxI0jai+rT9yLg\nKcBukl7f31rFREn6ALDd9ufH2nY6B0bJTYHTWn25/xXgbNtf73d9eugPgaMk3QB8AThE0mf7XKde\n2gJstv3Dev4rVAEyU7wUuMH2bfXjfM4BXtznOrVhm6QnA0iaD/yyz/XpOUlHUzUNFwX+dA6MB28K\nrEdorKC62W8m+Qxwne0z+l2RXrL9ftsLbT+N6u92ie0397tevVI3Y2yW9Ix60WHMrM79m4ADJT2y\nfoL0YcyMTv3Oq91zgbfU00cD0/1D20POT9Iy4H3AUbb/q6SAtu/0bs1INwX2uVo9I+kPgTcA10i6\niqq57f22v9HfmkWhd1I9qeARwA3An/e5Pj1j+3JJXwGuorqp9irg0/2t1c6R9HlgAHi8pJuAk4FT\ngC9LeitVSL62fzXcOSOc3/uBecBF9WNffmD7uFHLyY17ERFRYjo3SUVExCRKYERERJEERkREFElg\nREREkQRGREQUSWBERESRBEZERBRJYERERJH/D/bgKs4BxEh+AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f7cb85a0d10>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import scipy.stats as stats\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "\n",
    "N=10\n",
    "lt=1.5\n",
    "ps=[stats.poisson.pmf(i,lt) for i in xrange(N+1)]\n",
    "fig, ax = plt.subplots()\n",
    "rects1 = ax.bar(xrange(N+1),ps)\n",
    "plt.title('Poisson Distribution ($\\lambda T$='+str(lt)+')')\n",
    "plt.ylabel('Probability')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Beta Distribution\n",
    "\n",
    "$$ Beta(x|a,b) = \\frac{\\Gamma(a+b)}{\\Gamma(a)\\Gamma(b)} x^{a-1} (1-x)^{b-1}, $$\n",
    "\n",
    "where $X \\in [0,1]$\n",
    "\n",
    "$$ E[X] = \\frac{a}{a+b} $$\n",
    "$$ Var[X] = \\frac{ab}{(a+b)^2 (a+b+1)}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Note that $ \\int_{x=0}^{1} p(x|a,b) = 1 \\Rightarrow \\int_{x=0}^{1} x^{a-1} (1-x)^{b-1} = \\frac{\\Gamma(a) \\Gamma(b)}{\\Gamma(a+b)}$. This gives us a \"cute\" way to compute the expectation\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "E[X]&=\\int_{x=0}^1 x Beta(x|a,b) dx = \\frac{\\Gamma(a+b)}{\\Gamma(a)\\Gamma(b)}\\int_{x=0}^1 x^{a} (1-x)^{b-1} dx \\\\\n",
    "&= \\frac{\\Gamma(a+b)}{\\Gamma(a)\\Gamma(b)} \\frac{\\Gamma(a+1) \\Gamma(b)}{\\Gamma(a+b+1)} = \\frac{a}{a+b} \n",
    "\\end{align} $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.text.Text at 0x7f5441c85250>"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEKCAYAAAAfGVI8AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XmclXP/x/HXlDZaJLmVLVJUSKRNN4PKkkhlSbZSd7Il\nS7fljm4SoZSQhEIpRQqhREPa+ylL+6ay3EpatU0z5/fH5xxzZprlzMy5zvc657yfj8f1OOtc12eu\nputzfXcQERERERERERERERERERERERERERERERERERERERGR5FYDyAAWBbfvgblAswh+dhpwRDGP\n9x2wELgx+HlDYEIxjn122M+nAssKGV/OfUcaT2E8DbQqxs8PBNaTdQ7H5vKdEsBzwIxC7PdO4Efg\nB2ASUDX4fnngE6BsEeMVEZ+rAezM8d7VwMoIfjYTqBKF4x0PrALaFWI/kRw7laIlgqL8XpFqAkwu\n5j5mB/dTkC5Anwj3eTawDqgQfP0s8ErY5zcG3xORBFSDgy/MPYC0sNdtsFLCt8A32EVoJHbB/A44\nBrsDHRL83hJgKbmXKnI7HkBHYA5wPnZHCnYnOgG76/0/4FUgJcexbwg+zgIWY3faoZ9PBTYA44P7\nmAHUCvss9L3Q6x+BN8L2fWyO7/0r+HwxMDXHvmYBbwXP0ZLge7n5DLgsj89CLiLrbj98awmUAfYA\n7wfjeA84Lo/9vA2cF+F+AUoGH8tipYx+YT9bFvgVOKqA2EUkDtUADpB1UfgJ2AdcEvy8FlZdVDn4\nuh52QTgUu2CGqlCaAO+G7fdB4MM8jpdbIqgH7CJ7IrgR+DT4vASWCE4Kvg4dOzUYf+himEr2RHCA\nrLvnbliiyvm90Ovvc+w7/HsXYKWWUEnhZuyCH/pOOnBG8PW9ZE+kIYcHf8dDwt6rDTyBJYfRwOW5\n/Fy4GsDHZCWh+7Hkk5sNWKLsBNxTwH5D2gKbgY3AyTk+Gw/cEuF+RCSO1ODgC3NTYFvws9uxC0P4\nHeRG7KIXfsEEu6j1wKoQFgBfRHg8gLrAJrInghrYxWwG8HDwOyHhiWBd2PupZE8E4RfJ0lhiqEDu\niSD0OrdEMAC7YIfbBpwQ/M7qsPcvICuphGuIJZOQw7C7+sODr7/E7rhbkPude17tCtuDcYSrRfb2\ngY3kXSLIbb9dgTU53nsquEkcOqTgr4hkMwdYATTC7sS/AK4L+/x44JccP9MaGIw1UE4ClmN3o5E6\nh4Mvnj9hd6WpwIXAdOAurFok3K589puR43UAu3sPYNVMIaULiK9Eju8TfF0q+HxPjmPk/C5YgikZ\n9rodlmS2YVUv5bFkOB1okEccpwNnYtU+oRhSsN8pXHNgSvD5KcAO7N8xr/3WBKphVX9g1W+vYCXB\nrcH3SgL78/h58bkSrgOQuFM7uH2L3aW2wi4mYFVGi7G66gyyLqAtgI+A4Vh9/lVkv+gVdLz/YL1h\nwi+gPbAL0jSsqmkqVoVEjmPnpz524QToDswE9mKlnOOxnjEpWJVIIJ99TwWuBY4Mvu4M/EH2kkBB\n1mJ3/KF9H4m1RYCdv7lkVcnlJYC1xdQIvu4R3MevOb5XGWvzgMgaeqtj7QKhqq9OWJLaGvadmhSt\n8V18QCUCKUg5rIogpARWnx66yP0LGEfWnWcbYDcwEbuDvAK7e3wnuJ+tWM+Y+yI4XiZ2YX4Qaw9I\nJeuC/CZWVbQU+AvrMjkk+Fno2M+HfT8kEPa4DHgMa1v4H1a3T3Cfw7Guq79h9e4h4b9XILhNDx7r\ny+D52UT2+vy8Ygi3DUtEF2KNxmODv/elWELKJKuaKC8/YqWij7BEuxFraA+ZBdyNtdd0BY4O/n6j\nCtjvTOBJrG3jAFbiaxv2eRmsraVzAfsREZECNCV70om2PsCpHuz3FqydRCSbUlg95dfAPOwuMVwv\n7O5lRnCrHdPoRPxrIHCxR/v24o69AlY1Vs6DfUucuwUYFHxeGSu2h3ubvBumREQkARyG9XIAa2DK\n2dVsKTYYaCZWDyoiIgmqAtaAdl2O9/tgfbFLYXWirWMcl4iIBOXWnzlajsN6WLzEwb0SKmJ9l8G6\nuFUh+5B1atasGVizJmdBQkRECrCGg0d+58urcQT/wPp39+bgJFAJ64N8GJaILsS66WWzZs0aAoGA\ntkCAxx57zHkMftl0LnQudC7y37AxHYXi1TiCh7EL/qPBDWAEdvEfgbULzMDmrZmO9ZsWEREHvEoE\nPYNbXsaS+zzpIiISY5piIg6kpqa6DsE3dC6y6Fxk0bkoHi8bi4srEKzvEhGRCKWkpEAhr+0qEYiI\nJDklAhGRJKdEICKS5JQIRESSnBKBiEiSUyIQEUlySgQiIklOiUBEJMkpEYiIJDklAhGRJOfVpHMi\nEif27oXly2HTJtv+/BMqVICjjoKqVaFWLahc2XWU4iUlApEkEwjAnDnw6afw1Vfw7bdw4olw9NF2\n4T/iCNi1KysxrFwJtWvDBRdAq1Zw0UVQQnUJCUWTzokkia1b4e23YfhwyMiA9u3h/POhaVMrAeRl\n/36YPx9mzIBJk2D7drjrLujcGSpWjF38EpmiTDqnRCCS4HbuhAED4OWX4ZJLoHt3OO88SCnC//5A\nAGbPhiFDYPp0uP9+20qXjn7cUjSafVRE/nbgAAwbZtU6GzfC4sXwzjtWCihKEgD7uXPPhfHjYeFC\nSwr161tpQeKXSgQiCWj5crjxRqvyGTgQGjTw5jiBAEyeDD17QsuWMHQolCvnzbEkMioRiCS5zEx4\n4QX45z/h1lvhiy+8SwJgJYS2bWHJEti9G5o3h/XrvTueeEMlApEEsXUrXH99VqNwrVqxPX4gAM8/\nD888A2PGWO8iiT2VCESS1IoV0Lgx1KkD33wT+yQAVjq4915rh+jUyR4lPmgcgUicmzbN2gP697fq\nINcuvNB6FLVqBenpcPPNriOSgigRiMSxN9+EBx+ECROsS6hfnHYafPkltGhhyaBrV9cRSX6UCETi\n1Msvw9NPW9fNU091Hc3BTj3VYmvRwkYid+niOiLJixKBSBx69lkbI/DVVzY9hF/VqgWff26lleOO\nsy6m4j/qNSQSZ5580noFTZ8Oxx7rOprIzJxpU1p8+aVVG4l31GtIJMG9+CKMGgVpafGTBMDGNQwe\nDJdfDr/95joayUlVQyJxYswYmzNo5kybKTTeXH89rFkDV15pv0OZMq4jkhBVDYnEgY8+gm7drGql\nbl3X0RRdIADt2sEJJ1gJQaJPs4+KJKD586F1a5gyBRo1ch1N8W3datNeDB5s01NIdCkRiCSYDRts\nvYBhw+CKK1xHEz1z59rvs2CBlQ4ketRYLJJAdu6ENm1s2oZESgIATZpA795w7bU24EzcUolAxIcy\nMqzapFo1W1GsqOsH+FlmplV5NWsGffq4jiZxqGpIJEE8/LCtKzxtGpQq5Toa7/z8s7UXfPklnH66\n62gSg6qGRBLA5MkwerStApbISQBsLET//jb9xIEDrqNJXkoEIj6yerV1E50wAapWdR1NbHTtCpUq\n2Upq4oZXiaAU8DbwNTAPaJPj8zbAfGA2oHkJRbAVvtq3h//+19YWSBYpKTBihM2ftHy562iSk1dt\nBLcAZwD3ApWBxUCok1gpYCnQENgNzAIuBzbl2IfaCCSpdO5sjcRvvpmYjcMFeeEFmDjRZixNxt8/\nWvzURjABeDTsGOG1f3WA1cB2IB34BvDRTOoisTdmjPWtHzYseS+Ct99ug83Gj3cdSfLxKhH8BewC\nKmBJ4ZGwzypiSSBkJ1DJozhEfG/tWujVC8aNg8MOcx2NO4ccAkOHwgMPwF9/uY4muXg56dxxwETg\nJWBc2PvbsQQRUgHYmtsO+vbt+/fz1NRUUlNTox2jiFPp6dCxI/znP1C/vuto3DvvPGjeHJ56Cvr1\ncx1NfEhLSyMtLa1Y+/CqEPoPIA24HZiR47NSwBKgMVZymI01HuecnFZtBJLwHnoIfvjBJpVL1iqh\nnH75xZLivHlQs6braOKPnwaUDQGuBlaEvTcCOCz4eDnWhlACeB0Ylss+lAgkoX31lU3NvHhx8nQV\njdTTT8Ps2fDhh64jiT9+SgTRoEQgCWvnTjjjDFtopnVr19H4z759UKcOjBwJ55/vOpr4okQgEie6\ndbO5+V97zXUk/jVmjDUez5mjarPC8FP3URHJw5Qptt7woEGuI/G3jh1h716YNMl1JInPz3lWJQJJ\nOFu2WJXQmDGgTnAF++wz61r7ww/WvVQKphKBiM/dfTdcfbWSQKQuvtjWZ37zTdeRJDaVCERi5OOP\noWdPu7s99FDX0cSPefOgQwdYuRLKlXMdjf+pRCDiUzt22BQKr72mJFBYjRvDOefY9BviDZUIRGLg\ntttsRa5XX3UdSXz67ju49FJYs0algoKoRCDiQ2lpVi30zDOuI4lf9etbyUDdbb2hEoGIh/bssV5C\nAwcm3gL0sfbtt3YOV6+GsmVdR+NfKhGI+Ez//nY3qyRQfGedZesbjxzpOpLEoxKBiEeWLbPZNBcv\nhmOOcR1NYpg/37rfrloFpUu7jsafVCIQ8YlAwBqIH3tMSSCaGjWyOYg0riC6lAhEPDBqlK1B3KOH\n60gST58+NjtpRobrSBKHEoFIlP3xBzz4oHUVLVnSdTSJ59xzbbTxBx+4jiRxqI1AJMq6doXy5WHw\nYNeRJK5Jk2wVs7lzNTNpTpqGWsSxuXOhXTtrKK6klbg9k5FhbQWvvWYN8pJFjcUiDmVk2DQSzz6r\nJOC1kiXhvvvsXEvxKRGIRMnw4VCxoi0/Kd676SZYsACWLnUdSfxT1ZBIFGzaBKedBjNmQL16rqNJ\nHk88AT/9BK+/7joS/1AbgYgjXbpA5co2lYTEzpYtUKuWlQqOPtp1NP6gRCDiwPz50LYtLF9uVUMS\nW7fdBtWq2eA9USIQibnMTGjWzAaO3Xyz62iS05Il0KIFrF+vaSdAvYZEYu7tt206iRtvdB1J8qpX\nz7b33nMdSfxSiUCkiHbsgFNPtcFNjRq5jia5TZ6cNcAs2alEIBJD/frBJZcoCfjB5ZfD779be40U\nnkoEIkWwahU0bQo//qjeKn7x3HO2pOXbb7uOxC01FovEyJVX2uRnvXu7jkRCtm6Fk06y6T2SOTmr\nakgkBqZPt5JAz56uI5FwlSvDNdfYrK9SOCoRiBTCgQO2XOLjj8NVV7mORnL67jto0wbWroVDDnEd\njRsqEYh4bMQIOPJIG0Am/lO/PlSvDp9+6jqS+KISgUiEtm2DU06BadPsgiP+NGoUTJgAU6a4jsQN\nNRaLeOj++2H7disViH/t3g3HHw//939wwgmuo4k9JQIRj6xZA40bq7tovOjZEypUsLEeyUaJQMQj\nHTrAWWfBww+7jkQisXRp1vxDpUq5jia21Fgs4oFvvrERq716uY5EIlW3rk1P/eGHriOJD0oEIvnI\nzIR777V5bMqVcx2NFEb37vDKK66jiA9KBCL5GDvWHjt2dBuHFF779rBoEaxb5zoS//M6ETQGZuTy\nfi/gx+BnM4DaHschUmh79libwMCBUEK3THGnTBlbP3rUKNeR+J+XjcW9gRuAXUCzHJ+9DQwCFuXz\n82osFqcGDLC2gfffdx2JFFVopPG6dVCypOtoYsNvjcWrgXbkHtDZwMPATOBBD2MQKZLNm+HZZ+Hp\np11HIsVRvz4cdZTNDyV58zIRTAQO5PHZWKA7cCHQHGjtYRwihfb449Cpk/U8kfjWpQu88YbrKPzN\n1bRMQ4AdwedTgAbBx2z69u379/PU1FRSU1NjEJokuxUrYNw4m85Y4t/111tbz5YtUKWK62iiLy0t\njbS0tGLtI5J6pJeAEcDiIuy/Bnb33zTsvUrA90BdYDcwHngd+CzHz6qNQJy46ipbkP6BB1xHItFy\nww1wzjnJMXW4V20EHwOPALOBHkDFQsYVupp3BLoB27F2gRnA11jvoZxJQMSJmTOty+Fdd7mORKKp\nSxd4/XXQvWXuCpM1qmJVOlcCE4AngDVeBBWkEoHEVCAATZrA3Xdb+4AkjsxMOPlkGD8eGjZ0HY23\nvCoR1AUGYHfv27DG3RexZCCSMMaPt4VnNHgs8ZQoATffDG++6ToSf4oka3wDvIZd+P8Ke/9OLCF4\nRSUCiZl9+6BOHas+uOAC19GIF9atg0aN4JdfoHRp19F4x6sSwVRgFFlJ4Kngo5dJQCSmhg2zRKAk\nkLhOPNEmo/vkE9eR+E9+WeNWoCtWNbQ0+F4JoDTW3dNrKhFITGzbBrVrw4wZUK+e62jES2+8AR99\nBB984DoS70R7PYIyQDWsx1C/4HczgE3AvqKFWChKBBITvXvDn3/Ca6+5jkS8tmOHrV62erWtPZ2I\nop0IzgEWABeT1QU0Jfh8WhHiKywlAvHc+vW24MwPP9ii55L4OnWy3mGJ2kU42m0EFwYfO4Zt1wUf\nRRJCnz5wxx1KAslEvYcOVlDVUCCP76hqSOLeokVw6aWwapWtbyvJISPDFrWfOjUx24SiXTX0E1lV\nQuECwEmFOUgRKRGIp1q1grZt4fbbXUcisfbggzaAcMAA15FEnxavF4nQtGlWR/zjj8m3uLnAkiVw\n8cWwYUPiLToU7TaCl4KPc3Jss4sSnIhfZGTYhHJPPaUkkKzq1YOqVeGrr1xH4g/5TUP9ePCxI9mr\niPxcihAp0OjR1iZw1VWuIxGXbrjB/hY0iDCyi3pN4BngFOA7bAnKX7wMKkhVQxJ1e/bAKafAu+9C\n06YFf18S16+/wmmn2WPZsq6jiR6vppgYgc019E9s7QCt9SNxa8gQm29GSUCqV7cxJB9/7DoS9yJJ\nBAHgU2ArMBlQrarEpc2b4bnnrG1ABLKqh5JdfsWHs4OPD2ETz30FNAEuAm72OC5Q1ZBEWc+eNi/9\n0KGuIxG/2LEDjjvOZiY94gjX0URHtLuPjuLgcQShKSY6F+YgRaREIFGzerVNK7BsmfUWEQm57jpr\nMO7e3XUk0RGrcQTVgV+L8HOFpUQgUXP11dCggS1iLhLuo4/gmWdsmdJE4FUieAK4DZt++jBgIVZF\n5DUlAomKOXPgmmtgxQo49FDX0Yjf7N8PxxwDCxZAjRquoyk+r3oNXQEcB4wBTsUWmxeJC4EA3H8/\nPPGEkoDkrnRp6NABxo1zHYk7kSSC34C9QEVgNXCCpxGJRNEHH8CuXXDjja4jET/r2BHGjnUdhTuR\nJIKfsdXK/gKeBtTUJnEhPd0mF3v2WShZ0nU04mfNm9viREuWuI7EjUgSQXdgOnA/NqL4ek8jEomS\n4cNtndpWrVxHIn5XooT1HkrWUkEkDQpHYmMJQlNMPA3s9DKoIDUWS5Ft325TSUydCvXru45G4sGi\nRdC+PaxZAylxPKOaV43Fo7C2gT7YesVa20d876mnoHVrJQGJ3JlnQpkyMG+e60hiL7/ZR0MOBYYF\nny8C2noXjkjxrV8PI0bA99+7jkTiSUpKVqNxk1h0kPeR/EoERwBVsHaBltgcQ+cBK2MQl0iRPfII\n3Hmn9Q0XKYyOHW1m2gMHXEcSW/nVI83I57NYzOCtNgIptIUL4YorYOVKKF/edTQSj845B/r3h5Yt\nXUdSNF5OMVEFW5dgHbC5cGEVmRKBFEogYHPGdOoE3bq5jkbi1fPP2xKmr7/uOpKi8aqx+BpsicqH\ngbmAhuaIL334IfzxB3SOxZSIkrCuuQYmTYJ9+1xHEjuRJIJ7gbOwRuIzgZ6eRiRSBPv32zrEgwbB\nIZF0gRDJwzHH2Mpl06a5jiR2IkkEGcCu4POdwB7vwhEpmmHD4OSTNXhMouO665Jr7qFI6pFGA78D\nM7HlKqsAt3gYU4jaCCQif/4Jp54KaWlQt67raCQRbNoEtWvbesbxNlmhV20EXbBG4hbAWkDNcOIr\njz9us0cqCUi0HHUUNG4MU6a4jiQ2Iska0wAXBW6VCKRAK1dCs2awdKn95xWJlpEjbWH79993HUnh\neNV99F3gHWAFkBl8LxaDypQIpEBXXmmJ4N//dh2JJJpt2+CEE2DjRqhY0XU0kfMiEVTCZh7dleN9\nDSgT56ZPt3Vmly61OWJEou3KK20iuptuch1J5KLdRnAnsBioDDyDXfxDW6Qak/sI5TbAfGA20LUQ\n+xMBbAqAXr1srQElAfFKsvQeyi9rzAHOx1YmGw1cUsh99wZuwEoTzcLeLwUsBRoCu4FZwOXYzKbh\nVCKQPA0fbpODzZgR31MGi7/t2mXjCtauhSpVXEcTmWiXCPYA+4E/sIt3Ya0G2uUSUJ3gZ9uBdOAb\nbDI7kYhs2waPPWZTASgJiJfKl4eLL4aJE11H4q38EkH4f7FIupnmNBHIbQ6/ilgSCNmJtUWIRKRf\nP2jTBho0cB2JJINrr7UZSRNZfoPx62G9hVKAukBoEbcAxVuucjtQIex1BWBrbl/s27fv389TU1NJ\nTU0txmElEaxcCaNG2aRgIrFw2WVw663w++/wj3+4juZgaWlppKWlFWsf+RWsU7GLfs7vBICvItx/\nDSyBNA17rxSwBGtI/gtrMG4D/JbzOGojkJxat7YZRu+/33UkkkxuuAGaNoU77nAdScGK0kaQX4kg\nrTjBhAldzTsC5YER2ER2U7Eqp9c5OAmIHGTKFFi9Gj74wHUkkmyuuw6eeSY+EkFR+LmpTSUC+dv+\n/TYj5ODBVlQXiaX9+6FaNfjuOzj2WNfR5M+ruYZEnBsyBGrVUhIQN0qXtsFlEya4jsQbKhGI7/3v\nf1YamD3bZoQUcWHaNOjTB+bNcx1J/rxcqtIFJQIBbHh/tWowYIDrSCSZHTgA1avD3Llw0kmuo8mb\nqoYk4XzzDXz5pd2Jibh0yCE279D48a4jiT4lAvGtAwesl8Zzz9kITxHXEnVwmRKB+NYrr9j8Ltde\n6zoSEfPPf9rAshUrXEcSXWojEF/atAnq1bPlJ+vVcx2NSJaePe0G5dFHXUeSOzUWS8Lo3BkqV4ZB\ng1xHIpLdnDk25cSSJf6c9DDaI4tFnJg5Ez7/HJYtcx2JyMGaNIHdu22+q9NPdx1NdKiNQHwlPR1u\nv92mmK5QoeDvi8RaSgpcc01iLVijRCC+MniwLQTSoYPrSETyFlq5LFFqr31Yw/U3tREkmQ0b4Kyz\nbMDOySe7jkYkb4GAjXIfOxYaNnQdTXYaUCZxrWdPuPtuJQHxv5SUxFrPWCUC8YXJk+GBB+D776Fs\nWdfRiBRsyRK45BJYvx5K+OiWWiUCiUs7dsBdd8GrryoJSPyoV8+6OM+a5TqS4lMiEOf+8x9o2RK0\nEqnEm0SpHlLVkDg1bx60bWvF7COOcB2NSOGsXWvjCn791Sal8wNVDUlcSU+Hf/0LBg5UEpD4dNJJ\ncOKJNkNuPFMiEGeeew6OPho6dnQdiUjRdexo3UjjmaqGxInly6F5c1i4EGrUcB2NSNH9+qutoPfb\nb1CmjOtoVDUkcSIjA7p0gf/+V0lA4l/16lC/Pnz6qetIik6JQGLuxRetYa1HD9eRiERHvPceUtWQ\nxNTatdCokRail8SyZYs1HP/8s/vJElU1JL6WmQndukHv3koCkliqVLHVyyZPdh1J0SgRSMy8/LLN\n437vva4jEYm+Tp1gzBjXURSNqoYkJlatgmbNbDi+SgOSiP76y6ZQX7kSjjrKXRyqGhJfysiAW26x\nqSSUBCRRHXYYXH45jB/vOpLCUyIQzw0aBKVK2cRyIoksXquHVDUknvrhB7jwQpg/34biiySy9HSr\nHpo713oRuaCqIfGVvXvh+uvhmWeUBCQ5lCpl6xm/847rSApHJQLxTK9esHEjTJhgKzqJJIPZs+HW\nW2HpUjd/90UpEfhk4lRJNNOmwXvvweLFSgKSXJo2hX37YNEiW4M7HqhqSKLujz9sLqGRI22gjUgy\nSUmxRuPRo11HEjk/36upaigOBQK20MzJJ9s6AyLJaOVKOO88m3Ii1gvWqLFYnHvhBZuW96mnXEci\n4k7t2tZBYto015FERolAombhQnjySZuFsXRp19GIuHXTTfDWW66jiIxXVUMlgJeBM4B9QFdgTdjn\nvYBbgc3B192BlTn2oaqhOLJjhzWM9e9v3edEkt2ff1qpYP16OPzw2B3XT1VDbYHSQDPgQSBnbfFZ\nwI3ABcEtZxKQOBII2NrDLVooCYiEHHGE/Z947z3XkRTMq0RwLvBZ8Pk8oGGOz88GHgZmYolC4tiL\nL8KKFfD8864jEfGXeKke8ioRVAR2hL3OyHGssVh10IVAc6C1R3GIx2bPhn794P33oVw519GI+Mul\nl8KyZbYgk5951bFpBxC+Tk8JIDPs9RCyEsUUoEHwMZu+ffv+/Tw1NZXU1NQohynF8fvvVhX0xhvu\n5lUR8bPSpeHaa21MwaOPenOMtLQ00tLSirUPrxqL2wFtgM5AE6APWXf9lYDvgbrAbmA88DpZVUkh\naiz2sQMHoGVLW5Xp8cddRyPiXwsWWDJYvRpKxKCfpp8aiz8A9gKzsIbiXkBHoBuwHWsXmAF8DfzI\nwUlAfO6BB6BMGXjsMdeRiPhbw4ZQvjwU86bdUxpZLIX2xhswYADMmxfbbnEi8WrIEJuKPRZrFRSl\nRKBEIIUyaxZcdRXMnAmnnOI6GpH4sGUL1KwJ69ZB5creHstPVUOSgDZsgKuvhjffVBIQKYwqVeDi\ni/27ToESgURkxw644gq47z7rEicihXPrrfD6666jyJ2qhqRA6enQpg3UqAHDhml9AZGiyMy0KScm\nTYIGDbw7jqqGJOoCAbjjDuv29uKLSgIiRVWiBHTu7M9SgZ//W6tE4ANPPw3vvgtffw0VKhT8fRHJ\n2/r1Njnjzz97NxJfJQKJqlGjrCro44+VBESi4YQToHFju7nyEyUCydXkyfDQQ7awxjHHuI5GJHHc\nfju8/LLrKLJTIpCDpKVBt25WElA3UZHouvRS2LTJpp7wCyUCyWbhQptIbvx4OPts19GIJJ6SJaFH\nD3+VCtRYLH/79lu7WxkxwsYMiIg3Nm+2dY1Xr7bBZtGkxmIpssWL4bLL4JVXlAREvFa1qo3NGTnS\ndSRGJQLhu+9s+PtLL0H79q6jEUkOc+dCp06walV0p6dWiUAKbcECSwJDhyoJiMRS48ZQqRJMneo6\nEiWCpPafheqLAAAHTUlEQVTVV9C6tbUJXH2162hEkktKCtx9Nwwa5DoSVQ0lrU8+gVtugXHj4MIL\nXUcjkpz277dlXj/6KHrzD6lqSCLy9ts258mHHyoJiLhUurSVCp57zm0cKhEkkUAA+vWzFcamTIG6\ndV1HJCLbt1up4NtvbQqK4tIKZZKn9HTo3h2+/95GDB99tOuIRCTkgQfs/+jgwcXflxKB5GrzZrj2\nWltAe+xYOOww1xGJSLiff4YzzoA1a4q/lKXaCOQgixbBOedAkybwwQdKAiJ+dOyxNpBz2DA3x1eJ\nIIG98w707Glzmqh7qIi/LVkCF11k006UL1/0/ahqSADYvRvuuQe+/BImTrQip4j43/XXQ7168Mgj\nRd+HEoGwZIm1B9Svb8XMihVdRyQikVq1Cpo2hZUr4YgjirYPtREkscxMW1M4NRXuuw9Gj1YSEIk3\ntWpBu3bw7LOxPa5KBAlg3Tro0gX27bPZDLWYjEj82rgRzjzTSvdF6eatEkGSyciwUkCjRjZn0MyZ\nSgIi8e644+Dmm23wZ6yoRBCn5s2ztU/Ll7c1BOrUcR2RiETL5s1w6qkwfz7UrFm4n1WJIAls2mQj\nhNu2hV69bH1hJQGRxFK1Kvz733azF4v7YSWCOLFrFzz+uM0PVK4cLFsGN9xgU9mKSOLp1Qv+9z8b\nD+Q1JQKf27PH2gFq14bly62oOHgwHH6468hExEulStlaIffdB1u2eHssP99PJnUbwc6dVvc/aJCt\nZPToo3DWWa6jEpFYu/tuqxF4443Ivq8BZQlg3TqbEmLkSGjRAh5+WCODRZLZjh022vitt+CCCwr+\nvhqL49SBA7Y+QNu2NkEc2FrC48YpCYgku4oVYfhwuPFG+OUXb46hEoEjgYCtDTB6tG01atiqYZ06\naYZQETlY//42g/DXX1uHkbyoasjnMjPtTn/iRHj/fRsQdt11tnawBoKJSH4CAbtRBBgzJu8eg0oE\nPrRhg80COnUqfP45HHUUXHUVtG9vi1Wr+6eIRGrPHjj/fLjyyrxnKPVTG0EJ4BVgNjADyDk2rg0w\nP/h5V49iiLm//oLZs2HoUJtO9vjjrc5/yhSbZ/zbb2HpUnjySesBFGkSSEtL8zTueKJzkUXnIkuy\nnIty5WDSJOtBdM89sH9/dPbrVSJoC5QGmgEPAgPDPisFDAJaAucD/wKO8iiOqAsErE/vwoU20KNP\nH+jQwUb3Vq1qC8H88AO0bAlffGEDQiZMgK5dLTEURbL8kUdC5yKLzkWWZDoX1avb9WftWptt+Oef\ni7/PQ4q/i1ydC3wWfD4PaBj2WR1gNbA9+Pob4DzgPY9iyVdmpi3ksnMnbN9uXbX+/BP++MO2zZvh\n11/ht9/s8aefoGRJOPFEOPlkSwAdOlg//zp1bBCIiIiXKle2ksGAAVbrcM890KqVrUNSFF4lgorA\njrDXGVjpIzP42fawz3YClXLbSZs29pizqSD8dSCQfcvMtC0jwx7T0617Znq6FaP27bNt715LAHv3\nQtmyUKmSddOqWNEWhDjyyKyteXPLwtWqWe8ejeoVEddKlICHHrKxBWPGQMeOdhPrJwOB8FVyN4Y9\nPx2YEvZ6ENAul32sBgLatGnTpq1Q22p8oh0wMvi8Cdkv/KWAlUBlrB1hIVAtptGJiIjnUoBhwKzg\nVhvoCHQLfn451mtoIdDDRYAiIiIiIuITSTnmIA8FnYuOwFysp9Uw/D0gsLgKOhchrwJPxSooRwo6\nF+cAXwMzgXFYlWuiKuhcXAUswK4Zt8U2NCcaY+chp7i7brYDQhOsNgYmhX1WCliF9Soqhf1icTPm\noAjyOxflsEagssHX72D/2Ikqv3MR0h37Q+8fq6Acye9cpACLgJOCr7sBiTxhSUF/F+uAw8l+7UhU\nvYHvsf8D4Qp93fTD7KORjjlIJ2vMQaLK71zsBZoGH8G6/u6JXWgxl9+5ABus2AgYTmKXjCD/c1Eb\n2ALcC6RhF8EVsQwuxgr6u0jHzkE57O8iELvQYm41lhhz/v0X+rrph0SQ15iD0GcRjTlIEPmdiwCw\nOfj8LuAwYHrsQou5/M5FNeBR4E4SPwlA/ufiSCwpDgVaABcBEcxaH7fyOxdgXdf/D/gR+CjHdxPN\nROBALu8X+rrph0SwA6gQ9jo08Azslwn/rAKwNUZxuZDfuQi9fg77z94+hnG5kN+56IBdAD8B/g1c\nD9wU0+hiK79zsQW7+1uBXRQ+4+C75ESS37k4Hrs5OAGoAfwD+1tJNoW+bvohEcwCLgs+b4LVeYUs\nB2qRNebgPGBOTKOLrfzOBVg1SBmsQWwviS2/czEUu9hdADyNtZe8FdPoYiu/c7EWKE9Wo+k/sbvh\nRJXfuSiLlRD2YclhE1ZNlGzi8rqpMQdZ8jsXDbA/8hlhW1s3YcZEQX8XITeT+I3FBZ2LC7D68vnA\n8y4CjKGCzkUvrNfQTGxQq1fT6PhFDbIai5P1uikiIiIiIiIiIiIiIiIiIiIiIiIiIiIiIiIiIpH6\nf12NT9WRYyZuAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f5441c4e410>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "N=100\n",
    "a=5;b=3\n",
    "ps=[stats.beta.pdf(i,a,b) for i in np.linspace(0,1,N)]\n",
    "fig, ax = plt.subplots()\n",
    "rects1 = ax.plot(np.linspace(0,1,N),ps)\n",
    "plt.title('Beta Distribution ($a$='+str(a)+',$b$='+str(b)+')')\n",
    "plt.ylabel('Probability')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "maximum at 0.666666666667\n"
     ]
    }
   ],
   "source": [
    "# np.argmax(ps)\n",
    "xs=np.linspace(0,1,N)\n",
    "print (\"maximum at \" + str(xs[np.argmax(ps)]))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Conjugate Prior\n",
    "\n",
    "Note that \n",
    "$Beta(x|a,b)  \\propto x^{a-1} (1-x)^{b-1}$ has the same *form* as Bernoulli and Binomial Distributions. When $Beta$ is used as the prior distribution, the form of the posterior distribution remains the same. Thus we call $Beta$ a **conjugated prior** of Bernoulli and Binomial Distributions."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Posterior probability given *Bernoulli* likelihood and Beta Prior\n",
    "\n",
    "Upon observing $x$, we estimate $p$ by \n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "p(p|x,a,b) = & Const1 \\cdot Beta(p|a,b)Bern(x|p) \\\\\n",
    "= & Const2 \\cdot p^{a-1+x} (1-p)^{b-1+1-x} \\\\ = &Beta(p|\\tilde{a},\\tilde{b})\n",
    "\\end{align} \n",
    "$$\n",
    "\n",
    "where the \"constants\" do not depend on $p$ and \n",
    "\n",
    "<font color=\"red\"> $\\tilde{a} = a+x$ and $\\tilde{b}= b+1-x$</font> "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Posterior probability given *Binomial* likelihood and Beta Prior\n",
    "\n",
    "Upon observing $x$, we estimate $p$ by \n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "p(p|x,a,b) = & Const1 \\cdot Beta(p|a,b)Bin(x|p) \\\\\n",
    "= & Const2 \\cdot p^{a-1+x} (1-p)^{b-1+N-x} \\\\ = &Beta(p|\\tilde{a},\\tilde{b})\n",
    "\\end{align} \n",
    "$$\n",
    "\n",
    "where the \"constants\" do not depend on $p$ and \n",
    "\n",
    "<font color=\"red\"> $\\tilde{a} = a+x$ and $\\tilde{b} = b+N-x$ </font> "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Example:\n",
    "\n",
    "Say a coin with unknown probability of head which is modeled by a Beta distribution with $a=2$ and $b=2$. After flipping 4 times and getting 3 heads, what is the posterior probability distribution of getting a head?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "### Solution:\n",
    "\n",
    "$Beta(p|\\tilde{a},\\tilde{b})$ where the updated $\\tilde{a}= a+3 = 5$ and $\\tilde{b}=b+4-1=3$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "What is the probability of getting a head?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "### Solution:\n",
    "\n",
    "$E[p|\\tilde{a},\\tilde{b}] = \\frac{\\tilde{a}}{\\tilde{a} + \\tilde{b}} = 0.625$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Multinomial Distribution\n",
    "\n",
    "$$ Mult(x_1,\\cdots,x_n|p_1,\\cdots,p_n)={N \\choose x_1 x_2 \\cdots x_n} p_1^{x_1} p_2^{x_2} \\cdots p_n^{x_n},$$ \n",
    "\n",
    "where $p_1 + p_2 + \\cdots + p_n = 1$ and $x_1 + x_2 +\\cdots + x_n= N$.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Dirichlet Distribution\n",
    "\n",
    "   $$\n",
    "   \\begin{align*}\n",
    "  & Dir(x_1,\\cdots,x_n|\\alpha_1,\\cdots,\\alpha_n) \\\\ =& \\frac{\\Gamma( \\alpha_1 + \\cdots + \\alpha_n)}\n",
    "   {\\Gamma(\\alpha_1) \\Gamma(\\alpha_2) \\cdots \\Gamma(\\alpha_n)} x_1^{\\alpha_1-1} x_2^{\\alpha_2-1} \\cdots x_n^{\\alpha_n-1}\n",
    "      \\end{align*}\n",
    "$$\n",
    "   \n",
    "   Just as Beta distribution being the conjugate prior of binormial distribution. Dirichlet distribution is the conjugate prior of the multinomial distribution "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Posterior probability given Multinomial likelihood and *Dirichlet* prior\n",
    "\n",
    "Upon observing $x_1,\\cdots,x_n$, we estimate $p_1,\\cdots,p_n$ by \n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "& p(p_1,\\cdots,p_n|x_1,\\cdots,x_n,\\alpha_1,\\cdots,\\alpha_n) \\\\\n",
    "= & Const1 \\cdot Dir(p_1,\\cdots,p_n|\\alpha_1,\\cdots,\\alpha_n)Mult(x_1,\\cdots,x_n|p_1,\\cdots,p_n) \\\\\n",
    "= & Const2 \\cdot p_1^{x_1+\\alpha_1} \\cdots  p_n^{x_n+\\alpha_n} \\\\ \n",
    "= &Dir(p_1,\\cdots,p_n|\\tilde{p}_1,\\cdots,\\tilde{p}_n)\n",
    "\\end{align} \n",
    "$$\n",
    "\n",
    "where \n",
    "<font color=\"red\"> $\\tilde{p}_1 = x_1+\\alpha_1, \\cdots, \\tilde{p}_n = x_n + \\alpha_n$ </font> "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/usr/local/lib/python2.7/dist-packages/matplotlib/font_manager.py:273: UserWarning: Matplotlib is building the font cache using fc-list. This may take a moment.\n",
      "  warnings.warn('Matplotlib is building the font cache using fc-list. This may take a moment.')\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "from numpy import *\n",
    "from matplotlib import pyplot as plt\n",
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "from matplotlib.patches import FancyArrowPatch\n",
    "from mpl_toolkits.mplot3d import proj3d\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.tri as tri\n",
    "\n",
    "class Arrow3D(FancyArrowPatch):\n",
    "    def __init__(self, xs, ys, zs, *args, **kwargs):\n",
    "        FancyArrowPatch.__init__(self, (0,0), (0,0), *args, **kwargs)\n",
    "        self._verts3d = xs, ys, zs\n",
    "\n",
    "    def draw(self, renderer):\n",
    "        xs3d, ys3d, zs3d = self._verts3d\n",
    "        xs, ys, zs = proj3d.proj_transform(xs3d, ys3d, zs3d, renderer.M)\n",
    "        self.set_positions((xs[0],ys[0]),(xs[1],ys[1]))\n",
    "        FancyArrowPatch.draw(self, renderer)\n",
    "\n",
    "def mapto3dprob(x,y):\n",
    "    v1=np.array([1,-1,0])/(2**0.5)\n",
    "    vh2=np.array([-1,0,1])\n",
    "    v2=vh2-v1*(np.dot(v1,vh2))\n",
    "    v2=v2/np.sqrt(np.dot(v2,v2))\n",
    "    return (2**0.5)*v1*x+(2**0.5)*v2*y+np.array([0,1,0])\n",
    "\n",
    "def draw_pdf_contours(dist, alpha, nlevels=200, subdiv=8, **kwargs):\n",
    "    import math\n",
    "\n",
    "    refiner = tri.UniformTriRefiner(triangle)\n",
    "    trimesh = refiner.refine_triangulation(subdiv=subdiv)\n",
    "    pvals = [dist.pdf(mapto3dprob(xy[0],xy[1]),alpha) for xy in zip(trimesh.x, trimesh.y)]\n",
    "\n",
    "    plt.tricontourf(trimesh, pvals, nlevels, **kwargs)\n",
    "    plt.axis('equal')\n",
    "    plt.xlim(0, 1)\n",
    "    plt.ylim(0, 0.75**0.5)\n",
    "    plt.axis('off')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXEAAAD7CAYAAACc26SuAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztvbuWI0uWHGqoiMIrs07XnNM9RyAFKve/rkCFKgWqXKOO\nMOooI1C9X0ZS4HTzVBWQQAGdV/CwcPMd2yMC+UbmtrWwEA+PwCtgMNg2d1/c398jEAgEAteJT6/9\nBAKBQCDwcASJBwKBwBUjSDwQCASuGEHigUAgcMUIEg8EAoErRpB4IBAIXDGCxAMBAIvFP0XWNnCV\nCBIPfHiQwIPIA9eIIPFAIBC4YgSJBz40rPoONR64NgSJBwKBwBUjSDzwYVFT3aHGA9eE9rWfQCDw\nutg42/Yv/iwCgYcilHjgQyKpbY/AAWATajxwNVjEULSBj4Yhgesf0pMs73F//98XL/S0AoEHIeyU\nwAcECdy7/LnthLpSDwTeDsJOCXwoZJtkSr+0pn0g8DYRSjzwwbBBvuzHlPa+axdqPPC2EUo88GGw\nWPyzqGol589yG+4vjwsE3haCxAMfAqWNQoK2xG23qWoPBN4mgsQDHwQkZCVwdNvsTfen40KNB94q\nImIYePdIBGwJfI7CPgH42S3vAZxwf//fInIYeFOI/4qBDwYl8DmFTSATeSDw9hBKPPCuUapwEri1\nVDyUCjythxoPvD2EEg+8c0wReG3slM9IxL1BOZZKjKsSeFsIJR54t0gq/AtKAvfIWxW5Wid72aaK\n/Fuo8cCbQaRTAu8S2UapEbiXDYfZrqpd44aRVgm8HYSdEninUOuEpG1J+SHgAFlhqwTeBkKJB94d\nkkpWH5ykbRW2lxHXrLin3vP5Qo0H3gJCiQfeIayNYsnYdugZgyr2PcpCZ3x9Aq+PKGwG3hUWi3+5\nz8RN8v4FpfoeS6YQWtRUC+UP5ELnHlHkDLw2QkoE3hnU/+Y91bcl77HL36ZXGDfU9MoJKf0SCLwe\nQokH3g2SCv8FiZypvrluvfA5hU3GCrmc4oVZif/R7f8D9/f/NdR44FUQJB54F8iZcBI3lThz4nYc\n8Tl+OMkaSGRtifwPhK0SeG2EnRJ4J1C1XSPwSzv7cB/JnMd8QSLyX7r1E+KrFHgthBIPXD2yjULS\nphJXP9zLigO+raJFTa5rr8293P5A2CqB10TIh8A7gBYupwh8Tkcf7SRkO/WoIgcy0f8x47yBwNMj\nOvsErhpJhauNosTN+1+620N6ampMkT8Qen790fjcPZ9A4OUQSjxw5dgg2yjs4MNECgncIe85pkdP\nx3osffB9t/yHeexQ44GXRSjxwNVisfjX+1Ip2xELHQJfdLdGbrdy0+1s28N2HqKF8wWq1EONB14S\nUdgMXCXKYqZaHRozdAi8KTe5/0VPsrwHcIaocm7UiOE3lEXOb4giZ+ClEHZK4EqhnrQtZhoCt+TN\nq349cvo7eRjAkDk3asTwZ7fOad3CVgm8DMJOCVwd8vgoJGuSpnrjDoG3SJbJurvngIXr7qaDGNp2\nG2SLBUA5KiL/DeiPShQ5Ay+DUOKBK4SnvG3EEJnAObyJEjVkmwVVeIskstfdNiYOB4r8D3l8jq8S\nY6oEXgbhiQeuCrmY+TtKH1xtFSQCp4q25E3iHpMw9MXvZP3UrZ8AfIf45NYP3wP436Bnfn//X8Ib\nDzwbQokHrga5mGltFL2hTuBcBoBVd+8NofKz23/ozkPyBrIqv4UQOdU3h62lPw6UVdJA4OkRJB64\nImj6hOpbVThKArdEbom7VtgkUXPYlBWytUIyHxD5L4PT0FZZLP71PtR44LkQhc3AVSAXM+l7a8ee\nmQTOGuQawA3Gc+I3XTvWSlcoi6D6OD09a6cffZ6RHQ88H0KJB64EdnxwKvBOVnsEfotM3mskcgbG\n54RgPVKHESdxU5kDpRLvFbmNFdJWUXslEHhahBIPvHmkYiZVN8lbxg2fS+B0YajMba9NXafy/iL7\nSOQ3zuP18UObmsnxw/Q6AoGnRSjxwBVAB6DSXHanfBv4BP4FZUZ8SolrKlCVOPedunOcpR2VOJCC\nKSegVN0cY0VnCQoEng4RMQy8aST1+iuA35CY9FckUvw9NVC1TAJnRx1L4Dq5jxY1lcyVZzVeyHy4\nLn9DSrB8R44dftNzMGb4792O/wPg3yNyGHhSBIkH3ixSMfBXJMImefN+M07ganOQvG0+nGkVD4fu\nXvPilsDP3fbv3U2X74GcH//37vYHErH/e4yrEngyhJ0SeMP4BUMbpfObF90i0ydjBE67Y27M0MYL\nDxj22qS9UkNf6GTvTc4E9MvEgYHAZQgSD7xJZBvFplF+yYVMJW1L4LWIIVB2+vHADj42L86sOO/V\nD/fwDSjHH2eHoBaLxb/chxoPPAWCxANvFDqQlapx5EKmR+ArlINdaUIFKAucQH2KTbbR8VOUzDVa\nSCLXOZbv5Ni+MMu++1HkDDwdgsQDbw6lCtfemcYHXyORtCVwqnAlcLVXgPEr38uKnzG0XmyCRbcX\niRVrq6TsePTkDDwFgsQDbxCasRYCpw9O9X2DTOiWwEnyY8mUuQNg0f/mUOEnpB+HKTuF5/kO4N6z\nVSz7BwKXI0g88KawWPybqHC9oeweXyPwG2T1rdaKEremUsYGwGL7k1nem7Y1sFcn0PnjQ1tlsfi3\n+/v7/xxqPPBgBIkH3gxSpPA/YGijYKiuxwicPSyVwG1xc+zK5z4q7p8YJlTY7nv3uMS5e2w+D0D8\ncR5IFU5bJYqcgYcjSDzwhkDV7dgo7C6vA1VxJFpV5rRPONCVkvdD7RSOZ2UTKt+Qi5s3zjnsuf6G\nLnbIHpycpzPGVQk8HEHigTeBspgpN40TelaKJXBv7BRrp2iBUpMqhHarVztF44bEydx70DZ97HAv\nG6LIGXg4gsQDbwQyoJXaKDpU7A2mCdwmU7SjjzeS4Rj23fnZM5N2CrrtNmrozch2RlbpRbqQoxt+\nQRQ5A49BkHjg1VEWM504Icf/ZodNS+BeMkW9cB1HBRBbZWTIibtF2SuTA1/p5BD019nGEvke+UeE\nnvoXGFuFg63sQ40HHoQg8cCrYljMFBtFFTdVNYeI9QichU0tghYFzo60W/FLdJk4NZl8T510bxfl\nKIaqwtVeIZEfUBK6nauzt1XSPJxZmQcClyFIPPDK8IqZGHapv0EeXpakbgubJHqSvJI3ybq7X645\nwhXQCJGfT9kkP96tUnuS+qnJZE4wB66FTS8/rpbKwFbJ2fGIHAYuRZB44NXgFzM3pZJWH1yVuS10\ncjAsGy+87Ri3PWO5PvSErSReA9ueT00idABYN8lq4eOQ/++8MwjY7idyB6K/oesEREY/AfgjbJXA\nRQgSD7winGImbRQS8hh51wj8C7L6FvImcbck8k/H6jM7/n2J9naHU6fMm/acybxtsipnzHAM2g3/\njGS1DGwVzY5HkTMwH0HigVdBUuEcJ1yKmTWyph/uEfhgBMPTgLzb9tyTdiMZwsbJBp7RYvMpEel5\nmUi8bc84nZohmaMte3R6oNXNQifHJu9tFc2Op1uo8cBcBIkHXglOJlxjguqD17LhlsBvs/re3O4K\n8iZxk7RXyCpcSf2MBuj2HbBEg1NP6udlU5D5/vs2e+UY4VvNmaN77lrD/F8wtgonkwgEphEkHnhx\npEjhf8CgmKlxQvXBNYliCbxfTgS+vN316nu73KPBuSBuEvYKMzzxru0ZTU/ozadzr86Bzi//vgVu\nK0R+hzwOCx/yjNJL/wrgr0CZHf89ipyBWQgSD7woymKmDDerKtv64CRwErpH4OsjlusDNre7Qn03\nOPXkTeImOS8x4oljiS12nTJPxyQyT9gu970qB7okC5YAFmUv0FtnmVYKkP3xfshazY6HGg9MI0g8\n8MLYIE16bDLhY164DnjFGKEhcNonqr6VvBuce9L2LBXFAUtsujbn7itCEs9kvsyq/DYVPpOL3hE5\nJ1VWEqf61pEQ6Y//BcZWSUXOxeJ/3N/f/7+hxgNVBIkHXgy5mGlslK8obRS1TTxSX2NA4NvbXa++\nt9gNyNuzVIBhYfOMFltJh5CwN503bsl8hy22y9y+IHILPhRjhkD2xzlsbdEJKCVWosgZGEOQeOAF\noVOt/Y4iE642ivrgVYU+JPAlDljh2JM4yVuJe0qFQ7b3Prh4442Q+Q5bbLHDAUtgmaOLVSLXgbCs\nP05bpUir0GcJWyVQR5B44EXgFjMXSDaCR9K2kDlC4JtPezQ4YYt9r7432PXk/aSFzV6Jn/p9/Q/E\np3Nvn7hErlHDP3fL9Mdpq/SdgH5HFDkDcxAkHngh0D5hMRPDXpmaA/cKmRUCp/LeYleo72VH1M9R\n2KSSP8tXaIdNypfXiJwWCu9vkdU3bZWvqAxZG2o84CNIPPDsWCz+x31vn2gm/Ct8G4XbbnExgXN9\nTmETANpi8PCEE5p5hc1umftT+wb4BJ/IlbxT4+yHq61y1x14Iol/QypyhhoPDBEkHnhWlD0z/xHu\nCIXWRvG61lcInKRN+4S+eK2w2QqZE0rqZ7RYyr4jlj2pe4XNQ9da7RUAPpHvO/5Vb5z+uLVV7iBj\nq0SRM1BHkHjgmaHFTOla/xV1G8US+BdUCXyD7IdbX5zk3ToqfG5hk+Sd9uQenKrC03lZ1EzFTgBD\nIv/S9fpREmdXfGurcNu3DdKPYBQ5Az6CxAPPhlzMlPFRammUrwD+hNJS6ZX5aZLANVbI7S3q8cLl\nAwqbJHQlc7bbYdNHE9UnL4j81I21Ykc0rNkqfVqFajy9h6HGA4og8cAzgsXML0gdfFCOUMjbbxjx\nwU9YVghc/XBb2KQyXyEPekXibjGvsNlih5MpbKI7/xltKmRiX/rhQL+9xyfgtG5wvm1wxBY4dV87\nZsNt703aKZxgou+SzyJnTB4RyAgSDzwLsgqX4WZJ0r8h2yhfUfrgX2X59h7L2x02EwSuvriqbyrv\nJQ6upQLML2y2OOOEBkessMIRBwBb7HtVntrKmCpF33v0PTsB4Hj6gr7QqX2N/ozSVvkN0iU/ipwB\nH0HigWfCr8jFzN/LrvW0UJgR93xwLWQuhwSuHXusL+6R93hhsxzFcKywSTLXvDiJXm2UQ3EW9B2C\ncCvjrOwXQ3/8K0pb5QdMl3xO5xbeeCAhSDzw5EiRwv+IXMxE7lpPpc3CJkmdy7aQudxj2XfgyQTO\nAuamu1cbxZJ3qcLnd/Y5YIUN9r2qVi9cyZx++MBGMecCACyB4+2qLHQCpSdOBX5CJvVbmCLnPsZV\nCQAIEg88MdzxUbSYqWkUtVE2ui0VMr/+8reesKnEawTu+eKZyDNxq+oeK24esepTJqmwmQndkjnb\npOOWg3Mduo5DxPY2LReFTvXHGTukraLplShyBgyCxANPjF+Rc+G/lZnwP2PcRjGFTCXwTYXAmQ9X\nPzwr87KXpi1sKqFbsM2wsJmXlbBzO2+moFKdnz+1wG03FjkLnXa2H/6g0VZhl/y/Ir2v0ZMz0CFI\nPPBkcIuZ1jrRsVCsjdIVMr98/dYVMjOBJ/I+Dgic9gn9cBI5lTmA3l7x4oW1wiZx7GyQsrB5wAGr\nolCqccMpsFfn6bZT76cv2R+v2So3kC75Ou54jKvy0REkHnhCfEFRzFQb5Sv8NIraKJ0PXiZRdoWd\nMkbg2y7ulzv7ZPK+JF5IfW1jhkesnMImh6TNOfEpUNEzsXK8WyV/3LNVfiCJ7i/IXfK/I4qcgR5B\n4oEngVvMZCceTaSojfIVyWL5ityhZzmfwOl/W1/ckvdDC5s2/61q3JI5kP1wT90rqNqZWDnerpI/\n7tkqmlYhgX8F8Fdb5Aw1/lERJB54NEaLmTdIRG079RR58BM2X7+JD07CPo4SOH1w7tOOP2VC5WkL\nm5bMOa74HCuF5wNyYmV7u8v+OL+SaqtQmbvZcRY5f5312IH3hyDxwBOAqlCKmX9Btk5sGoVkLj54\n7pG5N1HCaQJn23bQ4We8sMntFomoE4GPFTaVwNVHT22GBU6CGXIq9/On1vfH2YOTnYC+wsmOs8j5\nLSKHHxRB4oFHoVrMZNd6Ejk98Rk+uEYJawRO/1s7/qj6Hitspm11T5w4ynCz9Ma9wia3A0MLxgMV\n+EaUftERiP44C5lqqxzgZMc5TvvPiBx+QASJBx4Jjo8ixUxaJXZGHs9GcXxwjRZOEfgwH14rbM4d\nxTAjd6fPoxjWCptU97xPCt0vdGabhj09c6Fz4I/vMbRVGEF0i5xhq3w0BIkHHoykwv8Titl6bCac\nyyM2itomW5TRQo/AlbjVPvELm2YUw3NW46tDncgPqyW22OPcpEGv7CiGXmETIHnvZNlkxNGWChxl\nobP3x8dsFebHOVBWUeT8I4qcHwxB4oFHgEPMOsVMzYdP2ijZBx8mU6YJ3C9s5kGwAKA5n3vSboxd\nrevn7huxPR2L9WZ17gm9MWRuC5vqjw8GwnL2aaGT/njVVtG8ON/nosj5K1D5BxB4nwgSDzwIOVLo\nFDPV//4HzLJRykKmRggPVQK3vjjVdz9Jcqe6V4djT9S8X1SCKa1sv1+l9uc2kfq5TQp91RxBMuck\nEFTlAAbKPD3HfaHCUzuJGiLZNWc047bKobt9Q86OnyBFzm+IIufHQpB44GKkSOF/QrZRTDFTu9Y7\noxMui3FRSh98gz1aQ9QtzpMEngqbMgRtp7ybUyLinrRVhY9FxVfAomvb9mocaAZkjl6VE56NssN2\nkE3nMVro5AQUo7bKXff+ana8L3JSif8RRc4PgiDxwANA8v6Copj5G7L/fYM0U4+OEd7ZKF++fsMS\nB/wD/jbwwVnIzF3tOVPPeaDSrX2yxc4nb510AZgm8pW0abs2Ham3bUnm2WZh8+l0Sn7oROJa6Dyj\nrdsqLGzeIk/hRktlUOT8HdGT82MgSDxwEcop135PG1nAZCacsUI71drtCX/6819HfXD1wkngXg9O\nJXDaJ9vdfkjeStxK2PUY95DAua37tlgy390ssW12hb2iBc4aGEs8iTdeDKxle3N+RSJveuI/kLPj\ngyJnjKvyURAkHrgQ7B0oxUw7Nriq739Ab6Nsvn6b4YOrIvcJXAuaDU7Ynve++rb3QCI6YqyDpQrp\nNTKJK5l3zbY/sipHk73uMSI/oXGjhikFk9Mr51+a0lbprROU2fFBkfMfkRg/1Ph7R5B4YDbK8VGk\nmEnSLoqXcutslDJOWPrgW0PYtEfGCFztk+XdCHmTuEnal9opP5BJ3ZC5qnLgCKyAXZMnh6gRuZ9S\nafPkERBSv92VtsoJ2RPX7HhR5OQoh7+FGn/nCBIPXAD2DPwVbjFTM+GMwlVtlFJ1e7nwOQS+/XHM\n6vsHhuStxP1QO4Wk3iL/IFgyR+p02ZyOwA0Kn9wSeVLhTbEOpEKn648v29JWUU9cs+M/IEqdc3Iy\ndhh4rwgSD8xCUuH/DwY2CouZmgn/ipQJ1049xkaxPriXC2c0j953lcB/oFTf37snfcbQTplL5NYP\n57YVMpmvpf2PpMqbm2Sv0CcnGDVscS56iybi3ndPben64wNb5U/dJMtfUc7+wwGyiiLnr4jI4ftG\nkHhgEjlSyKLZZljM1Ey49tiUNEppo1jizh12+kJlR945mdL9EHQFzFaVN5epvq2tAvi++BTUD4cs\nr7rHacr9CyR7ZYshkXugfUIy9/xx11bh62VWnNO5sXNQX+SMcVXeO4LEAzPASOE/ouiZSR9cJ3n4\nTZa7Tj2ejWJ9cC1k6iz2OhStS+Cqwq2douRtffE5RL6G74cfMCTzG3n8m/T8tkgFT2zrD0H7RBX4\nEUtssE9WSs1W+XOb44Z8fST0osj5O5L3EuOqvFcEiQdGMRgfRXtmkix+Q86Es2u9dOrxbBTrg2sh\nczlIpjgEroT9Hb76tr64JXLAt1T4rVACV0JfO8doW6AnctykLv+eIldvnMXNU6e6gUTw1lbZ326y\nrXK7yHFDfg62yPk/gazGY1yV94gg8cAEGCfsSNwrZtpM+Ff0Y6PUbBTPB1eLJZN89sAHBP4DmaRr\nRc1aMmWstyb3aUFTSdraKLy/QfpBEQ89CfGhtUJy1vV03xb58YGt8mmZOwHdrnKXfM2Oa5HzC7oi\nJwucETl8bwgSD1RRjo/SqXBbzLTjpNygt1G+LL+P2Cj7rsjnJ1NY2Gxw6ouYLoHTD7aqXMl7rMem\nZ6uspV0toWKtFnTrK1nvLBYl8qY59b4/ccQKG+xw6MYqdxMqsnxcrrKtsm/za/+GsshZ9OTkuCq/\nRZHznSFIPDAC/g3vxkfR7vMk7luU46R8QT/E7JiNkombkxsfizFTVjgmUj/vfQIfK2ra4iZkPzHW\n0edHd68EbQmdanyNrL69b5MhctyY8yKPoaL+OK0Wxg7VVtlhk8dW+fJLzotziFotchY9OfO4KoH3\ngyDxgIvF4v/rVLiMj6IKT4uZmgn/eihm6qnZKDYPvjWEPovAtah5163XyLtW0PTI3PbW1OO4bsmc\n+/jc9Nw3aRuLnUrkx27uztRc44V5hENrq/wDmjy2yu02Z8fpifPzcYucaVyVUOPvB0HigQHc8VF0\nmFlaKlrMZNd6yYTXbBQdmVAVeS5sHlNCpcuBgzePwNU+scVNS95TBU27j144kFV3rbjJH5Ab2cbE\niqw3pzQ07nmbiDpbJ03R3Z7q/GgTKt3yFrs0ZO3Xb9ifvubsuKaEWOSkT/5XoBxXJSKH7wFB4gEH\nZnyUL6gXM80Ihexa/wXfqjYK44TaoccWMplEWZC4pwhcvfEaec8dP8VT4kBJ4FzWouadc7xuR+7Z\nCeyL6CFVuPrjHF/l7Ngqe2z6IWv337c5O+4VOZkbP6Ercv4jInL4fhAkHiiQI4UyPoq1TXS2nt5i\nsV3r6zYK44QrlOOj9PNm7vZpLBQlb5KQJfCxVIoWNZ8iJ86OP+qHa5FTM+Pf7cnkPUYm8sN2WRQ6\njyZeCOTYYWmr7Pshaze3u1TkvJlb5GStYx+Rw3eAIPGAgUYKNzkyyDQKUyhazPyL17Xet1GmfHBG\nCXsCPyAr7TEC92wVwPfDx4qaiqmcONetjdLlw10ib9Nt0QJNC2zP+4E/zqJmVuD7ga1y7KbBKLLj\nU0XOGFflXSJIPNAjj4/SJVJqw8zaYqaTCa/ZKFSd7kiFWsikPfId4wSuWXEl8FpO/KFD0XpWii16\n0l6hFz6hyDWxcm5KfxygAvdtlR02w+z41xVwWtWLnNwW46q8KwSJBwBMjI+ixUySuFvMLLvTj9ko\nLGyyd+YKR6wOkkRhZ5UDSkIfI3CqdqAkdU+Bz7VTAL+jj1opkGWe1yPyToUrsbcAVu0R522OFzJa\nSLuJaZXklS+7Nls/O363HBY5mU5hJ6CvMOOqROTwmhEkHujgjI9ix0b5ItuYCf/6rRgnvOwq76dR\nuF8HuqIPDk2j0A4godcIvJYTr6lwm0xRclcF/h35G2LJWq0UkjptlTEi53Z5ncs7gP44RzjkwFia\nDyex77v30suOV4ucfD/pi0fk8N0gSDxgIoW/lJM9aHHMKWYyE57my8xk3uDUE3aDU0/aqsjbzloZ\n+OBeIsUjcM2FewSuRG6z28C4Gidpk6SV0D01XuvFqUTO527ff5T++Fnskw12ha3iFznTiOWzi5wk\n9v8JpH9dJ8RUbteLIPEABpFCHR9F1bh64qaY6WXC8wTGOjqhGdjK88GtIve621sC57otZJ4wVOJA\nSeYeLBlTaTeoq3G24/GE5sVVhd+ht1j6UQ9/aXCWzj+ereJlx/fcz3HHbZGT/2hiXJV3hyDxD45q\npJDjo6gvrp64KWaOZcJ1HwubjBNybsyBD16LF1oC9wqbNfK+JF5I1HpsUplbFa5JFUvkE4/brIDt\nLuXHj53CzoVNTasMs+M7bLIyv92VRU4drpbvzXdUxlUJNX5tCBL/8PgVxZRrqrzVPmEx809wi5m1\nTPiUjdLnwa0PTkV+Z9YvIXBb5OSyYmwALGAYM7Rkbu0U9cN5PFB+0xrkIidHQezspCWAw+qMTbMr\nbBXbCcjaKix0DoqcdrhaFjk5rdt3AN84rkpEDq8RQeIfGOWUazPGR6lMekyrpJYJ1zSKtVEWVnWr\nIidhazKlRuDqncPsA3wVXut6r/434KtxkjafT2v2AZnImdEGyq78Fl1+fPsj2yqaVlFSV1tlj22h\nzosipw5XqxNXc1wVTvF2isjhtSJI/INiGClEmT7R8VGKOTPtbD2ZzD0bZYlDHyO0NkofJ/R8cBL2\nnGSKVd81IlfSnvLE1ee2KRXriasPPiehcoNM/tYfx9BWAVAo8Abb/v0dWClekXPf5o4+f0Ka6Icq\nvB9XJaZyu1YEiX9YmEihHQdFx0eROTOHw8yqpXLsPW8uazHT2ig9gXk++FQyZYzAveImUBD3foYv\nvvFm8PEI3ObFgUzkto0Way26/VO2Cpf5HnNuTlvkLHpyagyT7zXHVLlDTOV2xQgS/4BwI4XWNrGR\nwsGcmfWemTYTbjv1DGwUzwdXsuM6980hcIe8Sdw/Z3a7//kD+Nyp8YLQvaJmjchr1omSPL1x8c89\nW6UR26RW5NSenLtP2zxcLSOH9MTtuCp/QUzldqUIEv+QYCFzJFKoM9fLnJnDSOF4MVPHRmlw6kcn\nHNgo6oN7hU0bI5xD4ELePy9U4T0MMRdk7hU1tb120wfyrEAecQPpR9OxVdgJaE52XHtycrja/e0O\nx798Ac6LiBy+QwSJfzCUU679NpxyrTZzfSVSyL/ztWJm9saP2J5leFlro6gPbguZd2ZZ21tyR96n\nyluJe2wocQset0dH4JbM1V65QUnkJGgWQOlLa6FT/XFt371P7ATEsVWouvlvR4uce6PUOVzt5naX\n5+SMyOG7Q5D4h8PMKdcuiBQ2OGGJQz/AVa2YWXTqmbJR1AdX++SMOoE76pskrMQ9nHe+DhnyO51L\nyHvjpVMskeuPi/batJ1+KrZK+yOPrWKz47RY+P5Toadi57YfrjYih+8bQeIfCE8ZKfSGmS3HCR8W\nM/tZemwHHrVRuM928CGB2yFnDYFb9U3yVuK+RInrcVsA3+66L01NlQMlkVN9k6h1P4lbwQmZaaus\n0tgqWuTXJI1AAAAgAElEQVRc4TAoctqenE8XOTxF5PCNI0j8g6AaKbSjFA6mXPMjhavuNtYzs1rM\n9JQ2VbnnfVtbpVLcVO+bBE4SVuLeX/C+bWSZ59qiosq75eq4KmeU5K1+eCuvievdj5hX5DwUBD7s\nvbnFzo8c3i2B/SL/uHijHBaRw5jK7a0jSPzDYOYohSZqODdSqNaJKnK3mGnz4dZGUcVtyVqJ3CFw\ntU92yOStxP3zge/gBult20GIHOiJeuN54SRtJXfrj3u2yo9yuVbkbDoaX3bjHNrIoar1/e3mAaMc\nRuTwrePTaz+BwPMjRQp/w2NGKVT/24sUrnAYRAurxUyvU48ldCVqtVZUiTsEztNaAv8pt9MFNx7D\n89jzn2B+QKxn7/Uy1XaawNG2J+Qp51jkPKUiZyLseWPW0CtvcMbm0x6b2x1we0qf76b7rFfmOvja\nXR8AtIaSrqPAW0Mo8Q8BxgmdSKHXsUdVeBcp1BTKsptKbSxSOOiZ6RUza93sSWw6TkrFUrEErvYJ\n1fdP2QbMt1M2coy3j4qcRF4ocrVStFenvg4qbU2rAMMip4x0uGqP2G1PWOGII1b9Z2Ejhzts+xEk\nJzsAReTwqhFK/J0jFTOppszEx5w3c6RjD8nYduxhhFBVuI0UFj0z1RLRSJtV5LpsSdvkwKcI3Crv\nPUpin7rt5WaVObftkFV5kUdX68e+jh8YKnRV4CRRfa+6ZVXj+k/Ifg6qxtN47t24NVaN8/O2apzb\nF0jXTTfKZbqeAm8JocTfMVIx83cMIoWc+JjEPbNjj86POVeFL6i6Sc6qyD0bRS0H64NXooQ1Atd1\nQv3wMUW+kbafpS198Z9mnY/F51QtdKo/3jj7qMBZCzCK3EYO02eyR+70c0EHIG/McXdiZY0c/jLy\nrgVeA6HE3zWYROlslLFI4cyxwlvYzj3Djj2D8VG8McE9FW7TKHYf8jbbC5O7SOCqvgFfXcPsswTP\nW80Xt4ocuMAft0Vb64erpXQql5d3SO9vX0A+VP4VlaNHsp6xxCGrcf33ZcfO+a3btgDyv7nfQ42/\nMQSJv1OkIlQ34bGd+LjWsadT4bZjj1oqjBg2XbcSr2PPxZHCO9TJ3C5XfHBL4EBJwGNWisLbZ8mc\n5/GI3C108h7IP0beD1atyHkolxcHYHU49rUKqnH/39Gu/1z4ea5wxHa5x/J2B9zeGxsN2V7RCCoY\nOUzXU/qXF3gLCDvl3aJT3xopnOrY06nw5frQZ8E19z2lwt2OPZr9HosUqoVCwgJK8uu2XULgXAaG\nSnsKnqWi22zhU4udLHRugDI+yIihfW1a5PQihzquitMBiN3sbXf8fVfanNUd/yf8DkDMkUfk8E0i\nlPg7xCBSCAzHCrcq/AbPq8LnRArHipliowAPJ3BrpXixQjhtvfNZS8YqcqBiq3ivb0qNa/FT1Piq\ny4arGl92wx+woMlbnh6v+xFessh5b/6NYajGv/CVRuTwrSGU+LuEiRR6Y4X/GaUK36DoXm9V+NQg\nV0+mwgG/mIlhhx4PlsAt+fL0Y9D9LXzVbhW590UqYofaFV8jhro8pca5fqEa5+BYjWzfY5N+hNcH\n7NfHend8/Xy+Aak7fkQO3xJCib8z5Eihmfi41rHHqPDlp6OrwjnI0oupcGBQzATGbZQpArdKeyxe\nCOcYT4WrAq/544C8D/q6zmb5gd74lBpnh59JNa4dgLTYrWOO90XOiBy+FQSJvyOkYhP/7nYTH491\n7Jkc5OoJVLguP1SFmzRKjcCJGoETXkHTwiNzPZ+1U2pEDlSKnDYuaWsCT5RU4WfGz6sfUbJrSzWO\n9TETuK2b8J7XUX+Nse4SeE0Eib8rmEjhXBW+LlV4HoNjWoVTAa5wrOfC56hwoK7CMRyV0INHrHKK\n0RihvcE5xhI5l8d+EDStUpzkidV4LalysTc+1h2/6ACk3nio8ddEkPg7QRkp7KJgXvaXKnyFWSrc\nJlLGemf2KpwkZJXknERKRYUTUzYKMIwBAj55j8Ejcz1f7TFrahyQ16GvkSfVH68xNe68v83psty4\np8a3t7ukxnmd3KAcimGgxnVi5YgcviaCxN8N+NdWxgpnl3otYg688NOoCreJlLwt+7ANTuUYKVSL\nuvzMKhwYEmuNwGHa1NIptr0lcm6b+kEY9cZtUkXbWDUOuO9vc0pqnJ8RRzWcyo2rGl9+Sh28sL5U\njfO6i8jhayFI/B0gqfCsigD4HXucQa7QnquJFFVv3tgcyVbJIxW6jEjC8lQ4UcuFY74K92wUoE7g\nHmHX9s05nz4Pq8b7Y89GjQNDNX6WA6wa12W5cYTDVTccLf9RzenFOarGGTHUa4nj7fQdgPK/v4gc\nvg6CxN8FzMTHXsceVeEtXBWuHndtjBSqPJ38eHU4JhLXqdbstGsebFFTlSny/RwVDvg2ClFLqIyh\nRuR231SR1E2q6L19/bX3TJW5eZ+bUxpvvMGpt0z035IdeXJpEi2z1bheVy2QO5Tx+gu8NILErxxl\npFA69jxChecxwocqXGftWaKLt51Sgc1VjHNihcDQWhErRTGlwpEPdTPiFmNFTT2XnmfMnvHUeLF/\nrMAJWZ4qcJr3eXHIapyfD++Hox0ezb+u0zCpUlPj1Q5AETl8LQSJXzFypFA69tS611sVvvZz4XaW\nGF3mOgm9wakeKwSmC5qEMqWQmWel1DBntp4xv3ts+5jtoo9fO99ogRN4WIETst4ts8C5wQ4aMZz6\nXO2YKqNqXAufhRrP12AUOV8WQeJXDXqSMlY4Z+xRFb7GUIWbMVJSUWw4a4+dwT4XM51YIVAWNHVd\nCZ3wrBSgILO5Vgrge9dT9sqcc+nxDx2HheeYbakoVI0TToHT6/wz/OxOZlz48rN/vBrnIFmBl0KQ\n+JWiGilcY7YK3y6T+tJZ6z21Zotk1VihvRkrpIdXwHPwGCtlDHNJd0477/GspTLYf6mlgsp2WzRG\nGTfkD7QWozWBtDSffd/rdkyNK4HrdVZ0MovI4UsiSPxqoZFC07HHDjU7osLpidIX5bLt3FOLFTYk\nEqptXa5ZKcTMVAowTs41K2UuWT8lvMe0KRVgpqXC/XMsla4mwc/DK3Darvj2c9eOQINenGNqPCKH\nr4og8StEjhQ6KtwSd4NRFd6rr84XpSrz4mkkAiAV0AD4PTQVSu6EtQxqFgOGfvJDMJYymUKNlD3M\nSakUqFkqXK6legi1VLp1a6logXOs888K+Zgt0vWxXB/8Xpx6fd3A6QAUkcOXRJD4VcJECm33+poK\n70YqXK6zt83OPFaFW+WtBU3XSgGGVopVlJ7K9JYfgFqW+ylh1fRDHmeWv+/9WxlLqQDFe6+WihK3\nLXDarvhMHPXXBnPjVo2vMVONR+TwJRAkfmUoI4WdJ24Huaqp8Nt7LNeHfu5MDmRFBWZjhbWC5qiV\nQlgrRTHmhwvLjQ05+9ZwicIvXteULz4Haqmg3oNTP8tVr8xLNV7ETNtzVuNr1NV4tTt+RA5fAkHi\nV4TRSCFVOImbiRTtZt+e+xnsvc4gtqCpf70BFNlwAPVUike+nh/uLXewRc2nwqW5CW3/ooPvz/XF\nK5YKkC0vErdaZFMx0gZnbD51Bc72nK8xT40PuuNzELaIHL4EgsSvClQ5pmOPp8KVvNcAZ7Bv23NR\nzNIOH7ag6VkpALKVAvhdILWDj0UtH273PQKfzfpzka8+zqU/Dj89W8mmVnTfHFRSKkrc9h+VV+DU\n5bY9p7k4W2Q1vkapxqnQezUOROTw5RAkfiXIkcIHqPBbAO0ZTZvUlRYxSdI65nT+0vtWCoDshwOl\ndeIRs2e1PDE8mqhRx1xKeezxj8ZYXtwqc33fTUol3eu/qTJ6yM9dC9vs/LP8dETTnoHbUxlhHVPj\n0QHoRREkfjXQ6JaocBK5VeGqmlpgebvD9nZXVVw5Tpj/WlsrBUDfwQdA2alnCrWK3gv53laNTxGx\n3f/q8xjOqYiqvQI/pVJ+pqWlYq8HFju3tFR4rakap72iatztABSRw+dCkPgVYBgpnKHCNfp1e0LT\nnrH8lL/UXkFTrRQAbioFQB4rhbB+uCpzmHbEE5K32hqXkPVm5DZ2jEIfz9o4TwavPuD54gKOpQKg\n/9xWXaOapeLFDanG06QRp3SdUY1vUAoHqvFOOGQ1HhMrPyeCxK8CpmMPN2mMkPYJVfhn5L+2Xcpg\nqqDpWSkAhNTTWCk9PKK2eKYCJVEjXI9QH6Km28pxNcJ+tNViOwF5qL3nNnYI9J9Xi2FKhdvzj7hf\n4Ow7/7TnfN3x+pqlxiNy+JwIEn/jyJFC6UhBFc7oIFU4v0hURA2AdYoV2oKmXQby32sAuYhp7gHj\nhz8TPjcz2ph1Emhr9tuEyVwyr6n6z2a/Je45inzO63NxQexQfXH7OXoplckC5/qQuuJrdJVqnPaK\n641HB6DnRJD4G0Y58XFFhfNLREVElcQv2Dr9FdaC5pSVwpuqcsD44YQtanok85DsM4DN2t/Ol6ao\nEacSuUfmYzc4x9Z+OLzt3o9Fi/rrApDI8CGo/Kjq56Xqm58xgN5S0X2epdK056TG+eKUrFns3Mg2\nV43/GkXOJ0aQ+JuGjRQ6Kpz2Cb9Arey7ranw0koB8t9qIH/Zuc8uAxgvao4RtVWg68pyhxbAFpm4\nPbLeOMs1r/oSu8MSf8171x8K+/w+I38s29oDecTtbav9QGrc0/HIV/I5U30DXn7ct1S0wLnkhMr0\nvbWwXlXjZXf8iBw+LYLE3yjKjj0t+u71ngq/QfmFks49XkGzlkoBcvFL/fAWpZWymJtIGUNNcV5g\nXCtp2mV7Kkvkc26oHO/ZKLXlSVxi1NsfwIl/NVrcTA+VD8j/rvIPtipze30M1Lh2/lHhMKnGI3L4\n1AgSf7NglFBUuO1Gb4uY3NerJL+gqVYKgN5KWZrthQ9ui5pPCUtOTd1yIFnXCNojcvXIP2N+isS2\n13N5NoolelXhFps13Nf9HNDPTW0UjRqmez+l4hY4tfMPr73PyGpcB2Mr1Hh0AHpqBIm/QeRIoXTs\n8VS4FjFb5C/PSEETgFHeR5e0C/X92IhJjZyU3RzS/tyRuVoq9nAlU0usqqitz/15xg3OsZ4/bv8F\neKTdW0LrSlFTX39rtj0BuWtxMz2ELVynayN31T8WhF8tcJK0dUYpJXaq8WJwLKrxKHI+BYLE3xgW\ni3+5z0WgigrXbs+fkb9MWtCcYaXoX+hSlecv7rJqfD8AJKO1s43LE4U96417RM79Y2TuFTDH9lvy\nto/FNqr6J71wfpb2PZiLB5C7fp6aWCkHPcvW2qilwhdIBd6gjB6qGm+A/G7xuv7SXe+BxyBI/M2B\nFzglzOeswgeeN8rOPf1sK/eulQLkiXSBrL6sH67bLAbDz7qNKsu1Npa4W/SWilXjxBiRA0PC5TbP\n764Rum1vz2UfzxK4fQx9Pa6V4v24cVtrtj8C5ec8/Mw1dkpLJW2vWCp886SH8Fw13ieuAg9GkPgb\nQir2jKhw7Q2nhSVHhVsrBaj5oV4SZai+3XjhJbAWgV0miVXUuGereESuKtwjYI/QxwqbtePsY3gE\nXiRrPBul9noviRlOkHpT+bFlpx9gmFgByhii9twdWCq8DvV67DuZoaLG9fr+EkXORyJI/E1BI1gj\nKlxtFY0VSkHTWim2Uw8AUVj5m946/ndNlbtQUrHE7Xm8LXwiMmpcN48RObcBQ6L1iHmOL67He2p8\nisAVvQq3O+x2Vev2/bsk0QKM/viWtY+hL87tuq+wVGpq3BJ7ocY/I8dYosj5WASJvxHkSKF2jKio\ncBYxdVkKmgAKKwXIX8yVUd5TBUyP1Lsd49tqRKPEZC0Vq8YdW4Wn9oicf0ZUlQMlKc+NF3rEbUmd\nj6Xndp+ftVHUC6+9D/b9qtkrdrmC1cgPsVfwtBZbacOJpeIVOO31yYJ8/9qUvNM1H2r84bjwNz3w\nfNALu1tWFb5C+gJoEbNS0FQrBagr72aCwKuplBbDjj66bQ3gh1m2RHU3saznRUeAd6nN/i4TJeet\nrGm5p9B4nr8N1BX5xQTuqXO1WawfvuruvUSLwVm2L3HEAUu3nXctqPo+YlV0z6elcrxbpjfgB4YF\nTi43KDul/Q3APdX4Hn6lIjAXQeJvAClq9SuyjdIVe/hF1kGutIjpWSlAb6UAmjDwI2XAiNp2cG6B\ndi6Br5AKoI1stwTfETO+y/l0O0m9KYkcSGTOYucUmWPmfg/2SzJG3lznv4YBgY89iPXIdbnmk88g\n8sOKczL5BJ4PH1or9rpZ4YAdNsZSaYeWygFZhZ+RP3u+xhNFy76/LRb/dn9//58Xo08yMEDYKa+M\nMlLoqHD+TfVUuBY2aaWw6IRhEsUWNYH5fvcBSxxWDgkoiaxkmaqxRUlKvLe2QWu2W1vFWCuaWqEq\n36LkEdKE2iFehHDqpsfrp1R9bPMcXa97yk5RFa5Wiqpyq8yBOtk/AJNZcloqfD5qqbDAyddzgzw8\ncuGN6yfVRuTwAQgl/uogzWjsCvlLQdJW5tCCplopgFgpQz8cmE6hKFgOPWKJjZcppNK2y0rgQFbj\n3OZZKGPbSEzd/o2uiyoHhsocKNX3BpdNalyzUvSx1DoBSv++sFC47tkotpBpVbgStp6vle2QZV2/\nAEscsZMwJ9dryny5PuC4XgKnRWmpnJDVOJX5pBr/Bc8+POY7RJD4KyIVc9i9nkxtVDgVjFcwAswX\nOqVSgGEGHCjzwWOFrincr4CFZsWbblmtHU6MzO3A0PP2yHoOuXfY2H3IZG7pQEmdw3mM0YX3xdhW\n9lfJmw1rBD1W5JTCbrEOsw3O9g73DybyZJkQtFBWOGKPbeGLp8c/pyeh/5q869VaLt+AoRr/icXi\nX+/v7/9L2CozEST+quCFO0OF38q6/nVVMui+VOqHE48hbcW5leyx/TE5oPTBIfs5wl6NuKeIHM4y\ngM0NgK7z4P4O+HyTJyHeC+E/VuPZL4oSt65XbSSPyMeKnN5xJESYbfoE9TNB+rzOzeU9hEjYus7i\nJpX6wBfna1sh1ThsgfMW5dDFexg1/rO7fbv4+X5kBIm/ErIK14BaVzJzs9+yrtm5XuXmaCFQkrZN\nmVyU++5wRluSgVfMJIFb5c1tStC6PkXkqGw3SZaNIXgldGJvjpkDOxjXgLiB8RqAV7AcI/C1s1/J\ne222aVrFknqH8xN81W1xE0jX2/77Nl1/zSJbfJ6lckZ5PQ/UOL8LX0KNX4Ag8VeDp8I36RPRThJa\n0NS/pqq4xA8HbFzscR7jEauiEHpYLdGcjimhMuWDk3iVSNVq4foYkbfdubh+I8ffyLml7Ua28/du\nL8ROWIJXeINUDUZWtK+5Rt66POWPq42i6+qFW2Xeyrm5HUmFazH6jKbrdznPZ0mRRL8tlTmAYdSw\nZqlo3JDXeVWNX1K5+NgIEn8FJBX+O0ov3FHhtqBJKwUYZoU7K6X3KVHLfs8j9SOWaHDGtnOT0/oJ\nW+yTpUJfnIpc1ThQKm1LiI3sA3wCZzsl/DEC9wgd6H9kNiRveUtmxw3t89dvjUfcdt3zx2tE3sqy\nJXBL7Nynj9s9t/tVzoh78cJj15H+Uuj1U1xf6ot/h2+pUJn/QCLwE0bU+CbU+EwEib8wUoSKmfAW\nAxWuX04taALZSqE3ruoO6IuagE/gU8hplBVaKQWe0PTno6VybrMNOlDjVOl8zkrEMMdY20SJ21Pf\nwDSB67o+BlBe8ZdWNmtKfIrIp4qb1i65gU/gltgh+7ita08CVwuMczopdEI+i1quvOjhKcX04jrd\nd/ffUP7o8L34KesLAPe+Gl8s/uX+/v6/BpGPIEj8xfGLuRkVbgua1kqBcw/0SnwKZ7S9quYQpA1O\n7pdYv9ycdRM44rBaYns65pSKqvE7DONtVnkTnOSZBH2WdmP2idoscwhcLJQHDY1u3xrPC/eIm9vH\n7BWrvj0LxaZZNIvtXB9MpbCTDz9lgj/WFoyTHrshscYwEAlU4jytRl/PKAveqs7XSKp8oMZJ6BE5\nnEKQ+Asi2Si2q0hH5AuUhUySglopgE/g67J/xGN8cHrgJzRYAjhg1dH3rif/VXPsld5AjVsPnPts\nYRLSvlbMVIKf4YWPEvhjJyWaQ+SXeONesXKKwHUwKfXI1QvvVLimUqiorR9Oj/yhGFxn63vg26L8\n98h77eZ6QPo8VJ27ajxd/KHGxxEk/qLQAfGNCte/y1RY1kqxfjgwMHYZL7wUQw88kfkG+0K5ndHi\nAKBZdZNMrLrOd3ol3SB3red6i+SR1ogc8FW5brfLU174FIGPJVWsdUIomc/1xr1lT31zuyXwG7Od\n97q9W6YXblW4WinWDz+ksQlH3owLwOImUPri1lLRAiffgxMQavxyBIm/EHKksKLCvYKmJW3aLUAm\nfAHTKZeApgoTKPzCN50az1/8Va/GNzjh3DR9UqVBR+RK3LVoMv9SW1tE7RhV3WfkIpglbT5GbR0o\nid0+jxtnu+IpffGa8tZ9nvq2+zwC764bJfBz0/TWGZBVuP4gW4/8jNb1zcew/FT28Oxfz3dZ5uer\npA6UBc9Nt++eBP4FpRqPImcNQeIvgFTMZBLFSaR4qhsoSVs/KZ0AssMlBK5ueLlNCfvcLSeFvsJh\noMbVVmnQ+ePlA5XETngWCsmc95bMG9kH+CrcU94eH11irXjHe1aKbrdE3pptNfLmsWMEvvK3k8Bp\nozCRUlPhaqXo8AoECX2OQu+z4vY90uuhNds9S6VX46nnZqjxeQgSfxFwdEJNpMgI1FrQVCsFGPrh\nD/zEDh0ZNzgX46Bk2+RUWCqn/kt+wBlN740DBxyx7PzybKsAQHPTKXLFDUrFbe9tjNDuh7SBaQfT\nBhgStCXh5yps6vbW7KvZKZeQd1vZd1MS+G67GdgoYyrcWileUfPQbbtEofew4kOva1oq7AdwQlfg\n9NR4dACqIUj8mZEjhVaFdwTOgqZaKUBJ5FyvoZJMYaLkjKYoQmlCZYt9YakwkUJy5/0Sx0EhbAlg\n12yAFUoiVwIGfB+coFWiZK6q+yBt9H2YInXiuUlcH3vKWiH5crlW7FyZ9h6BGw/c+uBHLHsCz59p\nXYVTeRPH7jwWeTrlzq75+xLHu4niqL5Wu26HllwAuAfS98Oq8c9R5HQQJP7sKIfazCrcFDT176YW\nK/n3eQynBse7FZr2nHK7n86TkcG0ni0Vkjoc+yTdJ3W2xc4l8vO2wXa371/SAsi+rbVJPLVtyRwY\nFiphjvEslMqkEuX75WyrwTveU9+63Spy2xXeK3bOJW+5rynwNHtmJuwdtj2pq72iKlwJXsm7liF/\nEJgPp1/O+OF3VAqcFDxZjaflgCJI/BlRqvANShUuBU0gWylA+lTUD38kaKWgU8sk7GWx79B/eYf2\nSbJVlOAtkR8AYAs05zNWhyOaFmitH04S8sgcuNxGuUEev4XH3KAk8ockU4g5CRVtM0ba2nbMVhkj\nb2Bgn7CIWSNwEjVVuWevqAqnl04//NAr9nbWxBKzoL44XzOviaLA+RN59h+q8U2ocYMg8WeFziNI\nZjYqXH1u7Vb/BCAR0wMnYW+x67/wS+S8cLJN0G0f2ifAAfvuV0eJ/IwGG+z7YudhtUz2yg3QnMRe\n0dqrJk64X0keGBJ6erCEMRvFXtWe8p5KptRgz71y9l1qq1iLZUSZ1wh8j23hgVsFTq/b2iu8DqwK\n139p1g+nR/5o8PXams+oGud9gAgSfyb4KpxELheh5mctWnPYTGWeFNgS224QIfXF9Quc7ZMTVjj2\ntgm3A0P7hER+RiP70rYGZ5zRomlOOG+bQpUXZO4pb1vY9JaBOqlDzrGCr8CfcixafQ7A5bZKzRuv\nKPMp9a0EvusIvUbgJblnm8WqcABVYue+0+mJVAe/A+yqDyCTNr8IocY9BIk/AxaLf5ZMOBMpXO5k\nh1opQPa+mQ9/wCdzPqUv1bITSepnUoUfgN4S0S8ntydC3mKDHTwiT6o772u6YmiDM5ZdkmU2mavy\nVpIes07UNpmyUtITHuKS4Wi9H86xQmfNUhmzWGp+OC4jb9onSuB7bKvqXPfxnBo/VcWtVsrANz81\nwBSZ0//WQmaDsnOQtv2OrsDJDSeUavwXLBb/fH9//98+PJEHiT8LOEW9JlK43kGtlLne9x2SDTHY\nvgBaSRb8fYnmU/LAaanQNgFKzzvbJ9oN+yBfUkvkh352lxMaME+e7ZXErOkLv6ySOWAIXUnaEnbN\nOtGr1yN3OPvZ5iHwvi2enTLHYhmxVHQEwjHyBjCqvsfUuRL4Dhs5L331bK8oaSuxD5IpdyN8Ojay\nLAuenCzCtVRaWeb6F+dkHw9B4k+M3LHHU+HGSnkK7NFfy31C5XbX/z1e9d52tk287Z7qTvlxX3Vb\nKyUr8KzKW5yrNguAQp0DHaEDQ2XuqW4daOts9pGka/2fHptOAXyiBsZz4546N4obQDEOeE15A6iq\nb2ufTKlzJW4+Bq8ftVnsgGhEQeR7XPZPZxJU3mT5cj1slSDxZ4COWKUVG1mnlaJ+uNMLs4CQNYBc\nELxBp2AaoD33lgraUo2zmEnbRG0Va5/kTHnat+z88i32OGCFNLPLtrBPqMotmQMY2Cxo0mN4hA50\npH4DLJS4LZFbxa3ErN448dgBsNIbM4RnpehjVwicIw1OETcAeLYJAJe8PXtlj01vk2hbVeD0wdVG\nYSLF2ivMh59OTbJSgHT9eZ+BYurHU/+ZniCWSu27RIH0sREk/oQYqnDxwOdgarwRgrUeQNIa2VI5\n3q0GaryfFxFZRXGKLUvky25fScRpm0fUNTJPL6lsd8Cy/4NuCR1AQepYZ5VeJXZgSOZql8wZjWAs\n2eLB9m3xyFu2e4QNlKSdnqpP3Gm9VNVsM0beSvAka9oqZSRxVRA4yd2qdKvCewV+arKVoj+WcyPd\nTKhMflYk8FDjiiDxJ4VV4RUr5ZLC5U+Uao8FQfqHQFbpXaef5fpQqPE09lC2TwD0oxVaIi9VN/of\nALu03dwAAA8CSURBVJIz90+ROZW6ttthixUOA0JPLyGdm6S+26bZ1YE6sQOZ3AEheH2vPP/7knHC\nap0RzeenM8uf2+GyR9gABqSdnnbZo/Igy0rc3DaHvMs2y95SUZtmh22hvFkELe2VUoUf71a5qEnf\nW69NXQbmE3vhi1tLhT05Q40DQeJPhroKt9aK4NJ3P7sc6XS0VPj3s1Pj9MaxPuC8bPpDckedIZHX\nVXcLRhNVaY+ReTpHqc7Ty1WFnggd3fOypN7g3P+oeMQOZHIHMsH3+5y/7v22BxQ2z85nZbcV81nK\njDoeYad17RGb1XY6Zqi4tY1H9HPIW7fX1LklcCX3/d83OHXk3RP43SJdfxw2gVaK/liqQp9jq7jQ\n71OocSJI/AmQIoWaSBlRCPYyG0umUHUDZSeXW2RLxVHjALJKArBdJolEst52Q8oCiRyt6s7KWol6\nmsz32GCJo0Po2UJRhc726XmUpJ6OK4k9vUzuO/YWEwmeUKJXKOlfCiVlhe30otE7S9ZASdi6z5I2\n23qKm+tK3Gx/CXnrOVgItUVUq85J4H2skF64qnAgK271xuf2kv0u6/1YKp9RRlyG37WPGjkMEn8S\nqDqwCRQnkTLn358tEnFm8C8oLRVPja+zrWKJnOLIqu5kt5REfQmZp3OV6pyPo8vWFyepA+h7gyqx\nAxiQe2qTZZ6SfDrPkKyV9B8Kr8v52Zy0HERqNWh3MOStBKzHeH44jx8jbm0/l7xz4TK3d9X5cdMT\n+PH7NqvwPUoVzt65KjAImzIC6kVn5sWLDWqpWJvlY9LZx3zVTwhfhQPluJvALOY+m3sgWyfqMdqJ\nDrSDzC26IlOiQiXytj1j88lX3STnS8mcPwAeoTc4DRT6boaNosQODMk9tcvjV69QmtyD+R87LHGJ\nGV6iNo2ZJXE7/Vk+viRrPWfNC+f5PFtFJ3ywqlv3X0renjqnB+4SuNoo9L9VhXO/3Q5cmBgiYQP5\nu8YvxcdW40Hij4aqcN3G9QvSKRacOIHXLtU4O7DZdhzGdYNM5KcGy9tdIvL1Aft2g+WnOlHPJfPU\nWeg0SujpnSjbAxjYKMBQbS87MlKVXVPie/MDqWRf4kuV4MdgiVph1fmUEtc2HmHruiVtu22KuNn+\nIeTNNkrge5K3Evh3ZLK+Q7pWqcLVZiFxT/nhk5aL2iqtrPNL8fEo7eO94idEqcIB30oBRolcPe41\nsmXC6xIYqnEeRz/9s6z/QFI4VOTrTjl1GfKmPeO0btC25weROZALpDVCT20kXdIpbiCTcibxTOzl\n9pLceR5iacjTWipzyHpMmc+ZQHjMSrHn0LZTlootgPoKvYwjjhG3tn8IeRcplDECVxvlDtlmoQpX\nK4U8zPtZqRUlcDsoFotEH0+NB4k/CqrCrZXyAHB8bI+seZGT6Dey7zs60kY5AUMD4LQA2jaROcRe\nWR+AdtkReTNJ5gBw7FR2jdCBbKMAZRGy6aOMmdjTvnFy/4bbgqCzyh6qamur6OP4eFi3bU+ZHwzp\ne1YKUBK9R9ZpfUjYbO8RvhZQx4hb24+RN4Bx9a0euBL4d7N+lnWYZe4D5qv0HhukUbKobPjd22No\nY75/fKxX+4RIkULmwoEsm1tzPwIdxU9TKDU1bmGJnMVOJfSe8LO9YlU5gKoy12x3jdDVRhl63tZK\nsWTtk3tqUxK8HkfojEWW8C3qNss0xsbRtgrczojjkbhf9PSjiJ5C9ywYVedjxG3bW+UNoK6+lcBJ\n2mqhkMD3sj6mwvlx1XrfVmG/Z1TkHOXw46jxIPEHINko2lee917x0tmmXeg1+w1ki0Qr+2ODZCmR\nE5wphT8SZwxVeUfm++9bLNcH12YBUCV0ZrhpuRyx6gl4Z6wUwPO8rZ0yJOpG3rsiH26IumaLjBH6\nY2GJm7A2TGmjlBMRp/ZDgq/lya0NM0XaebtP3AB82wSYr749AldFbverxQJpk9+EjLFBswroVG5Z\njX8UIg8SfxC0cOn53jOSKJ5lssJQQbOtJWqiJ2iUxL1ytlkyRwus7wfqfIdE7CR0EnRNdZe2SVbq\nAApiB4aK2/O8yxQKbZUhKau1Ytt7qOXH52BsEgSrvIHxgudUltwjaz2u5otznyVtHusRN1BR3elJ\nX0beStAegavFoioc0kbXXWhKRZd1/7fufvavwFUjSPxClCrcWikXgsTM06igVCK3HjhMOzhtpsgc\nGKrzuyXQnntCB4CmPWO57myNEVJPpzv3pDy0TUpyZ/tv+NLvq6ntqQhhzSaxhP7NbfUweMQNDMl7\nTgTRU+iegh9T6B7hW9IGkDvqdMs9cQND1Q2Mk7fun0Pg3zH8YbBWiyXzexjUiJkiSgcW2n8INR4k\nfjFsT0xbxLTr+3LbPUrfG6grbRI5UFfc/Be5dtoAPpl/R/7x6Al90Vn8SaEf1x2hA9i3257ILakD\nwPFTJvGa4k7bDviOuuL+1rcbEj3bWHjq+tIYof5I2ALlFPwip9cpqJXlcYsFmLZSeKwlfS1MAhXS\nBuqKGyiJm/stKdvjxshb23gEbokbGCr12bAJlvevxhf394OfukAFWYXb8VGcSZD7dnY/UlfiRprp\nkLQ60wvJWdt6x8G0ZztMtOX62izbYwFg3V0nnVIHUBC7rpPcl59KIlcSrpE8MaW25xK6xVPnxIkp\n4gbGVfqlCt0jawAlYQNDpQ2UxUUSKrd7xK377LEeeQOXE/jJeYwzRInvkQuXvD8B+EMOsPv3fbv3\nrMZDic/EYvFPXSb8wuFlPagatzbJlCLnNFequNU+UWUO0xYox1vhj4Ums1Sl/5DzUKmj7dQ6kloH\nCsUOZDIH6gT/HZnkgXHFbcnZI2Kq/BoeQt4WU2TuqfihnVIfa6WIJf5dCp4VotZthcIGSsIGhkqb\n+yxppyc5JGM9hy1KziFv4BEE/lDwe/qURtrbQ5D4bMz1vbXoolaKbsd8Ij+gJGbtnFYjaJKxqvMD\nskJn3UdJm8crqXMbkDsSASWxs1FP8F1CQ5U7gD2+9MuAT/R2O5BIX8lZiV8xRdJjBc+5qPngeX8l\nsfJ3U+Q081HqzDguQQNDkgZ8os5PpmuD3MbmsWukzTZWcXN/TaHzOZyd9bHjJgm8psLnWiXvO6kS\nJD4DSYU/dMxiQ96KGpHzwrbRQpI5yZunLQqVGPrhJHRgnNSBkrR1XY/32ulkt2sAWHTPsZWX3/bH\n9lG8tXxb23Mie1m3sCQPlD8CY/COnYuCUEdwPg2J3D3WtvPIGahH7ixJa9sxsh5rp2rdiwBaAuay\nPa6mvPVYXuPAMyhwH4vFP93f3//3d0fkQeKzYD3uS942Ty0IoVsiZ8ab5Ggz36rMgWlCB4Y2iSV1\nIH/JWjlen2qN3OGcy+5fA/i/le06Nm9j3teNE91rK2S6nv7m95/ETNIHMCTbMdQmCvb+BNjLwj6l\nWtyutt0Sre73yFof07abIm2eY4y4eZ6aYq8peDeNYlX4XOj1o5nd94UobE4gq3CPxL3CJffp+sZs\nswNm8cGQC5dASeYsfAJ+kRLS1p6D2+Bsn1LcHmnbc9n9U/NO6vm9dt55xtqNtX8pjA3c5JH41ByU\nNv48NZ+o3W6JWvePkbXunyJ9z2f3iJvLYz57NUqoBM7tJ7OtVtjkvtKCeW9qPJT4JJSYp8Bf+9as\n816/mZTTRpUzT0uiVHtk7WyzapukXhQl5aUcUBL7N9kHDB9Lt1kS/r8jbfVx9fxeO3uOWhtijrPx\nyNrzLMwZtKnm4MwhdtvGs1NsOz2HPj9vHO8psh7bV7NZ7HmnvHZXQ1rS9ghcn+B8b/w9qvFQ4iMo\nVThQV+K6rGyo+2rHAKPKHJinoueQr6fYCXvu2n7bZo7q1sf2jvOO956DxUO+jw+RLQ+piY7xiufm\nzCH1p1Ln9jnMKYzOVe1z1fws8uYDq6rW7SdzXO2Y4TnekxoPEq+gLGZaS8Tb5i1ba0W31Y7VffZJ\ndfc1UgeG9sfYPsD3s5Vg5xC5bWfbeu3t44y1A+apb+BlFLjF3Ml/H6PKvXZjvvpcFe8p9zFfXbeN\n7ZtF2tqwprztSX+ObPNUO889PNd7IfKwU2bDK4zoNs8y0fZajeT4Di3KC8/GEQ2p84vA6/AbhsSu\nSROgVOh3ZptdrqlqIJHo38zT8o4jpki99ji14y2e4l/x2NX/+ETi9L/8Wn219thTxVBgXL3bH5E5\n6n3MftHnpM9lVBd6pO2t218Dj7y94+zy+0cocQdJhQPj6nlse229NdvHHsM7vnZsBaozxkh1zJ/2\niO6hfvaYmp5S0JcWLS8IlTwYl/YfGit+AuOK3lPxD/HV7fpY+7EfjUna0IN/Otu9bZ6vM5e8az8C\n4z8O70GNhxKvYo7UY4VRVbidlds7r1YjVaED5beIH88f8Incbtd9AO5ln/3yep3Y7OXsEaH3ttSu\nohrxzrnqLrkyXzOVYjFF1Io5ar/WpvY4c+KMwPAH6GItVyPpsX22IDnW3tt/6Y/Bx0AocYOhCgdK\nRqmp4jElbbeP7bPsNXYO23bsMWptLmn3QExpnYeo5msIGVz6r36Osn/Sr6v3BL2/A7bdnDZjkt4e\nP0XoU+fwyHvqxyNvu3Y1Hkq8CvW7NdPHD1/9b50A89tEW7YfUzKWVNV3J7yP7hIyrn30Y+w45XlU\njh0QjznPLOFkzn2Vw2E8xKudWzWdOvfYeWrHzjXna+ef87dgTMXP2W8fY4z436dCDyUuyCqcmOsd\nzFW3c4+tHf+YY8aOI+b8pl8a/3isZH6POuOxZHLpj8GcH4Kp5/TQH4mx4+YoesJ7fnOPn/4xuWY1\nHiTeYUjgiseS39Txc4jxEjK7lDgfk8t7Dl/jNXKCr4W5SvsSPCadcenzueSx5vx4TT3+1OM9/Mfo\nWon8PcqcZ8CcC3WMzC797//cH8tbNpQ/EoETz0HkT4Xnjus9pcXxsaKFRCjxQCAQuGJ8eu0nEAgE\nAoGHI0g8EAgErhhB4oFAIHDFCBIPBAKBK0aQeCAQCFwxgsQDgUDgihEkHggEAleMIPFAIBC4YgSJ\nBwKBwBUjSDwQCASuGEHigUAgcMUIEg8EAoErRpB4IBAIXDGCxAOBQOCKESQeCAQCV4wg8UAgELhi\nBIkHAoHAFSNIPBAIBK4YQeKBQCBwxQgSDwQCgSvG/w/LZmo3q0TkEwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fb69691fe50>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import scipy.stats as stats\n",
    "%matplotlib inline\n",
    "\n",
    "delta = 5e-3 # workaround for scipy.stats.dirichlet does not accept input equal to zero\n",
    "corners = np.array([[0+delta*0.75**0.5, 0+delta*0.5], [1-delta*0.75**0.5, 0+delta*0.5], [0.5, 0.75**0.5-delta]])\n",
    "triangle = tri.Triangulation(corners[:, 0], corners[:, 1])\n",
    "\n",
    "draw_pdf_contours(stats.dirichlet,[2,2,2])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Normal (Gaussian) Distribution\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "Norm(x|\\mu,\\sigma^2) =& \\frac{1}{\\sqrt{2 \\pi \\sigma^2}} exp\\left(- \\frac{(x-\\mu)^2}{2 \\sigma^2} \\right) \\\\\n",
    "=&  \n",
    "\\sqrt{\\frac{\\lambda}{2 \\pi}} exp\\left( - \\frac{\\lambda (x-\\mu)^2}{2} \\right) \n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "$$ E[X|\\mu, \\sigma^2] = \\mu $$ \n",
    "$$ E[X|\\mu, \\sigma^2] = \\sigma^2, $$\n",
    "\n",
    "where $\\lambda = \\frac{1}{\\sigma^2}$ is known as the *precision* parameter that simplifies computations in many cases"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Consider $\\sigma^2$ fixed and $\\mu$ as the model parameter, then the posterior probability is given by\n",
    "\n",
    "$$\n",
    "\\begin{align*}\n",
    "& p(\\mu| x; \\sigma^2) \\propto p(\\mu,x; \\sigma^2) \\\\\n",
    "= & p(\\mu) Norm(x|\\mu; \\sigma^2) \\\\ \n",
    "\\propto & p(\\mu) exp\\left(- \\frac{(x-\\mu)^2}{2 \\sigma^2} \\right) \n",
    "\\end{align*}\n",
    "$$\n",
    "\n",
    "It is apparent that the posterior will keep the same form if $p(\\mu)$ is also normal. Therefore, normal distribution is the conjugate prior of itself for fixed variance"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Given prior $p(\\mu) = Norm(\\mu| \\mu_0, \\sigma_0^2) $ and likelihood  $Norm(x|\\mu ; \\sigma^2)$. Let's find the posterior probability,\n",
    "$$ \n",
    "\\begin{align}\n",
    "&p(\\mu|x; \\sigma^2, \\mu_0, \\sigma_0^2) \\\\=& Const \\cdot Norm(\\mu| \\mu_0, \\sigma_0^2) Norm(x|\\mu ; \\sigma^2) \\\\\n",
    "=& Const2 \\cdot exp\\left( -\\frac{(x-\\mu)^2}{2 \\sigma^2} - \\frac{(\\mu - \\mu_0)^2}{2 \\sigma_0^2} \\right) \\\\\n",
    "=& Norm\\left(\\mu; \\tilde{\\mu}, \\tilde{\\sigma}^2 \\right),\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "where <font color = \"red\">\n",
    " $\\tilde{\\mu}=\\frac{\\sigma_0^2 x + \\mu_0 \\sigma^2}{\\sigma_0^2 + \\sigma^2}$ and $ \\tilde{\\sigma}^2=\\frac{\\sigma_0^2 \\sigma^2}{\\sigma_0^2 + \\sigma^2}$. </font>\n",
    " \n",
    " Alternatively, $\\tilde{\\lambda} = \\lambda_0 + \\lambda$ and \n",
    " $\\tilde{\\mu} = \\frac{\\lambda}{\\tilde{\\lambda}} x + \\frac{\\lambda_0}{\\tilde{\\lambda}} \\mu_0$\n",
    "\n",
    "<!--\n",
    "Clean it up, \n",
    "\n",
    " $$\\mu \\sim Norm\\left(\\frac{\\sigma_0^2 x + \\mu_0 \\sigma^2}{\\sigma_0^2 + \\sigma^2}, \\frac{\\sigma_0^2 \\sigma^2}\n",
    " {\\sigma_0^2 + \\sigma^2}\\right) $$\n",
    "-->"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Consider $\\mu$ fixed and $\\lambda$ as the model parameter \n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "p(x|\\lambda; \\mu) \\propto & p(x,\\lambda; \\mu) = p(\\lambda) Norm(x|\\lambda;\\mu) \\\\ \\propto & p(\\lambda) \\sqrt{\\lambda} \\exp \\left( -\\frac{\\lambda (x-\\mu)^2}{2} \\right) \n",
    "\\end{align}$$\n",
    "\n",
    "More generally, when we have $N$ *observations* from the same source, \n",
    "$$ \n",
    "\\begin{align}\n",
    "p(x_1,\\cdots,x_N,\\lambda; \\mu) = & p(\\lambda) \\prod_{i=1}^N Norm(x_i|\\lambda;\\mu) \\\\ \\propto & p(\\lambda) \\lambda^{\\frac{N}{2}} \\exp \\left(- \\lambda  \\sum_{i=1}^N \\frac{(x_i-\\mu)^2}{2} \\right)\n",
    "\\end{align}$$\n",
    "\n",
    "From inspection, the conjugate prior should have a form $\\lambda^a \\exp (-b \\lambda  )$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Gamma Distribution\n",
    "\n",
    "$$ Gamma(\\lambda|a,b) = \\frac{1}{\\Gamma(a)} b^a \\lambda^{a-1} exp(- b \\lambda) $$\n",
    "\n",
    "$$ E[\\lambda] = \\frac{a}{b} $$\n",
    "\n",
    "$$ Var[\\lambda] = \\frac{a}{b^2}, $$\n",
    "\n",
    "where $a,b$ > 0 and $\\lambda \\ge 0$\n",
    "\n",
    "N.B. when $a=1$, Gamma reduces to the *exponential* distribution. When $a$ is integer, it reduces to *Erlang* distribution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYYAAAEKCAYAAAAW8vJGAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl8VNX5+PHPk30PYUnCDoLKqgiyiUuUVhEVra0LihXU\nr7utVutWFLRq3ZcquFWtgkprpaL+EBUxKmURJCKC7BASQkIC2dfJzPn9cSdhkkySmWQmyYTn/Xrd\n19y598y5Z4Ywz5xz7jlHjDEopZRSNYLauwBKKaU6Fg0MSiml6tDAoJRSqg4NDEopperQwKCUUqoO\nDQxKKaXq0MCglI+ISLGIDPBRXveJyGvO/f4i4hARn/x/FZG+IlIkIuKL/FTno4FB+Y2IXC4ia0Sk\nRESyRWS1iNzU3uXyloicISJ255dpkYjsE5F/icjJrumMMbHGmL0e5JXR3DWNMX8zxlzveqhlpQcR\n2SMiZ7nknWGMiTM6iEk1QgOD8gsRuRN4DngCSDLGJAM3AqeISGi7Fq5l9ju/TOOACcBW4DsROdPL\nfIRmvuRFJLiFZVTKJzQwKJ8TkTjgIeAmY8x/jTGlAMaYjcaYq4wxNme6qSKyQUQKRSRdROa45FHT\nfDLT+Qv9kIjcICIni8hGETksIi+6pL9aRFaKyLMiki8iO0VkovP4PmeN5fcu6Ru9dnOMMVnGmDnA\nP7ACX02eDhE5xiX/zc4aRoaI/ElEooClQC9ns1ORiCSLyBwR+UBEFohIAXC189gC148VuFZE9ju3\nP7lc9y0RedjleW2tRETeAfoBnzivd1f9pikR6SkiS5yf8XYRuc4lrznO2tHbztdvEpHRnn5WKjBp\nYFD+MBEIAz5uJl0JcJUxJh44D7hRRKbVSzMOGAxcBjwP3A+cBYwALhWR0+ql/RHoCrwPLAJOBgYB\nVwEvOb+cPb12cxYDo0Uk0vnctSbwD+D/nDWMEcAKY0wZcC6Q5Wx2ijPGZDvTTwP+bYzpArznJj+A\nFOd7OQe417V5yA0DYIz5PbAPON95vafd5L3ImSYZuAR4rF5N6AJnmeKBT4B5TVxXdQIaGJQ/dAfy\njDGOmgMi8j/nL/kyETkVwBjzrTFms3P/Z6wvqDNc8jHAw8aYKmPMcqAUeN8Yc8gYkwV8B5zkkn6P\nMeYdZ9v5v4A+wEPGGJsx5kugCivIeHJtT2Rh/ZLvUvM2Xc5VAcNFJNYYU2iM+bGZvFYbYz5xlqei\nkTRzjTEVzvK+BUz3oqxuO5pFpC9wCnCP83PaiBXUrnJJttIY87nzc10AnODFdVUA0sCg/OEQ0N31\nLhpjzCRjTAKQh/PvTkTGi8gKETnobEK5ASuouDrosl8O5NR7HuPyvP45jDF57tJ7eO3m9MYKXgVu\nzv0WqyaSLiJfi8iEZvJqrkPaAJkuz9OBXp4WtAk9gcPO2oxr3r1dnme77JcBEb66Q0p1TPqPq/xh\nNVAJXOjmnOsv13eBj4DeziaUV2nkl60f+OLaFwMbjDHl9U8YY34wxlwE9ACWAP+uOdVIXp7cIdTX\nZb8fVo0FrJpUlMu5nl7knQV0FZHoennv96A8qpPSwKB8zhhTCDwMzBeR34pItFhGUfcLLAbIN8bY\nRGQccEW9rFobJJp6fXPXbjQvEenl7Ky+BrivQUKRUBG5QkTijDF2oBiodp7OAbo5O+i9IcADIhIp\nIsOBWVjNX2D1q0wVkQQRSQb+WO+12cAx7t6PMSYTWAX8TUTCReQE4FpgYTNlUZ2YBgblF8aYp4A/\nAXdjfRlmAy87n69yJrsZ+KuIFAKzsfoF6mTj5fMGxWjieXPXrq+n866cYuB7YDhwhjHmq0byvwrY\n42ymuh6YAWCM2YbVMb7beWdVcjPXdc37G2An8CXwpMu1FwA/AXuBZRwJGDUexwoqh13uZnIt63Rg\nIFbt4UPgAWPMimbKojox8ecYFxF5AzgfyDHGNOiwEpErgHuw/tBKsG5v3OS3AimllGqWv2sMb2Hd\nWteY3cDpxphRwCPA634uj1JKqWaE+DNzY8xKEenfxPk1Lk/XUPdOCKWUUu2gI/UxXAd81t6FUEqp\no51fawyeco6ynAWc2t5lUUqpo127Bwbn7XGvAVOMMflNpNM7IZRSqgWMMV7dYtwWTUlC48Px+2Hd\nHneVMWZXcxkZY3Tz0TZnzpx2L0Nn2vTz1M+yo24t4dcag4i8hzXxVzcR2QfMwZpczRhjXgMewJrw\nbL6ICGAzxozzZ5mUUko1zd93JTU5mtQY83/A//mzDEoppbzTke5KUm0oJSWlvYvQqejn6Tv6WbY/\nv4589iURMYFSVqWU6ihEBONl53O735WklGpbAwYMID09vb2LoXysf//+7N271yd5aY1BqaOM8xdk\nexdD+Vhj/64tqTFoH4NSSqk6NDAopZSqQwODUkqpOjQwKKWUqkMDg1JKqTo0MCilOpV58+YxduxY\nIiIiuOaaa5pMm5+fz29+8xtiYmIYOHAg77//vt/K5c21UlJSiIyMJC4ujtjYWIYOHeq3crmjgUEp\n1an07t2bBx54gGuvvbbZtDfffDMRERHk5uaycOFCbrrpJn755Re3aR0OB3/84x+ZPHlyi8rlzbVE\nhPnz51NUVERxcXGj6fxFA4NSqlO56KKLmDZtGl27dm0yXVlZGYsXL+aRRx4hMjKSSZMmMW3aNBYs\nWOA2fVBQEOPHj+eUU07xukzeXgto17EmGhiUUgHhggsuICEhga5duzZ4nDZtmtf5bd++nZCQEAYN\nGlR77MQTT2Tz5s2Nvuarr75qUGPwpFwtudZ9991HYmIip512Gt98843X7681dEoMpVQd4tUY2ca1\n9Afvzp07+ec//8mpp57KBx98wLnnnsvvfvc7PvnkE98UzKmkpIT4+Pg6x+Lj4ykuLm70NStWrODM\nM8/kvffe4+DBg9x+++0elcvbaz355JMMGzaMsLAw3n//fS644AI2btzIwIEDPXhnrac1BqVUHcb4\nZmuJsrIyLr74Yu68806mTJlCdnY2w4cP9+0bdIqJiaGoqKjOsaKiImJjY92m37VrFwMHDmTGjBlc\nccUVPPPMM3671tixY4mOjiY0NJTf//73TJo0iaVLl3p8vdbSwKCU6jAWL17MiBEjSEhIwG63s3fv\n3to7cqZOnUpsbCxxcXENtvPOO8/rax133HFUV1eza9eRxSM3btzYaCD69ttvmTp1KgDbtm2rrQF4\nUi5vr1VfW89vpU1JSqkOIzc3l9GjRwOQmprKuHHjWL58OZMnT/b4F7Pdbsdms2G326murqayspKQ\nkBCCg4PrpIuKiuLiiy/mwQcf5PXXXyctLY2PP/6YVatW1aaZNWsWIsKbb75Jfn4+I0eOBGDhwoXc\nddddAB6Vy5Nr1SgsLGTt2rWcccYZhISEsGjRIr777jteeOEFj96/L2iNQSnVYUyfPp39+/ezbNky\n9uzZQ2xsLHl5eYgXHR+PPPIIUVFRPPHEE7z77rtERUXx6KOP1p6fOnUqjz/+OGCNeSgrKyMxMZEr\nr7ySV155pc6YgYyMDE499VQALrvsMtauXcvbb79NcnIyM2fO9Oq9NXetmnLZbDZmz55NYmIiPXr0\nYN68eSxZsoRjjz3Wq+u1hk67rdRRRqfd9ozNZmPUqFH89NNPDWobHZEvp93WwKDUUUYDQ+ek6zEo\npZTyGw0MSiml6tDAoJRSqg4NDEopperQwKCUUqoODQxKKaXq0MCglFKqjoAMDM9nZFBut7d3MZRS\nqlPya2AQkTdEJEdEfmoizd9FZIeI/CgiozzJ94mMDAqqq31XUKWUUrX8XWN4CzinsZMici4wyBhz\nLHAD8IonmUYHBVGqNQallBtXXXUVvXr1Ij4+niFDhvDGG280mlbXfHbPr4HBGLMSyG8iyYXAO860\na4F4EUlqLt/o4GBKHQ7fFFIp1ancf//9pKenU1hYyMcff8zs2bNJS0tzm1bXfHavvfsYegMZLs/3\nO481KTo4WGsMSim3hg4dSmhoKGCtmywiddZBqKFrPjeuvQODu4mdmv00ooKCKNPAoNRRxZs1n2+5\n5Raio6MZOnQovXr1ql1gx5Wu+dy49l6oJxPo6/K8D5DVWOK5c+cCkJGTw5opU/jVhRf6tXBKHY3k\nId8s+mzmtOwXry/WfJ43bx4vvfQSq1evJjU1lfDw8AZpOuuaz6mpqaSmpjabrknGGL9uwABgUyPn\npgL/z7k/AVjTRD6mxvTNm8272dlGKeU91/9LHU1paakZOXKkOXz4sDHGmKlTp5otW7a0Ks8bb7zR\nvPjiiw2Op6Wlmejo6DrHnnnmGTNt2jS3+ezcudOceeaZtc/79OnjcRm8vVZ9U6ZMMS+99FKTaRr7\nd3Ue9+p729+3q74HrAKOE5F9IjJLRG4Qkeud3/RLgT0ishN4FbjZk3yj9K4kpTolf6z5XH+t5Rq6\n5nMTvI0k7bXhEg3/sH27eT4jo8noqZRyjw5cY3j22WfNU089ZYwxZvny5WbmzJnmyy+/NA6Hw6PX\nHzx40CxatMiUlJQYu91uli1bZmJiYswnn3ziNv306dPNFVdcYUpLS83KlStNly5d6tRQZs6caWbN\nmmWMsX7hL1u2zBhjzOzZs81bb73l1Xtr7lo1CgoKzOeff24qKipMdXW1WbhwoYmJiTHbt29vMv/G\n/l3paDUGf9G7kpTqnFq75rOI8PLLL9O3b1+6du3K3XffzQsvvMD5559fm0bXfG5eQC7t+Wh6OmV2\nO48ec0w7l0qpwKNLe3rmaF7zOSBrDNrHoJTyt9DQUDZv3hwQQcHXAjIwRAcHU6Yjn5VSyi8CNjBo\njUEppfwjMAODNiUppZTfBGRgiNJJ9JRSym8CMjBEBwfrXElKKeUngRkYtClJKaX8JjADgzYlKaWU\n3wRkYIjSu5KUUspvAjIwROt6DEop5TeBGRi0KUkp1YitW7cyefJkunTpwnHHHcdHH33UaNq2WvO5\nqqqK6667jgEDBhAfH8+YMWNYtmxZu5erMQEZGEKDghCgSoODUsqF3W7nwgsvZNq0aeTn5/Pqq68y\nY8YMdu7c6TZ9W635XF1dTb9+/fjuu+8oLCzk4Ycf5tJLL2Xfvn2tLpdfeDsda3tt1JtSNv7bb83h\nqqomp6FVSjVU//9SZ/Lzzz+b2NjYOsfOPvts8+CDDzZIW1paasLCwszOnTtrj1111VXmvvvuazT/\nd99918yePdsnZT3hhBPM4sWLfVIuY3TabUDnS1LqaOPJ2srGzeyixhh+/vnnBsfbcs3n+nJyctix\nY4fbhXpaUi5fa+81n1tM50tSyk88XPugWS2c2rs1az4PGTKExMREnn76aW6//XZWrFjBN998w1ln\nndUgbVuu+eyqurqaGTNmMHPmTI477jiflMvXNDAopepqx7UaysrKuPjii/nmm29ISEjgxRdf9Hj5\nS4CQkBA++ugjbr31Vp544glOPvlkLrvsMsLDwxukjYmJoaioqM6xoqIiYmNj3ea9a9cuBg4cyIwZ\nMwDo27cvt99+uxfvzqq9zJgxg/DwcF588UW3abwtlz8EbFOSrsmgVOfjizWfR4wYQWpqKrm5uXz2\n2Wfs2rWLcePGNbhWW675XOPaa68lLy+PxYsXN7rOQ2vXh/aFgK4xaB+DUp1Lbm4uo0ePBiA1NZVx\n48axfPlyJk+ezNKlSz3KY9OmTRx33HHY7Xbmz59Pdna222U4o6KiuPjii3nwwQd5/fXXSUtL4+OP\nP2bVqlW1aWbNmoWI8Oabb5Kfn8/IkSMBWLhwIXfddReAx+W68cYb2bp1K8uXLycsLKzRdJ6Uy98C\ntsagTUlKdT6tXfMZYMGCBfTs2ZPk5GS+/vprvvzyS0JDQ2vPt8eaz/v27eO1117jxx9/JCkpqbaG\n4To+wZty+VtArvkMMGPLFqZ07cqM5OR2LJVSgUfXfPbM0bzmc8A2JemaDEopf6pZ8/loFNBNSTpf\nklJK+V7gBga9K0kppfwicAODNiUppZRfBGxg0DUZlFLKPwI2MOiaDEop5R+BGxi0KUkppfzC74FB\nRKaIyFYR2S4i97g531dEVojIBhH5UUTO9SRfHeCmlFL+4dfAICJBwEvAOcBwYLqIDKmXbDbwL2PM\naGA6MN+TvHWuJKWU8g9/1xjGATuMMenGGBuwCLiwXhoHEOfc7wLs9yRjnStJKaX8w9+BoTeQ4fI8\n03nM1UPAVSKSAXwK3OZJxtqUpJRyx5v1klNSUoiMjCQuLo7Y2Fi/zkc0b948xo4dS0REBNdcc02T\nadt7zWd/T4nhbn6O+pN5TAfeMsY8JyITgIVYzU4NzJ07t3b/2IkTKa23mIVSSrmul7xhwwbOO+88\nRo0a5fZLX0SYP38+s2bNajZfh8PBHXfcwc8//8xXX33ldbl69+7NAw88wOeff055ebnP3kN9qamp\npKamel0+V/4ODJlAP5fnfYCsemmuxeqDwBizRkQiRKS7MSavfmaugeFAZSV3rl/v8wIrpQJXWVkZ\nixcvZsuWLURGRjJp0iSmTZvGggULeOyxx9y+xtMJBYOCghg/fjxxcXHNJ3bjoosuAmDdunXs3994\ni3lL3oOrlJQUUlJSap8/9NBDXpfV301J64DBItJfRMKAy4GP66VJB34FICJDgXB3QaE+7WNQ6uji\nydrKLVkv+b777iMxMZHTTjuNb775psky+GrN56Z0+jWfjTF2EbkV+AIrCL1hjPlFRB4C1hljPgXu\nAl4XkTuwOqKv9iTvmkn0jDFezdWulGqatLIZooZx+dXqjdas+ezteslPPvkkw4YNIywsjPfff58L\nLriAjRs3MnDgQLfpfbHms6/fgz/4fdptY8wy4Ph6x+a47P8CnOptvsEihIhQ6XAQEQBzpSsVKFr6\nhe4LrV3z2dv1kseOHVu7//vf/57333+fpUuXcssttzRI64s1nz2haz63kq7JoFTn0to1n1u7XnJT\nixj5Ys1nT+iaz61U05zUzWXZPqVU4Grtms/erJdcWFjI2rVrOeOMMwgJCWHRokV89913vPDCC7Vp\nfLnms91ux2azYbfbqa6uprKykpCQkAarw+maz62kazIo1bn4Ys3n5tZLrllb2WazMXv2bBITE+nR\nowfz5s1jyZIlHHvssbVpfbXmM8AjjzxCVFQUTzzxBO+++y5RUVE8+uijDcrlyXvwt4Bd8xlgzPr1\nvHb88Yxpw7Y3pQKdrvnsGV3zOUDpmgxKKX/RNZ8DlK7JoJRSvhfYgUHvSlJKKZ8L/MCgNQallPIp\njwKDiHwoIuc511foMHRNBqWU8j1Pv+hfBq4AdojI424W22kXOl+SUkr5nkd3JRljlgPLRSQea5rs\nL53rJ7wOLHQuwtPmtClJKe/1799f5xfrhPr37++zvDy+XVVEugEzgKuANOBdrDmOrgZSfFYiL0QF\nBZFfXd0el1YqYO3du7e9i6A6OI8Cg4gsBoYAC4ALjDEHnKf+JSLttihCdHAwmZWV7XV5pZTqlDyt\nMfzDGFNnQhARCTfGVBpjTvZDuTyifQxKKeV7nnY+P+Lm2GpfFqQltI9BKaV8r8kag4gkA72BSBE5\niSNrOMcBUX4uW7P0dlWllPK95pqSzgFmYq3V/KzL8WLgfj+VyWM68lkppXyvycBgjHkbeFtEfmuM\n+bCNyuSxmvUYlFJK+U5zTUkzjDELgQEi8qf6540xz7p5WZvR9RiUUsr3mmtKinY+xvi7IC2hS3sq\npZTvNdeU9Krz8aG2KY539K4kpZTyveaakv7e1HljzB98Wxzv6HoMSinle801Jf3QJqVooajgYMod\nDowxOveLUkr5iCd3JXVYQSKEBQVR7nAQFQBrsiqlVCBorinpeWPM7SLyCdBglWljzDS/lcxDNXcm\naWBQSinfaK4paYHz8Wl/F6SldL4kpZTyreaakn5wPn4jImFYM6waYJsxpqoNytcsvTNJKaV8y9Np\nt88DXgF2Yc2XNFBEbjDGfObPwnlC50tSSinf8nR21WeAM40xKcaYM4Azgec8eaGITBGRrSKyXUTu\naSTNpSKyWUQ2ichCD8sEaI1BKaV8zdP1GIqNMTtdnu/GmkivSSISBLwETAaygHUissQYs9UlzWDg\nHmCiMaZIRLp7XHq0j0EppXytubuSLnburheRpcC/sfoYLgHWeZD/OGCHMSbdmd8i4EJgq0ua/wPm\nGWOKAIwxed68Aa0xKKWUbzVXY7jAZT8HOMO5nwtEepB/byDD5XkmVrBwdRyAiKzEatp6yBjzuQd5\nA9rHoJRSvtbcXUmzWpm/u+HI9cdDhACDgdOBfsB3IjK8pgbhau7cubX7KSkppKSk6JoMSinlIjU1\nldTU1FblIcY0GLfWMJFIBHAtMByIqDlujLmmmddNAOYaY6Y4n99rvcw84ZLmZWC1MeYd5/PlwD01\nt8q6pDPuynr3rl10Cw3lnn79mn0fSil1tBERjDFezRnk6V1JC4BkrBXdvsFa0a3ZzmesfojBItLf\nOQ7icuDjemk+As4CcHY8H4vVue2RHqGh5FZ1iCEVSinVKXgaGAYbYx4ASp3zJ50HjG/uRcYYO3Ar\n8AWwGVhkjPlFRB4SkfOdaT4HDonIZuAr4C5jTL6nbyA5LIxsDQxKKeUznt6uanM+FojICCAbSPTk\nhcaYZcDx9Y7Nqff8TuBOD8tSR1JYGDk2W/MJlVJKecTTwPCaiCQAD2A1BcU499ud1hiUUsq3PAoM\nxph/OHe/AY7xX3G8l6SBQSmlfMqjPgYR6SYiL4rIBhH5QUSeF5Fu/i6cJ7qHhlJQXY1Nb1lVSimf\n8LTzeRFwEPgt8DsgD/iXvwrljWARuoeGkqv9DEop5ROeBoaexpi/GmP2OLdHgCR/FswbSaGh2pyk\nlFI+4mlg+EJELheRIOd2KeDxtBX+lhwWRo4GBqWU8onmJtErxprCQoDbgZopsYOAEuAuv5bOQ9oB\nrZRSvtPcXEmxbVWQ1tAag1JK+Y6n4xgQkWlYE90BpBpjPvVPkbyXHBZGekVFexdDKaU6BU9vV30c\n+COwxbn90XmsQ9CmJKWU8h1PawxTgVHGGAeAiLwNpAH3+qtg3kjWaTGUUspnPL0rCaCLy368rwvS\nGlpjUEop3/G0xvA3IE1Evsa6Q+l04D6/lcpLOl+SUkr5TrOBQUQEWAlMAMZiBYZ7jDHZfi6bxxJC\nQii126l0OAgP8qYSpJRSqr5mA4MxxojIUmPMSBoustOmSkogJqbh8SAREkNDyamqol9ERMMESiml\nPObpz+sNIjLWryXxwAcfNH5OxzIopZRveBoYxgNrRGSXiPwkIptE5Cd/Fsydt95q/Jx2QCullG94\n2vl8jl9L4aGtW2HnThg8uOE57YBWSinfaLLGICIRInI78GdgCrDfGJNes7VJCV1ceSX885/uzyVp\nU5JSSvlEc01JbwMnA5uAc4Fn/F6iJsyaBe+8A3Z7w3NaY1BKKd9oLjAMM8bMMMa8irVAz2ltUKZG\nnXACdO8OK1Y0PKeBQSmlfKO5wFA7z4QxptrPZfHIrFnuO6GTdFoMpZTyieYCw4kiUuTcioETavZF\npKgtCljfFVfA0qWQn1/3uNYYlFLKN5oMDMaYYGNMnHOLNcaEuOzHtVUhXXXrBued17DWkOQc4KaU\nUqp1AnL+iFtugfnzweE4ciw+JIQqh4NSdz3TSimlPBaQgWHiRIiLgy++OHJMRHT0s1JK+UBABgYR\nq9Ywb17d4zqWQSmlWi8gAwPA9OmwejXs2XPkmHZAK6VU6/k9MIjIFBHZKiLbReSeJtL9TkQcIjLa\nk3yjomDmTHj55SPHtMaglFKt59fAICJBwEtYcy0NB6aLyBA36WKA24A13uR/003W3Unl5dZzrTEo\npVTr+bvGMA7Y4ZxbyQYsAi50k+6vwBNApTeZDxoEY8fCokXWc51hVSmlWs/fgaE3kOHyPNN5rJaI\njAL6GGOWtuQCd9wBzzxj3bqarKOflVKq1TyddrulxM0xU3vSWjb0OeDqZl4DwNy5c2v3U1JSSElJ\n4Ve/grAwazR08mlaY1BKHd1SU1NJTU1tVR5ijGk+VUszF5kAzDXGTHE+vxdrtdAnnM/jgJ1ACVZA\nSAYOAdOMMRvq5WUaK+uiRdatqwu+KOfUtDQyTznFb+9JKaUCiYhgjGn0B7c7/m5KWgcMFpH+IhIG\nXI7LutHGmCJjTKIx5hhjzECszucL6geF5vzud7B/P2RsiOBwdTXF1R1ivj+llApIfg0Mxhg7cCvw\nBbAZWGSM+UVEHhKR8929hCaakhoTEgJ33gnPPCUcFxnJtrKy1hVcKaWOYn5tSvKlppqSAMrKYOBA\nGLNkM1cM6MaM5OQ2LJ1SSnVMHbEpqc1ERVnTZBxYG8VWrTEopVSLdZrAAFZg2LE8mvV5GhiUUqql\nOlVg6NYNrjg1irVZGhiUUqqlOlVgAJh7XSQFkRVs2eZoPrFSSqkGOl1g6NUtmAQTxj3PV7R3UZRS\nKiB1usAAcHJiFN/uLWPz5vYuiVJKBZ5OGRhGxEUxaXoZLjNoKKWU8lCnDAxDoqLoPrqMlSthg1dj\nqJVSSnXKwDA0KoqdVWXMmQN/+hMEyBg+pZTqEDplYBgSFcUvZWVce60hLw8++qi9S6SUUoGjUwaG\n7qGhBAGHjY3nnoM//xkqvVoCSCmljl6dMjCICEOjrKkxfv1rGDIEXnqpvUullFKBoVMGBrCak2rm\nTHr6aXj8ccjNbedCKaVUAOjUgeGX0lJrfwhMnw4PPNDOhVJKqQDQaQPD0OjoOrOsPvQQfPwxrF7d\njoVSSqkA0GkDg2tTEkBCAjz3HFx/Pdhs7VgwpZTq4DptYBgQEcFBm40yu7322KWXQp8+8Mwz7Vgw\npZTq4ELauwD+EizC4MhItpaVMTo2FgARmD8fxo6FSy6BQYPauZCestkgNRX27YPsbGtLToZx4+Dk\nk63qkFJK+UinrTEAjIuNZXVRUZ1jAwfCPffATTcFwIjonBz4619hwACYMwdWroTSUiuiFRTAo49C\nv34wZgwsWgQutSOllGqpTrPmsztvZ2ez9NAh/jV8eJ3jNhtMmAA33GD1OXQ41dVWIJg/36ra3HYb\njBzpPq2w4ktLAAAb5ElEQVTdDp9/Do89ZtUk7rkHrrkGgoPbtsxKqQ6pJWs+d+rAsKe8nElpaeyf\nOBGRup/Lli1wxhnWXUqDB/uypK2Uk2PdWxscDAsXQlKS56/97jv4y1+sYPHOOwHUVqaU8peWBIZO\n3ZQ0ICKCYBF2lZc3ODdsGMyeDb//vfUDvUNYudJqFjr1VFi2zLugAHDaaVZfxO9+B+PHw+uvB0B7\nmVKqo+nUNQaAK7dsYXJCAtf07NngnMMBZ58NZ55p/dBuV99+C7/9rfVL/9xzW5/f5s0wY4Y1uu+t\ntyAiovV5KqUCjtYY3Di9Sxe+LShwey4oyPrOfOEFWL++jQvmauNG61f+++/7JigADB9utZMZA5Mn\n63wgSimPdfrAcFp8PN8WFjZ6vm9fmDfPGuOQn9+GBauxezdMnWrN8verX/k274gIeO89SEmBiRNh\n2zbf5q+U6pQ6fWAYGhVFsd1OZkVFo2kuuQQuuACuvtpqXmozeXlwzjlWZ8ell/rnGkFB1m2tf/mL\n1duuS9oppZrR6QODiHBafDzfNVFrAHjqKau15amn2qhgxsCsWXDRRdagCn+bNQteftlqqvr+e/9f\nTykVsDp9YAA4vZnmJICwMPj3v635lFJT26BQL71k3Zr66KNtcDGn3/wG3ngDzj8fVq1qu+sqpQKK\n3wODiEwRka0isl1E7nFz/g4R2SwiP4rIlyLS19dlaKoD2lXfvrBggTWMYM8eX5fCxcaN8PDDVmdz\nWJgfL+TG+edbb/Kii+B//2vbayulAoJfA4OIBAEvAecAw4HpIjKkXrINwBhjzCjgQ8DnjTknxsSw\nv7KS3KqqZtP++tdw//3W92czlYyWKS2Fyy+H559vvwFo55xjDZ77zW/ghx/apwxKqQ7L3zWGccAO\nY0y6McYGLAIudE1gjPnGGFPTM7wG6O3rQgSLMDE+npUeftPfdps1tuHSS/0w+O3Pf7Zm8bvySh9n\n7KWzz4bXXoPzzoOff27fsiilOhR/B4beQIbL80ya/uK/FvjMHwU5PT6eVA+ak2o8/7x1Q89tt/lw\n8PCqVbBkCfz97z7KsJUuugiefdaqQWzf3t6lUUp1EP4ODO5G27n9mhWRGcAY/NCUBHBBt24szsvD\n4eG3fEgI/OtfVjP8Y4/5oABVVdaMfc89B126+CBDH7niCmt5u7PPhoyM5tMrpTo9f6/HkAn0c3ne\nB8iqn0hEfgXcB5zubHJya+7cubX7KSkppKSkeFyQETExJISEsLKwkNM9/GKOi7MmLj3tNGv/tts8\nvlxDzzwD/ftbgyY6muuus6bxPvtsayK+7t3bu0RKqRZKTU0ltZW3Vvp1riQRCQa2AZOBA8D3wHRj\nzC8uaU4CPgDOMcbsaiKvFs2V5Orx9HT2VlTwyvHHe/W6vXvh9NPhkUesSfe8tmuXNand+vXW2god\n1f33w5dfwldfWZFQKRXwOtxcScYYO3Ar8AWwGVhkjPlFRB4SkfOdyZ4EooEPRCRNRD7yV3kuT0zk\nw7w8qrwc3jxgAHzxhbXUwX/+4+VFjbEGsN17b8cOCmCNqTj5ZLjwQmhipLhSqnPr9LOr1jdpwwbu\n79+f87p18/q1P/5oDRx+5hmrad4jH34Ic+dCWprVcdHR2e3WrKxlZVYUDA1t7xIppVqhw9UYOqIr\nkpJ4PyenRa8dNQqWL7fuOH3zTQ9eUF4Od91l3YUUCEEBrAWC3n7bWubummvaePIopVRHcNQFhkt6\n9ODTQ4cobeH6yMOHw9dfW5WAefOaSfz009bCO2ee2aJrtZuwMKu2kJ4Of/yjLvaj1FHmqGtKApiy\ncSMzk5O53NsV0lzs2WPdxHPJJVandFD9EJuRYVUx1q+HgQNbV+D2UlhoBbVzz23bOZ2UUj6jTUke\nuiIpifcPHmxVHgMHWuPVUlOtJvnKynoJ7r4bbr45cIMCQHy8db/uRx9Z0U8pdVQ4KgPDRd27821h\nYZNrNHiiRw/rzk6bzZpj6dAh54mVK63t3ntbX9j2VvMm33nH6nVXSnV6R2VgiAsJ4drkZJ70wUjf\nyEhrhPQpp1h3em5Y77Da5Z94AqKjfVDaDiA52QoO8+ZZ04UrpTq1o7KPASC7spJh69axZexYksPD\nfZLnBx/Aymve5C9J/yBxx/9AvGrW6/j27LH6HO68s5XDwJVSbaUlfQxHbWAA+MOOHYSJ8PTgwb7J\nsLgY26DjmR65hK7njOW55zyrNFTZq9idv5sdh3aQWZTJgZIDHCg+wKHyQxRXFVNcWUxJVQnVjmqq\nHdXYjZ2QoBDCgsMICw4jMiSS+Ih4ukR0oUt4F3pE9yAxOpGk6CR6xfaib3xfesX2IiTIB7fM7t0L\nZ51lBYY77mh9fkopv9LA4KXMigpOWL+ebePG0cMXC+bcdx9kZVH04tvccgusWwfvvQejRx9JkleW\nx7r96/jhwA9sOLCBjTkbySzKpF98PwZ3HUzfuL70jOlJz9iedI/qTlx4HHHhccSExRASFEKwBBMc\nFIzdYafSXkmVvYoyWxmFFYUUVhZSUFFAbmkuOaU55JTmkFWcRUZhBgdLD5IUk8SghEHW1nUQx3c7\nniHdhzC462DCQ7yoNe3bZwWHG2+0xmkopTosDQwtcNP27XQJCeFvxxzTuox274Zx4+Cnn6BXL8AK\nCn/8cwlTb/2K+BNXkJr+NemF6Zzc62TG9BzDmJ5jGJU8imMSjiE02L8jjG12G5lFmezO382u/F3s\nOLSD7Ye3szVvK+kF6fSL78fIpJGM6DGCEYkjGJU8ikFdBxEkjXRDZWZawWH6dGtQR2drNlOqk9DA\n0ALpFRWMXr+eHePH07U10z9cfLHV+3z//eSW5vLfrf9lybYlfLv3O0JzxhN98Fc8fmMKl502xjdN\nOj5UZa9ix6Ed/HzwZ34++DObDm5iY85G8sryGJk4kpOST2JMrzGM7jma4T2GHwliBw9aYxzGjbM6\npYOD2/eNKKUa0MDQQjds20aICPOOO65lGXz2GY4/3MZH/36Yf25dxDfp33Du4HO5aMhFnDv4XGLD\n4nnjDWvy0ptush4jInz7HvwhvzyfjTkbSTuQxobsDfyQ9QPphemMTBzJ2F5jGdt7LBPihnHsdfcg\nXbtay4X6qCNfKeUbGhhaKN9mY+S6dSwYOpQzExK8em169jaix0zglnMd5J1+MledcBW/HfpbYsNj\nG6TNyoJbb7Vam555BqZNC7wWmJKqEjYc2MD6rPV8v/971mWto7goj/98HEHfynC2vfoYJ434NUkx\nLR9VrpTyHQ0MrbD00CFu3bGDn04+mRgPJrxblbGKJ//3JOPf/ILzSnsT9clnDO7q2d1NX3wBt98O\nffpYS4gOG9ba0rev3NJcvs9YQ9xDjzNoRRoXXxlCTv9uTOgzgQm9JzCx70RGJY8iLNgHHfxKKa9o\nYGilWVu3EhUU1GiTkjGG5buX89jKx9hbsJe/9p/Flde9gGxIg3793L6mMTYbzJ9vzTRxwQUwZ461\nwFvA++c/MXffTeaLj7Li+DBWZ65mdeZqdh7eyYlJJzKxz0Qm9LGCRZ+4Pu1dWqU6PQ0MrVRgszFy\n/XreHjKEs+o1KX2b/i33Lr+X/Ip87jv1PqaPmE7oRRfDpEmtmvqioMBqVpo/35pz6Z57am9qClwr\nV1qzC958s9WhEhxMcWUx67PWsyZzTW2wCA8OZ2LfiUzoPYEJfSYwuudoIkMj27v0SnUqGhh8YNmh\nQ1y3bRv/Gz2a/hER/Jj9I/d/dT+/5P3CX8/8K9NHTCc4KNiaB+Ohh6zVe3wwBiInBx5/3FoK4ZJL\nrDn4Bg3ywRtqL1lZ1q2s4eFWp3RiYp3Txhh25+9mdeZq1mauZc3+NWzJ3cLQ7kOZ0GcC43uPZ0Kf\nCQzuOhgJtI4YpToQDQw+8vfMTF7K3Me43IUs37aY2afP5vox1x9pIz9wwJpS+9NPYexYn147N9da\n1+fll2HyZGvapYkTA6+TGoDqanjwQWsCvrfftt5QE8pt5aRlp7Emcw1rMtewdv9aiiuLGdd7XO02\nttdY7dhWygsaGHygyl7F39f+nQf2ZpCQfAY/jD+TntEuzUrGWJ0Co0fDww/7rRxFRfDPf8KLL0Jc\nnDUDxaWXQlSU3y7pP198Addea31uTz4JMTEevzS7JJvv93/P2sy1rMtax7qsdcSGxTK291jG9hpb\nO1gwIdK7u8mUOlpoYGilr/d8zc1Lb+aYhGN49uxneSzXQX51NR8OH05ozUo8b75pDeZas8YnTUjN\ncThg2TJrYtPVq+Gyy6zv2DFjAqwWUVBgza2UmgpvvGGNmm4BYww7D+9kfdZ61mWtY33WetKy00iM\nTqwdTT6652hO6nkS3aO6+/Y9KBWANDC0UE5JDnd9eRffpn/LC1Ne4MLjL0REsDkcXLJ5M0V2O/8Z\nPpyuWVnW6OYVK2DkSL+UpSmZmVYt4s03rem+p0+3toDqi/j0U7jlFqt97OmnrXt2W8nusLP90HZ+\nOPADP2T9QFp2GmnZacSHx3NSz5MYlTSKUcmjODH5RAZ0GdD4NB9KdUIaGLxkjOHNtDe576v7mDlq\nJg+e8SAxYXWbOezGcM+uXXycl8enDz/McWeeCX/+s0/L4S1jrNrDe+/Bv/8NAwbAb35jbUOGtGvR\nPFNWZvW0z59vTeF9++1WpPMhh3GwJ38PadlpbMzeyI85P5J2II2iyiJGJo3khMQTOCHpBEYkjmBk\n0ki6RHTx6fWV6ig0MHhh+6HtXP/J9ZTZynj9gtc5MfnExhMbwxsPP8z9Y8fy1rhxTO3ecZooqqut\n1pn//tdagTMuDs47D6ZOhVNPbZPWrpbbvdsKsmvWwF/+YrWR+XlKjcPlh9mUY80FtSlnE5sObmJz\n7ma6RHRhWI9hDO8xnOE9hjO0x1CGdh+qfRcq4Glg8ECVvYqn/vcUz615jgdOf4Bbx91q3X7alL/9\nDRYv5ttPP+XqPXuYFBfHc4MH+2aqbh9yOGD9evjsM1i6FLZuhTPOsJrzzzoLRoyAoI7YivLDD9bd\nS5s3WwM5rr66TXvZHcZBekE6W3K3sDl3M5tzN7Mldwtb87YSHRrN0B5Da6coP77b8RzX7Tj6d+nf\n4SZDVModDQzNWJu5lus+uY5+8f2YP3U+/bt4MNR4yRKrTXztWujdm1K7nQf37GFhTg5PDhrEjKQk\ngjtoL3BurtUdsmKFtTJnYaE1Hu/UU63tpJM62Jx3q1bBU09ZA+Suv9763NtxtJ8xhsyiTLbmbWXb\noW1sy9vG1kNb2XFoB9kl2QzoMoBjux3L4ITBDO5qbYO6DqJffD+d/kN1GBoYGlFYUcgDXz/AB1s+\n4Nmzn+XyEZd7Nmjq//0/mDnT6jAdP77OqXVFRfxh504KqquZ3b8/l/XoQUiH/Dl+REYG/O9/1vfu\nypWwY4dVixg3zupTP+kkGDoUWjP7uE/s3GkN5liwAE47DWbNstrHOlANrdxWzq78Xew8vLPOtqdg\nD5lFmfSM6cnAhIEM7DKQAV0GMKDLAPrH96dffD/6xPXx+/obStXQwFCPMYYPtnzAHZ/fwdTBU3n8\nV4/TLaqbZy9+5x1r+PGSJQ2Cgmv+y/PzeTg9neyqKm7t3ZsrExPp3oG+wJpSWgobNsD331tNUGlp\n1uJsw4ZZN12NGGE9DhsGvXu3w+2xJSXwn//AW2/BL79Ya1789reQktIBolfjbHYbGUUZ7Mnfw56C\nPezJ38Pewr2kF6STXphOTkkOidGJ9IvvR9/4vvSN60ufuD70ju1Nn7g+9IrtRc/YnlrrUD6hgcHF\n5oOb+dMXfyKrOItXznuFSf0meX6xZ5+1pj39/HPrJ3QzjDF8V1jIa1lZfHroEJMTErgqKYlfJSR4\nNFNrR1JSAps2wc8/W4+bNll9FSUlcPzx1jZ48JHtmGOs2S78HjR27YIPP4TFi62qzrnnwtlnw69/\nDT17+vnivmWz26wlV4syyCjMILMo09qKrces4ixySnLoEtGFnrE9a5d67RnTk+SYZJKik0iOSSYx\nOpHE6EQSIhP0FlzVqA4ZGERkCvA8EAS8YYx5ot75MOAdYAyQB1xmjNnnJh+PAkNuaS4Pfv0gH/7y\nIfefdj+3jL3F82r7gQPWrZObNlmjyrycMRWgsLqafx08yL8OHuT74mImxcUxtVs3zujShRHR0R22\nP6I5BQWwbZu17dxpbTt2wN69Vs1jwABrdth+/aBvX+uxd2+ri6B3b4iN9WHwyMiweti//NLqPOnV\ny+o0OeUUaxs0KMBG/zVkd9jJLcvlQPEBDpQcqA0W2SXZZJdmk12STW5pLgdLD1JSVULXyK70iO5B\nj6ge9IjuQffI7nSP6k63qG50i+xW+9g1sisJkQnEh8c3f9OF6hQ6XGAQkSBgOzAZyALWAZcbY7a6\npLkJGGmMuVlELgN+Y4y53E1eTQaGvLI8XljzAq/88ApXjLiCOSlz6BrZ1bOC2u3w6qvW3NfXXw+z\nZ/vkvvqi6mq+zM/ns0OHWFlYyIGqKsbHxTE2NpYTYmI4MTqawZGR7dI3kZqaSkpKik/yKi62AsS+\nfdaWkWE9ZmXB/v3WZgwkJ1s/7pOSrFpGjx7WY/fuR7Zu3SAhAaKjPfxut9utNrBVq6wOlFWrrOrN\nqFFWp8kJJ1iDO44/3srYT3z5eXqryl5FXlkeuaW55Jblkluay6HyQ+SV5ZFXlsfh8sMcKj/EobJD\n5Ffkc7j8MEWVRcSFx9ElogsJEQl0iehSZ4sPjycuPI74COuxZosNiyU2PJbYsFhiwmL80lfSnp9l\nZ9QRA8MEYI4x5lzn83sB41prEJFlzjRrRSQYyDbG9HCTl9vAsDt/N39f+3fe2fgOlwy7hLsn3c2g\nrh4OBT5wwJr58623rG+kV16B4cNb9F49kVdVxaqiItJKSvippISfSkvJrKykf3g4gyMjGRwZSf+I\nCPqEh9M3PJxe4eEkhoYS4Ye1lOfOncvcuXN9nq87xljf1dnZR7bcXGvJ6IMH4dAhyMuztsOHrc1m\ns77HExKgSxeIj6+7xcVZtZCaLSbmyBZXnkN8+kZidvxI+LafCNq5Ddm61Qr2gwZZ1ZuBA+tWa3r1\nsiJVC5v+2vLz9AW7w05BRUHtll+RT0FFAYUVhRRUFFBUWURhZWHtY3FlMUWVRRRVFlFcVUxxZTHF\nVcWEBoUSExZDTFgM0WHR1mNoNNFh0dZjaDRRoVFEh1mPUaFRRIZEEhkaWWc/IiSCyBDr8fVnX+fu\nv9xNREgE4cHhhIeEExoUqrPstlBLAoO/G8B7AxkuzzOBcY2lMcbYRaRARLoaYw67y7DaUc2uw7v4\naOtHfLDlA/YV7uOqE65i002b6B3Xu/GSOBzWz9iNG62psletgnXrrA7NV1+1miL8/IfXPSyMad27\nM81lgFy53c7uigp2lpezs7ycvRUVrCwsJKOykqzKSnJtNiKCgkgMC6NrSAhdQ0PpGhJCfEgIccHB\nxIWEEBscTHRwMDHOx8igoNotIiiIcJctTIRQ6w/Fr+/VlciRL/Bjj/XsNZWVVoAoLLSasfLzrf2i\noiOP6elWbaWoyGrOKimxttLSJEpLz3ZuVqUiMsJwjD2LY3fvYdDePQz4eg+9HT+QaP+ERNt+ulVm\nEWPLpyI0juLIRMoiu1IR2ZXKqASqohOojoqjOioOR3QsjqgYTFS0Va2JioKoKHavzeWrf+whKCqC\noMhwgiLDCYkKIygshJBQISQEgoNp8NjYFhRU91HEt3+ewUHBVvOSpzdjuGGMoby6nNKqUkptpRRX\nFlNqK6W0qpSSqhLKbGWU2cpqj5VXl5Ndkl27X15dTrnNeqyorqCiuoJyWzlZP2Xx3zf+S7mtnEp7\nJRXVFdgddsJDwmsDRVhwWJ19d1toUKj1GBx6ZD8olNDgUEKCQmr3Q4Ocz53H3W3BEmw9BgU3eB4s\nwU0+BklQnf2a57X7Lscb2wRp08Do78Dg7p3U/0aqn0bcpAHgqxFRVFRXEB4cwbnRPbgmthfdIsYS\n9M02+Pv11k9Th8P6uVlVZW0FBdZP0oICqw1j1Chru+EGa7hwO09XGhkczPDoaIZHR7s9b4yhsLqa\ngzYbh202DldXc9hmo8hup7C6miK7nf2VlZTa7ZQ6HJTa7ZQ7HJQ7HyscDiodDiqNodLhwGYMVQ4H\nVenpPJKaSogIIc5gUbMf7NyCoM5+kMu+OJ/XecT6dSLO867nav6Ra/646xxzOefuj4FIa5N6Qxpc\n/6MIEOPcGn6IYHeAww759jjWOk5ktf1EHM5jNeccDqDKhtgqwWZDqquQahtSbSPIXo3YqwlyVCMO\nO+LIJ8iRhxg7QcbOnurdHDi8BDnkIMg4EKzHIOPAIUEYxLm57rvbqN0HwZia/wzifCvi3BXXt1fn\n+ZHdep+muDnW4HzTxFkGqX+wSdHOremMw7C2IvsS+tkvtP6gggBna5Wp/Vow4DDgAGxHjtkxlAPl\n4pKuziuNS1kbprFUY8TmUrD6548cNw3SAOLuq6t+yoZ5Ncpdfh79+G9dEGmLpqS5xpgpzufumpI+\nc6apaUo6YIxJdJNXYNw+pZRSHUxHa0paBwwWkf7AAeByYHq9NJ8AVwNrgUuAFe4y8vaNKaWUahm/\nBgZnn8GtwBccuV31FxF5CFhnjPkUeANYICI7gENYwUMppVQ7CZgBbkoppdpGQAyXFJEpIrJVRLaL\nyD3tXZ5AJyJ7RWSjiKSJyPftXZ5AIiJviEiOiPzkcixBRL4QkW0i8rmIxLdnGQNJI5/nHBHJFJEN\nzm1Ke5YxUIhIHxFZISJbRGSTiPzBedzrv88OHxicg+ReAs4BhgPTRSQQlqPpyBxAijHmJGNM/duH\nVdPewvpbdHUvsNwYczxWH9l9bV6qwOXu8wR41hgz2rkta+tCBahq4E/GmGHAROAW53el13+fHT4w\nYI172GGMSTfG2IBFwIXtXKZAV3N3qfKSMWYlkF/v8IXA2879t4GL2rRQAayRzxNae7/lUcgYk22M\n+dG5XwL8AvShBX+fgfDl4G6QXBMj2ZQHDPC5iKwTkf9r78J0AonGmByw/nMCDUbuK6/dIiI/isg/\ntGnOeyIyABgFrAGSvP37DITA4MkgOeWdU4wxJwNTsf4DntreBVLKxXxgkDFmFJANPNvO5QkoIhID\n/Af4o7Pm4PX3ZSAEhkzAdZrTPlgT8qkWcv5qwBiTC/yXhtOUKO/kiEgSgIgkAwfbuTwBzRiT6zIx\n2uvA2PYsTyARkRCsoLDAGLPEedjrv89ACAy1g+ScU3RfDnzczmUKWCIS5fxFgYhEA2cDP7dvqQKO\n64weYP09znTuXw0sqf8C1aQ6n6fzy6vGxejfpzfeBLYYY15wOeb132dAjGNw3q72AkcGyT3ezkUK\nWCIyEKuWYLAGOL6rn6fnROQ9IAXoBuQAc4CPgA+AvsA+4BJjTEF7lTGQNPJ5nonVPu4A9gI31LSR\nq8aJyCTgW2ATOKfdgvuB74F/48XfZ0AEBqWUUm0nEJqSlFJKtSENDEopperQwKCUUqoODQxKKaXq\n0MCglFKqDg0MSiml6tDAoJRSqg4NDEopper4/zgEl3XNVlMtAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fcef6004910>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "N=100; r=20.\n",
    "parms=np.array([[1,0.5],[3,0.5],[9,2],[0.5,1]])\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "for i in xrange(parms.shape[0]):\n",
    "    a=parms[i][0]; b=parms[i][1]\n",
    "    ps=[stats.gamma.pdf(i,a,0,1./b) for i in np.linspace(0.,r,N)]\n",
    "    ax.plot(np.linspace(0.,r,N),ps,label='$a$='+str(a)+',$b$='+str(b))\n",
    "\n",
    "handles, labels = ax.get_legend_handles_labels()\n",
    "ax.legend(handles, labels)\n",
    "plt.title('Gamma Distribution')\n",
    "plt.ylabel('Probability')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Posterior probability given Normal likelihood (fixed mean) and Gamma prior\n",
    "\n",
    "$$ \n",
    "\\begin{align}\n",
    "p(\\lambda| x, a, b; \\mu) = & Const1 \\cdot Gamma(\\lambda| a,b) Norm(x|\\lambda; \\mu) \\\\\n",
    "= & Const2 \\cdot \\lambda^{a-1} \\exp(-b \\lambda) \\sqrt{\\lambda} \\exp \\left( -\\lambda \\frac{(x-\\mu)^2}{2} \\right) \\\\\n",
    "= & Gamma\\left( \\lambda; \\tilde{a}, \\tilde{b}  \\right),\n",
    "\\end{align}\n",
    "$$\n",
    "where \n",
    "\n",
    "\n",
    "<font color = \"red\"> $\\tilde{a} = a + \\frac{1}{2}$ and $\\tilde{b} = b + \\frac{(x-\\mu)^2}{2}$ </font>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "| Distribution |Fixed Parameters | Model Parameters | Conjugate Prior |\n",
    "|---     | ---   | ---    | --- |\n",
    "| Binormial | - | $p$  | Beta |\n",
    "| Bernoulli | - |$p$ | Beta |\n",
    "| Multinomial | - | $p_1,\\cdots, p_n$| Dirichlet|\n",
    "| Normal | $\\sigma^2$ | $\\mu$ | Normal |\n",
    "| Normal | $\\mu$      | $\\sigma^2$ | Inverse Gamma|\n",
    "| Normal | $\\mu$     | $ \\lambda$  | Gamma |\n",
    "| Multivariate Normal | $\\Sigma$ | $\\mu$ | Multivariate Normal|\n",
    "| Multivariate Normal | $\\bf \\mu$ | $\\Sigma$ | Inverse Wishart|\n",
    "| Multivariate Normal | $\\bf \\mu$ | $\\Lambda$ | Wishart|\n",
    "| Multivariate Normal | - | $\\Lambda, \\bf \\mu$ | Normal-Wishart|\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "celltoolbar": "Slideshow",
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}