{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A $\\Pi$-Day Notebook for NLP\n",
    "========\n",
    "\n",
    "- Jacob Eisenstein\n",
    "- 3/14/2018\n",
    "- For CS4650/7650 at Georgia Tech\n",
    "\n",
    "To celebrate $\\pi$-day, this notebook describes one way in which $\\pi$ arises in natural language processing."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy.stats import beta, binom\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Generative models\n",
    "\n",
    "Classical NLP relies heavily on *generative models* of the form $p(x \\mid y)$, where $y$ is a label, and $x$ is a word or vector of word counts:\n",
    "\n",
    "- In a Naive Bayes classifier, we compute $p(\\vec{x} \\mid y)$, where $y$ is a label and $\\vec{x}$ is a bag-of-words vector.\n",
    "- In a Hidden Markov Model, we compute $p(x \\mid y)$, where $y$ is a tag (such as a part-of-speech), and $x$ is a word. \n",
    "- In a trigram language model, we compute $p(w_m \\mid w_{m-1}, w_{m-2})$, where $w_m$ is token $m$ in a sequence.\n",
    "\n",
    "Each of these models involves estimating a parameter vector $\\vec{\\theta}_y$, where $\\theta_{y,w}$ is the probability of word $w$ given $y$, which could be a label, a tag, or a bigram context. A persistent problem with all these models is *sparsity*: if the word $w$ is rare, then the parameter $\\theta_{y,w}$ will be difficult to estimate accurately from finite data. \n",
    "\n",
    "For example, the word \"bicycle\" appears only once in the [Universal Dependencies](http://universaldependencies.org/) English treebank; it is tagged as a noun, among 35000 other nouns. \n",
    "\n",
    "- Suppose that the true probability $P(\\text{bicycle} \\mid N) = \\frac{1}{35000}.$ \n",
    "- Suppose that the nouns in the UD treebank data are a random sample. \n",
    "- If we drew another random sample, how often would we get a count of $1$ for bicycle as a noun?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD8CAYAAAB5Pm/hAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAEPZJREFUeJzt3X+sX3V9x/Hna+WHBn+Acse6tqxE\n2RY0sZAbhsEsDqIiOIuJMsimxJDUJZhgMNPiP+oyEkymOJONpQ5m2ZxIREIDTGWAMSQC3mJBoDLv\nsIY2hV7lhxAjS+G9P+6HecPa3u+93/u93/XT5yP55vs5n/M557xPCK978un5npOqQpLUr98adwGS\npNEy6CWpcwa9JHXOoJekzhn0ktQ5g16SOmfQS1LnDHpJ6tzAQZ9kRZIfJrm5LZ+Q5J4k00m+nuSI\n1n9kW55u69eOpnRJ0iAOW8DYS4DtwGva8ueAK6vquiT/CFwEXNW+n6qqNyY5v437swPt+Nhjj621\na9cutHZJOqRt3br151U1Md+4DPIIhCSrgc3A5cClwJ8CM8DvVNXeJG8FPlNV70ry7db+fpLDgMeB\niTrAgSYnJ2tqamqgE5MkzUqytaom5xs36NTNF4FPAC+25dcDT1fV3ra8E1jV2quAxwDa+mfaeEnS\nGMwb9EneA+ypqq1LeeAkG5JMJZmamZlZyl1LkuYY5Ir+dOC9SXYA1wFnAH8HHN2mZgBWA7taexew\nBqCtfy3wi5fvtKo2VdVkVU1OTMw7xSRJWqR5g76qLquq1VW1FjgfuKOq/hy4E3h/G3YhcFNrb2nL\ntPV3HGh+XpI0WsPcR/9J4NIk08zOwV/d+q8GXt/6LwU2DleiJGkYC7m9kqr6LvDd1n4UOHUfY34N\nfGAJapMkLQF/GStJnTPoJalzBr0kdW5Bc/T6/2HtxlvGXcKy23HFOeMuQTpoeUUvSZ0z6CWpcwa9\nJHXOoJekzhn0ktQ5g16SOmfQS1LnDHpJ6pxBL0mdM+glqXMGvSR1zqCXpM4Z9JLUOYNekjo3b9An\neUWSe5Pcn+ShJJ9t/V9J8tMk29pnXetPki8lmU7yQJJTRn0SkqT9G+R59M8DZ1TVc0kOB+5K8u9t\n3V9V1TdeNv7dwInt80fAVe1bkjQG817R16zn2uLh7VMH2GQ9cG3b7m7g6CQrhy9VkrQYA83RJ1mR\nZBuwB7itqu5pqy5v0zNXJjmy9a0CHpuz+c7WJ0kag4FeJVhVLwDrkhwN3JjkzcBlwOPAEcAm4JPA\nXw964CQbgA0Axx9//ALL/o1D8bV6krQQC7rrpqqeBu4Ezqqq3W165nngn4FT27BdwJo5m61ufS/f\n16aqmqyqyYmJicVVL0ma1yB33Uy0K3mSvBJ4B/Djl+bdkwQ4F3iwbbIF+FC7++Y04Jmq2j2S6iVJ\n8xpk6mYlsDnJCmb/MFxfVTcnuSPJBBBgG/CXbfytwNnANPAr4MNLX7YkaVDzBn1VPQCcvI/+M/Yz\nvoCLhy9NkrQU/GWsJHXOoJekzhn0ktQ5g16SOmfQS1LnDHpJ6pxBL0mdM+glqXMGvSR1zqCXpM4Z\n9JLUOYNekjpn0EtS5wx6SeqcQS9JnTPoJalzBr0kdc6gl6TODfJy8FckuTfJ/UkeSvLZ1n9CknuS\nTCf5epIjWv+RbXm6rV872lOQJB3IIFf0zwNnVNVbgHXAWUlOAz4HXFlVbwSeAi5q4y8Cnmr9V7Zx\nkqQxmTfoa9ZzbfHw9ingDOAbrX8zcG5rr2/LtPVnJsmSVSxJWpCB5uiTrEiyDdgD3Ab8F/B0Ve1t\nQ3YCq1p7FfAYQFv/DPD6fexzQ5KpJFMzMzPDnYUkab8GCvqqeqGq1gGrgVOBPxz2wFW1qaomq2py\nYmJi2N1JkvZjQXfdVNXTwJ3AW4GjkxzWVq0GdrX2LmANQFv/WuAXS1KtJGnBBrnrZiLJ0a39SuAd\nwHZmA//9bdiFwE2tvaUt09bfUVW1lEVLkgZ32PxDWAlsTrKC2T8M11fVzUkeBq5L8jfAD4Gr2/ir\ngX9JMg08CZw/grolSQOaN+ir6gHg5H30P8rsfP3L+38NfGBJqpMkDc1fxkpS5wx6SeqcQS9JnTPo\nJalzBr0kdc6gl6TOGfSS1DmDXpI6Z9BLUucMeknqnEEvSZ0z6CWpcwa9JHXOoJekzhn0ktQ5g16S\nOmfQS1LnDHpJ6twgLwdfk+TOJA8neSjJJa3/M0l2JdnWPmfP2eayJNNJHknyrlGegCTpwAZ5Ofhe\n4ONVdV+SVwNbk9zW1l1ZVX87d3CSk5h9IfibgN8F/iPJ71fVC0tZuCRpMPNe0VfV7qq6r7WfBbYD\nqw6wyXrguqp6vqp+Ckyzj5eIS5KWx4Lm6JOsBU4G7mldH03yQJJrkhzT+lYBj83ZbCf7+MOQZEOS\nqSRTMzMzCy5ckjSYgYM+yauAG4CPVdUvgauANwDrgN3A5xdy4KraVFWTVTU5MTGxkE0lSQswUNAn\nOZzZkP9qVX0ToKqeqKoXqupF4Mv8ZnpmF7BmzuarW58kaQwGuesmwNXA9qr6wpz+lXOGvQ94sLW3\nAOcnOTLJCcCJwL1LV7IkaSEGuevmdOCDwI+SbGt9nwIuSLIOKGAH8BGAqnooyfXAw8zesXOxd9xI\n0vjMG/RVdReQfay69QDbXA5cPkRdkqQl4i9jJalzBr0kdc6gl6TOGfSS1DmDXpI6Z9BLUucMeknq\nnEEvSZ0z6CWpcwa9JHXOoJekzhn0ktQ5g16SOmfQS1LnDHpJ6pxBL0mdM+glqXMGvSR1bpCXg69J\ncmeSh5M8lOSS1v+6JLcl+Un7Pqb1J8mXkkwneSDJKaM+CUnS/g1yRb8X+HhVnQScBlyc5CRgI3B7\nVZ0I3N6WAd4NnNg+G4CrlrxqSdLA5g36qtpdVfe19rPAdmAVsB7Y3IZtBs5t7fXAtTXrbuDoJCuX\nvHJJ0kAWNEefZC1wMnAPcFxV7W6rHgeOa+1VwGNzNtvZ+l6+rw1JppJMzczMLLBsSdKgBg76JK8C\nbgA+VlW/nLuuqgqohRy4qjZV1WRVTU5MTCxkU0nSAgwU9EkOZzbkv1pV32zdT7w0JdO+97T+XcCa\nOZuvbn2SpDEY5K6bAFcD26vqC3NWbQEubO0LgZvm9H+o3X1zGvDMnCkeSdIyO2yAMacDHwR+lGRb\n6/sUcAVwfZKLgJ8B57V1twJnA9PAr4APL2nFkqQFmTfoq+ouIPtZfeY+xhdw8ZB1SZKWyCBX9NLY\nrd14y7hLWHY7rjhn3CWoEz4CQZI6Z9BLUucMeknqnEEvSZ0z6CWpcwa9JHXOoJekzhn0ktQ5g16S\nOmfQS1LnDHpJ6pxBL0mdM+glqXMGvSR1zqCXpM4Z9JLUuUHeGXtNkj1JHpzT95kku5Jsa5+z56y7\nLMl0kkeSvGtUhUuSBjPIFf1XgLP20X9lVa1rn1sBkpwEnA+8qW3zD0lWLFWxkqSFmzfoq+p7wJMD\n7m89cF1VPV9VP2X2BeGnDlGfJGlIw8zRfzTJA21q55jWtwp4bM6Yna1PkjQmiw36q4A3AOuA3cDn\nF7qDJBuSTCWZmpmZWWQZkqT5LCroq+qJqnqhql4Evsxvpmd2AWvmDF3d+va1j01VNVlVkxMTE4sp\nQ5I0gEUFfZKVcxbfB7x0R84W4PwkRyY5ATgRuHe4EiVJwzhsvgFJvga8HTg2yU7g08Dbk6wDCtgB\nfASgqh5Kcj3wMLAXuLiqXhhN6ZKkQcwb9FV1wT66rz7A+MuBy4cpSpK0dPxlrCR1zqCXpM4Z9JLU\nOYNekjpn0EtS5wx6SeqcQS9JnTPoJalzBr0kdc6gl6TOGfSS1DmDXpI6Z9BLUucMeknqnEEvSZ0z\n6CWpcwa9JHXOoJekzs0b9EmuSbInyYNz+l6X5LYkP2nfx7T+JPlSkukkDyQ5ZZTFS5LmN8gV/VeA\ns17WtxG4vapOBG5vywDvBk5snw3AVUtTpiRpseYN+qr6HvDky7rXA5tbezNw7pz+a2vW3cDRSVYu\nVbGSpIVb7Bz9cVW1u7UfB45r7VXAY3PG7Wx9kqQxGfofY6uqgFrodkk2JJlKMjUzMzNsGZKk/Vhs\n0D/x0pRM+97T+ncBa+aMW936/o+q2lRVk1U1OTExscgyJEnzWWzQbwEubO0LgZvm9H+o3X1zGvDM\nnCkeSdIYHDbfgCRfA94OHJtkJ/Bp4Arg+iQXAT8DzmvDbwXOBqaBXwEfHkHNkqQFmDfoq+qC/aw6\ncx9jC7h42KIkSUvHX8ZKUucMeknqnEEvSZ0z6CWpcwa9JHXOoJekzhn0ktQ5g16SOmfQS1LnDHpJ\n6pxBL0mdM+glqXMGvSR1zqCXpM4Z9JLUOYNekjpn0EtS5wx6SercvK8SPJAkO4BngReAvVU1meR1\nwNeBtcAO4Lyqemq4MiVJizVU0Dd/UlU/n7O8Ebi9qq5IsrEtf3IJjiMdUtZuvGXcJSy7HVecM+4S\nujSKqZv1wObW3gycO4JjSJIGNGzQF/CdJFuTbGh9x1XV7tZ+HDhuyGNIkoYw7NTN26pqV5LfBm5L\n8uO5K6uqktS+Nmx/GDYAHH/88UOWIUnan6Gu6KtqV/veA9wInAo8kWQlQPves59tN1XVZFVNTkxM\nDFOGJOkAFh30SY5K8uqX2sA7gQeBLcCFbdiFwE3DFilJWrxhpm6OA25M8tJ+/q2qvpXkB8D1SS4C\nfgacN3yZkqTFWnTQV9WjwFv20f8L4MxhipIkLR1/GStJnTPoJalzBr0kdc6gl6TOGfSS1DmDXpI6\nZ9BLUucMeknqnEEvSZ0z6CWpcwa9JHVuKV4lKElLwtcnjoZX9JLUOYNekjpn0EtS5wx6SeqcQS9J\nnTPoJalzIwv6JGcleSTJdJKNozqOJOnARhL0SVYAfw+8GzgJuCDJSaM4liTpwEZ1RX8qMF1Vj1bV\nfwPXAetHdCxJ0gGMKuhXAY/NWd7Z+iRJy2xsj0BIsgHY0BafS/LIInd1LPDzpanqoOE5Hxo850NA\nPjfUOf/eIINGFfS7gDVzlle3vv9VVZuATcMeKMlUVU0Ou5+Died8aPCcDw3Lcc6jmrr5AXBikhOS\nHAGcD2wZ0bEkSQcwkiv6qtqb5KPAt4EVwDVV9dAojiVJOrCRzdFX1a3AraPa/xxDT/8chDznQ4Pn\nfGgY+TmnqkZ9DEnSGPkIBEnq3EEd9IfaYxaSXJNkT5IHx13LckmyJsmdSR5O8lCSS8Zd06gleUWS\ne5Pc3875s+OuaTkkWZHkh0luHnctyyHJjiQ/SrItydRIj3WwTt20xyz8J/AOZn+Q9QPggqp6eKyF\njVCSPwaeA66tqjePu57lkGQlsLKq7kvyamArcG7n/50DHFVVzyU5HLgLuKSq7h5zaSOV5FJgEnhN\nVb1n3PWMWpIdwGRVjfx3AwfzFf0h95iFqvoe8OS461hOVbW7qu5r7WeB7XT+K+ua9VxbPLx9Ds4r\nsgElWQ2cA/zTuGvp0cEc9D5m4RCTZC1wMnDPeCsZvTaNsQ3YA9xWVb2f8xeBTwAvjruQZVTAd5Js\nbU8KGJmDOeh1CEnyKuAG4GNV9ctx1zNqVfVCVa1j9lflpybpdqouyXuAPVW1ddy1LLO3VdUpzD7l\n9+I2NTsSB3PQz/uYBfWhzVPfAHy1qr457nqWU1U9DdwJnDXuWkbodOC9bc76OuCMJP863pJGr6p2\nte89wI3MTkePxMEc9D5m4RDQ/mHyamB7VX1h3PUshyQTSY5u7Vcye8PBj8db1ehU1WVVtbqq1jL7\n//EdVfUXYy5rpJIc1W4uIMlRwDuBkd1Nd9AGfVXtBV56zMJ24PreH7OQ5GvA94E/SLIzyUXjrmkZ\nnA58kNmrvG3tc/a4ixqxlcCdSR5g9oLmtqo6JG45PIQcB9yV5H7gXuCWqvrWqA520N5eKUkazEF7\nRS9JGoxBL0mdM+glqXMGvSR1zqCXpM4Z9JLUOYNekjpn0EtS5/4HzJaxAdweM98AAAAASUVORK5C\nYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f09f38a80f0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.hist((np.random.rand(35000,1000) < 1/35000.).sum(axis=0),[0,1,2,3,4,5]);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "According to this simulation:\n",
    "\n",
    "- In 35% of the simulated UD treebanks, \"bicycle\" never appears\n",
    "- In 38% of the treebanks, it appears once\n",
    "- In 19% of the treebanks, it appears twice\n",
    "- In the remaining 8% of treebanks, it appears more than twice"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Smoothing\n",
    "\n",
    "Clearly, we cannot naively trust the empirical frequencies that we observe from a treebank. The usual solution is *smoothing*: artificially smooth out the empirical probabilities so that they will be more uniform, and less dependent on the data:\n",
    "\\begin{equation}\n",
    "\\hat{\\theta}_{y,w} \\triangleq \\frac{n_{y,w} + \\alpha}{N_y + V\\alpha}\n",
    "\\end{equation}\n",
    "where,\n",
    "- $n_{y,w}$ is the count of word $w$ with label $y$\n",
    "- $N_y$ is the total count of label $y$\n",
    "- $\\alpha>0$ is some arbitrary smoothing constant\n",
    "- $V$ is the size of the vocabulary. \n",
    "\n",
    "This method of computing probabilities increases bias (since real data is probably not uniformly distributed), but it decreases variance (since the estimates are less dependent on the finite sample that we use to compute them)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 99,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZQAAAEOCAYAAACuOOGFAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAIABJREFUeJzt3Xl8VNXh/vHPmUkmK4EsBCUBwi6b\nbAEKVrQqiloUtbVABREEN6xWW6tWq1at4gruoiiggv22ahVc6k+tKyiryCaCECDIloVAErLNnN8f\nE0LAAFkmuZPM827nlZk7d24eYpInd8699xhrLSIiInXlcjqAiIg0DSoUEREJCBWKiIgEhApFREQC\nQoUiIiIBoUIREZGAUKGIiEhAqFBERCQgVCgiIhIQKhQREQmIMKcDNKSkpCSblpbmdAwRkUZl2bJl\nWdbalsdbL6QKJS0tjaVLlzodQ0SkUTHGbKnOenrLS0REAkKFIiIiAaFCERGRgGjUYyjGmJHA+UAc\nMNNa+6HDkUTEQaWlpWRmZlJUVOR0lEYpMjKS1NRUwsPDa/V6xwrFGPMS8Gtgt7W2Z6Xlw4HpgBt4\n0Vr74NG2Ya39D/AfY0w88AigQhEJYZmZmTRr1oy0tDSMMU7HaVSstWRnZ5OZmUn79u1rtQ0n3/Ka\nBQyvvMAY4waeBs4FugOjjTHdjTG9jDELjrglV3rpHeWvE5EQVlRURGJiosqkFowxJCYm1mnvzrE9\nFGvt58aYtCMWDwQ2Wms3ARhjXgcutNY+gH9v5jDG/13zIPC+tXZ5fWUt3bGDonXr6mvzDaSaP2BH\nrnbED+ZhP6hH/tD+7IfYHFpmyl9rKi07+LyptN3D1vEvMy4DLhcYF8btArcbjMG43eByYVz+ZcYY\n/3OVl7lchx67XP7XhIfrF04Tpv+2tVfXr12wjaGkANsqPc4EBh1j/euBs4DmxphO1trnjlzBGDMZ\nmAzQtm3bWoUq+Pobdtx2W61eK8HJRERgIiJwlX/03zy4PJXuR0RgPEc8Fxl56LmISMISEwhLSsKd\nlERYy2RcMdH6hSYhK9gKpUastU8ATxxnnRnADID09HRbm88Te/pppL3x79q8NDhU919tj1zRHv35\nI9c94rG1ttLLrf/5g7fKz1t76POUP3/Ya60F6/Mv9/rA58X6fODzYb1e8Fn/Mq8P7KFl1ueFimWV\nXuf1YktK8BUXY4tLsMXF2JJifAfvFxfjKynGl1tIWcXjys+VQGnpUb+EJiqKsJYtCUtKqvQx6bBl\n7qQkwhISMGGN+sdP5GeC7Tt6O9Cm0uPU8mWOCouPJyw+3ukYEiRsWZm/lA4cwJuTQ9mePZRlZfk/\n7smqeFy8YQMFixbh27fv5xtxuXAnJBwqnuSWRPftS8yQIYS3bt3w/ygJqI8//phZs2bxyiuv1Gk7\nEyZMYMGCBSQnJ7N69eoApas/wVYoS4DOxpj2+ItkFDCmrhs1xowARnTq1KmumxLBhIVhwsJwRUcT\nlphIROfOx1zfV1REWVY23qw9lO7Zg7eK8ilavZq8N94EwNOuHdFDBhMzZAgxAwfibt68If5ZEkAr\nV66kb9++dd7O+PHjmTJlCuPGjQtAqvrn5GHD84DTgSRjTCZwl7V2pjFmCvBf/IcNv2StXVPXz2Wt\nnQ/MT09Pn1TXbYnUlCsyEk9qCqSmEHWUday1lGzcSMGiRRR8tZC8t99h77zXweUismdPYgb7Cyaq\nbx9cHk+D5peaW7lyJePHj6e4uJirrrqK1q1bc//999d4fG3o0KFkZGTUT8h64ORRXqOPsvw94L0G\njiPiKGMMEZ07E9G5MwnjxmFLSjiwahUFCxdRsHAh2S++SPbzz2MiI4lOT/fvvQwZTESXLv4j2ORn\n7pm/hrU/VfF2Yx10bx3HXSN6HHe97777juTkZM455xyuvPJKLrvssornTj31VPbv3/+z1zzyyCOc\nddZZAc3b0ILtLS8RAYzHQ3T//kT370/L66fgzc+ncPESChYupGDRInY/9BAA7oSE8r2XwcQMHqzx\nlyBQWlrKpk2bGD16NM8//zyDBw8+7PkvvvjCoWT1LyQKRWMo0ti5Y2NpdsavaHbGrwAo3bWLgkWL\nKFy0iIKFi9j37ruAf/wl5pQhJIwfj6eWh8k3FdXZk6gP69atY8CAAeTk5OB2u3/2vPZQGjmNoUhT\nE96qFS1GjqTFyJGHj78sXMTeN99i7xtvknTN1SROmIDRmEuDWrlyJUOGDOGyyy7joosu4pNPPqFV\nq1YVzzflPRS9+SrSyB0cf0kYN442zz1Lx/9+QOyvfsWeadPZdNHFFCxe7HTEkLJy5Up69uxJly5d\nmDp1Kpdeeimlxzh36VhGjx7N4MGDWb9+PampqcycOTPAaQPL2J+dzNZ0paenW83YKKEi/7PP2Hnv\nfZRmZtJ85EiS/3JLkz+fat26dXTr1s3pGI1aVV9DY8wya2368V4bEnsoxpgRxpgZeXl5TkcRaTCx\np51Gh/nvkDh5MnkLFrBp+LnsfeMN/xUDROpBSBSKtXa+tXZyc50gJiHGFRVF8k1/pMNbb+Lp3Ikd\nf72DLWPHUbxhg9PRpAkKiUIRCXURnTvTbs4cTrz/fkp+/JFNF13M7kcfw3fggNPRpAlRoYiECONy\n0eKSi+nw/ns0HzGC7BdeYNOvR5D/2WdOR5MmIiQKRWMoIoeExcfT+oF/0HbObExkJNuuuprMG26k\ndNcup6NJIxcShaIxFJGfixk4kA5vvUnLG28k/9NP2XTueeTMmYMtK3M6mjRSIVEoIlI14/GQdPVV\ndFgwn6j+/dn1jwfIuPR3HFi1yulo0gipUEQET5s2tJnxPCnTHqdszx4yLv0dO++9D28VlwgRORoV\niogA/jPu44YPp8N77xI/Zgy5c+ey6bzzObBypdPRGp2PP/6YsWPH1nk7EyZMIDk5mZ49ex62/IMP\nPqBr16506tSJBx98sM6fJ1BColA0KC9Sfe5mzTjhzjtI+7//w0RGsnXSZIrW/+B0rEYlkBNsffDB\nB4ct83q9XHfddbz//vusXbuWefPmsXbt2jp/rkAIiULRoLxIzUX16knbl1/GFRXF1isnUrJli9OR\nGo2DhVJcXMz48eO5/fbbqc1lroYOHUpCQsJhyxYvXkynTp3o0KEDHo+HUaNG8fbbbwcqep2ExNWG\nRaR2PKkptH1pJlt+fxlbJ0yk3dzXCK905dyg9v6tsDPABxec0AvOPf5bTPU5wdb27dtp06ZNxePU\n1FS++eabav4D6pcKRUSOKaJjR9q88AJbx49n68SJtHvllSZ/kcm60ARbIiLHENWrJ6nPPMO2SZPY\nNvkq2r78Mu7YGKdjHVs19iTqQ31PsJWSksK2bdsqHmdmZpKSklK30AGiQhGRaokZNJCUadPIvP56\nMq+7jjYznscVEeF0rKBT3xNsDRgwgA0bNrB582ZSUlJ4/fXXmTt3bl1jB0RIDMqLSGA0O+NXtH7w\nAQoXL2b7H2/C1nLiqKasvifYCgsL46mnnuKcc86hW7duXHrppfTo4cx0x0cKiQm2Ks0pP2mDLtst\nUmc5c+ey6+/3EnfBCFo/+CDGFRx/m2qCrbrTBFvHocOGRQIrYcwYWt54A/vemc+u++6v1SGx0vRo\nDEVEaiXxqqvw5u0j5+WXcTWPI/mGG5yOJA5ToYhIrRhjSL7lz3j37yP72edwxzUn8YrxTscSB6lQ\nRKTWjDGceM89+PIL2D11Ku64ZrS45BKnY4lDVCgiUifG7Sbloalsy89nx51/wxUTS9zwc5yOJQ4I\niUF5EalfxuMh9YnpRPXpw/Y//5n8L79yOpI4QIUiIgHhio6mzXPPEtGxI5nXX0/h8hVOR5IGFhKF\nosvXizQMd1wcbV98gfDkZLZdfTVF69c7HUkaUEgUis5DEWk4YUlJtH1pJq7oaLZOvJKSjAynIzU4\npybYcnrirZAoFBFpWOEp/sve4/WyZcIESnfudDpSg3Jigq1gmHhLhSIi9SKiQwfavPgCvn372Tph\nImU5OU5HajBOTLAVDBNv6bBhEak3UT160ObZZ9h65SS2XTmJtnNm446NbZDPPXXxVL7P+T6g2zwp\n4ST+MvAvx13PiQm2gmHiLRWKiNSr6AEDSH1iOtuuvY4dd95JymOPYYxxOla90QRbIiL1KPa002h5\nww3seewx9g4ZQvxvf1vvn7M6exL1wakJtoJh4i0Viog0iMQrJ1L49SJ23f8Povv2JaJTJ6cj1Qun\nJtjq2rWr4xNvaVBeRBqEcbloPXUqrpgYtv/xJnxFRU5HqhdOTbAVDBNvhcQEWwelp6fbpUuXOh1D\nJKTlf/El2yZNosWo33Hi3XcHdNuaYKvuNMGWiDQasaf+koSJE9j7+j/Z998PnY4jARQShaJLr4gE\nl+QbbiCyVy923Hknpdu3Ox1HAiQkCkWXXhEJLsbjIeWxR8HnY/vNf8LWcoxBgktIFIqIBB9Pmzac\ncM/dHPj2W/Y89bTTcSQAVCgi4pjm559P899cQvaMGRQsWuR0HKkjFYqIOOqE22/H074922+5hbLs\nbKfjSB2oUETEUa7oaFIefwxf3j5+uu02rM/ndCSpJRWKiDgusmtXkm/9CwWff0HOrNlOx5FaUqGI\nSFCIHz2aZsPOYvfjj3Ng1Sqn49RJME6w1RCTb6lQRCQoGGM48d57CUtKYvtNN+PNz3c6Uq0F2wRb\nDTX5lgpFRIKGu0ULUh59hNKffmLnXXfXalKqYBBsE2w11ORbutqwiASV6H79aDnlOvZMf4KYIUNo\nccnFtdrOzn/8g+J1gZ1gK6LbSZxw++3HXS/YJthqqMm3VCgiEnQSJ0+m4Otv2HnffUT17UNEhw5O\nR6o2TbAlIhJEjNtN64ceYvPIkWz/402k/d8/cUVE1Ggb1dmTqA/BOMFWg02+Za0NmVv//v2tiDQe\n+/73P7u260l2x9/vrdb6a9euredExzdnzhx755132vXr19vu3bvbnTt31ml7mzdvtj169Kh4XFpa\natu3b283bdpki4uL7cknn2xXr1591OXHek1VqvoaAkttNX7HalBeRIJWs9NPJ+Hyy8l97TX2f/SR\n03GqJRgn2Gqoybc0wZaIBDVfSQlbRo+hJDOTDv95i/ATTzzquppgq+5CdoItY0w3Y8xzxph/G2Ou\ncTqPiASe6+Cl7ktL2f6nP2PLypyOJEfhWKEYY14yxuw2xqw+YvlwY8x6Y8xGY8ytx9qGtXadtfZq\n4FLglPrMKyLO8bRr57/U/bJlZD3zrNNx5Cic3EOZBQyvvMAY4waeBs4FugOjjTHdjTG9jDELjrgl\nl7/mAuBd4L2GjS8iDan5iBE0HzmSrGefpeDrwJ9DIXXn2GHD1trPjTFpRyweCGy01m4CMMa8Dlxo\nrX0A+PVRtvMO8I4x5l1gbn1k3ZS3icU7FtfHphuMwVRvPXP89apap6rtG0zFugefN8Yc9X7Fa8pf\nV/k5l3HhwoUxBrdxH1p28FbpOZc5Yj0Orec2bsJd4XjcnsNuYSasWv92cdYJd97BgW+/5adbbqH9\nf94i7IizyMF/5Kr+W9ZOXcfUg+08lBRgW6XHmcCgo61sjDkduBiI4Ch7KMaYycBkgLZt29Yq1Ko9\nq7j/m/tr9VppHAzmUMG4/B8j3BGEu8PxuA7dj3BH4HF5CHeHE+mOpEVkCxIjE0mITPjZLdwd7vQ/\nq8lxxcSQ8vhjZFz6O3bc/ldSn33msPKIjIwkOzubxMRElUoNWWvJzs4mMjKy1tsItkKpEWvtp8Cn\nx1lnBjAD/Ed51ebznJN2Dr9M+WVtXhoULLX/q+PIv1iq2lZVf9VUXs9aS8X/yte1WPz/txXrVl7v\n4MstFp/1VbzWa71Y61/mtd6K5yvfKtaj0nqVXlPqK6XYW0yJt4RSb/l93+H3S7yHbsW+Yv9zZcXs\nK95Hqa+UEm8JRWVF5BbnUuqr+pDQOE9cRbkkRlVdOgeXx3ni9AuwmiK7dSP5llvYdf/95L7yCgnj\nxlU8l5qaSmZmJnv27HEwYeMVGRlJampqrV8fbIWyHWhT6XFq+TJHRYZFEhlW+9aWpstaS35pPjlF\nOf7bgRyyi7IPPS6/bdq7iaVFS9lbvLfKUg53hdM9sTv9WvWjf3J/+iT3oXlEcwf+RY1D/GW/p2DR\nInY9/AhR/fsTVX5ORXh4OO3bt3c4Xehy9DyU8jGUBdbanuWPw4AfgDPxF8kSYIy1dk0dP88IYESn\nTp0mbdiwoU6ZReqizFfG3uK9ZB84vHR2Fuxk5Z6VrMleQ5mvDIOhU3wn+iX3o3+r/vRL7kermFZO\nxw8qZbm5bB55Ea7ISNLeeAN3bIzTkZqs6p6H4lihGGPmAacDScAu4C5r7UxjzHnANMANvGStDdjg\nhU5slGBXVFbEqqxVLNu1jBW7V/Dt7m8pLCsEICU2paJc+rXqR1pcWsi/TVa4ZAlbLh9P8xG/pvXU\nqU7HabKCvlCcoEKRxqbMV8b63PUs37W8omRyinIASIhMqCiXfq360TW+K2GuYHsXu/7teeppsp56\nihMffIAWI0c6HadJUqFUQYUijZ21ls37NrNi1wqW7/aXzPZ8/zBjdFg0vVv2Jv2EdC7qdBEto1s6\nnLZhWK+XrZeP58DatbR/499EaAwl4FQolWgMRZqynQU7WbF7Bct2LWP57uVszN1IhDuCy7pfxhU9\nryDOE+d0xHpXunMnmy8cSVhKa9Jefx2Xx+N0pCZFhVIF7aFIKNi2bxtPfvsk729+nzhPHBN7TWTM\nSWOa/JGK+z/5hMxrryN+3FjH5kJpqkLi4pAi8nNt4trw0NCH+NeIf3Fyy5N5fNnjnP/W+bzxwxuU\n+ZruhRWbnXEG8WPHkjvnFfZ/8j+n44Qk7aGINHFLdi5h2vJpfLfnO9Li0ri+7/UMazesSR4h5isp\nIeN3oyjbsYP2b/+H8FY61DoQtIdSiTFmhDFmRl5entNRRBrcgBMG8Oq5rzLtV9NwGRc3f3Yzo98d\nzdc7vnY6WsC5PB5SHn0UX0kJP/3pz1iv1+lIISUkCsVaO99aO7l5c515LKHJGMOZbc/kzQve5N5T\n7iW7KJtJH05i8oeTWZNdp/OGg05Eh/accOedFC5ZQtZzzzkdJ6SERKGIiJ/b5WZkp5EsuGgBf07/\nM+ty1jFqwSj+9NmfyMjLcDpewDQfeSFxI0aQ9fQzFOpt7gajMRSRELa/ZD+z18xmzto5lHhLuKjz\nRVzT+xqSo5OdjlZn3vwCNl9yMba4hPZvvUlYfLzTkRotHTZcic5DETm2rANZzPhuBv/64V+EmTDG\ndBvDhJ4TGv0FKg+sXkPG6NHEDh1K6lNPNskDERqCBuUr0RiKyLElRSVx+6DbeWfkO5zZ7kxeXv0y\n5755LjNXzaTEW+J0vFqL6tmD5JtvIv/jj8l9rV7m35NKQqJQRKR62jRrw4OnPsi/RvyLPi37MG35\nNCZ9OInsA9lOR6u1hMsvJ+a0oeyeOpWideucjtOkqVBE5Ge6JnTlmbOe4eGhD7Mmew1j3h3DD7k/\nOB2rVowxtH7gAdwtWrD9ppvxFRY6HanJUqGIyFENbz+cWcNnUeorZex7Y/l026dOR6qVsIQEWj/8\nMCUZGey8T9N515eQKBSd2ChSez2TejLv/HmkNU/jD5/8gZdWv1TltM/BLuYXg0i8+iry3nyTvPkL\nnI7TJIVEoWhQXqRuWsW0YtbwWQxrN4zHlz3OHV/d0SgH61tedx1R/fqx8+67Kdmyxek4TU5IFIqI\n1F1UWBSPnPYI1/a+lnd+fIeJ/53Y6AbrTVgYKY88DG4322/+E7ak8ZViMFOhiEi1GWO4ps81PHza\nw6zLWcfod0ezPme907FqJLx1a068/z6KVq9m9+PTnI7TpKhQRKTGhqcNZ/bw2Xh9Xsa+P5ZPtn7i\ndKQaiRs2jPgxo8l5+WXyP//c6ThNhgpFRGqlR1IP5p4/lw7NO3Dj/25k5qqZjWqwPvmWW4jo0oWf\n/nIrpbt2Ox2nSVChiEittYppxcvDX+bstLOZtnwad3x1B8XeYqdjVYsrMpKUxx/DV1zMtkmT8Ooo\n0DoLiULRYcMi9ScqLIqHhz7MtX0ODdZnHchyOla1RHTsSOqTT1CyeTPbrrpaJz3WUUgUig4bFqlf\nxhiu6X0Nj5z2COtz1jPm3TGNZrA+9pRTaP3oIxz47jsyp1yPT0d+1VpIFIqINIxz0s5h1rmzKgbr\nP976sdORqiXu7LM58d57KVi40D/TY1mZ05EaJRWKiARUj8QezPv1PDo278gf//dHXlz1YqMYrG9x\nycUk3/oX9n/4ITvuuqtRZA42YTVZ2RiTBlwHdARygG+B+dZanXIqIhWSo5N5efjL/O2rvzF9+XR+\n3Psjdw+5mwh3hNPRjilx/Hh8+/aR9cyzuJvFkfyXWzSHSg0ccw/FGNPDGPNapUVvA98DTwPDgN7A\n58aYp40xwf2dIiINKjIskqlDp3Jdn+tYsGkBE/47oVEM1iddfz3xl11GzqxZZGtO+ho53lteHwF/\nrfTYba2daa39GMix1k7Cv7eSAcyon4gi0lgZY7i699U8etqj/JDzA6PfHc33Od87HeuYjDG0uv02\nml94AXumP0HOq68d/0UCHL9QzgYqX+v5I2PMlPL7FsBaW2atfRgYXA/5RKQJODvtbGafOxuf9THu\n/XF8vCW4B+uNy8WJ991H7BlnsOu++8h75x2nIzUKxywUa+0qa+3vKy26CWhujFkKtDbGTDbGXGaM\neRoI2qvE6TwUEed1T+zO6+e/TsfmHbnx0xt54bsXgnrg24SHk/L4Y0QPGsRPt93O/k8a1+VlnGBq\n8x/UGBMNnAX0AeKBH4B51tq9gY0XWOnp6Xbp0qVOxxAJaUVlRfztq7/xfsb7nN/hfO4Zck9QD9Z7\n8wvYesUVFK9fT5sZM4j5xSCnIzU4Y8wya236cdcL5r8QAk2FIhIcrLXM+G4GT337FCe3PJnpv5pO\nUlSS07GOqiw3l63jxlG6/Sfazp5FVK9eTkdqUNUtFJ2HIiINzhjDVb2v4rHTH2ND7oagH6wPi4+n\nzYszcScksO3KSRRv2OB0pKCkQhERxwxrN4zZw2djrWXc++P4aMtHTkc6qvBWybR9aSZ4wtk68UpK\nMjOdjhR0VCgi4qhuid2Yd/48OrfozB8//SMzvpsRtIP1nrZtaTtzJr7iYrZeMYHS3brsfWUqFBFx\nXMvolsw8ZybntT+PJ1c8ya1f3EpRWZHTsaoU2aULbWc8T1l2NtsmXol3b1Afi9SgVCgiEhQiwyJ5\n8NQH+UPfP/De5veC+jL4Ub170+bppyjJyGDrVVfhKyhwOlJQUKGISNAwxjDp5ElMO30aG/ZuYNSC\nUazLXud0rCrFDB5M68cepWjVarZNmaLL3qNCEZEgdGa7M5lz7hyMMVz+weVBO1gfN2wYJ95/P4WL\nvuanm28O+cveq1BEJCidlHBSoxisb3HRSFrdfhv7/99H7Ljzb1ifz+lIjgmJQtGlV0Qap6SoJF4a\n/hLndzifJ1c8yV+++EtQDtYnjBtH0pQp5L31FtuuupqSbducjuQInSkvIkHPWsvM1TOZvnw6vZJ6\nMf1X02kZ3dLpWIex1pL76mvsefxxrNdL0rXXknjFeIzH43S0OtOZ8iLSZBhjuLLXlUz71TQ27t3I\n6HdH8+X2L4PqLTBjDAljL6PDe+8SO3Qoex5/nE0XX0xhCP0Rq0IRkUbjzLZn8sq5r+Bxe7jmo2uY\n+OFEVu5Z6XSsw4SfcAKpTz5B6rPPYAsPsOWysfx0+18py811Olq901teItLolHpL+dcP/+L5754n\npyiHM9qcwQ39bqBDiw5ORzuMr7CQrGefJfvlWbhjY0m+5RaaXzSy0U0rrKsNV0GFItK0FJYWMmft\nHGatmcWBsgNc0PECru19LSfGnuh0tMMUrf+BnXffzYEVK4hOT+eEu+8iolMnp2NVmwqlCioUkaYp\ntyiXF1a9wOvfv47BMOqkUUzqNYkWkS2cjlbB+nzsfeMNdj/yKL7CQhInTCDpmqtxRUY6He24VChV\nUKGING078nfw9LdPM3/TfKLDohnfYzxju48lOjza6WgVynJy2D31IfLefpvwNm044W93EnvqqU7H\nOiYVShVUKCKhYWPuRp5c8SSfbPuExMhErup9Fb/p/BvC3eFOR6tQ8PU37LznHko2b6bZucNpdett\nhLdKdjpWlVQoVVChiISWb3d/y7Tl01i2axkpsSlM6TuF89qfh8sExwGuvpISsl98keznnsd4PLS8\n8UbiR4/CuN1ORzuMCqUKKhSR0GOt5cvtXzJ9+XTW566nS3wXbuh3A6emnBo0R1uVbNnCznv+TsHC\nhUT27MkJ99xNVI8eTseqoEKpggpFJHT5rI8PNn/AkyueJDM/k/6t+nNjvxvpk9zH6WiAv/j2vfce\nux54EG9ODvG//z2JV4wnvHVrp6OpUKqiQhGRUm8pb2x4g+dWPkd2UTantzmdCzpeQL/kfiRGJTod\nD+++feyZNo3cea+DtXjS0ogZMpjowYOJGTQId1xcg2dSoVRBhSIiBxWWFvLquleZtXoW+0v3A5AW\nl0a/Vv3ol9yPfq36kRqb6tjbYsWbNlPwxecULFxEwZIl2MJCcLmI7NWTmMGDiRk8hKi+fXA1wLXC\nVChVUKGIyJFKvaWsyV7D8t3LWb5rOct3L2d/ib9gkqOS/QVTXjKd4zs7MqBvS0o48N13/nJZuJAD\nq1aB14uJiiI6Pd1fMEMGE9GlC8YV+HwqlCqoUETkeHzWx8a9G/3lsms5y3YvY3fhbgCaeZrRN7kv\n/ZL70b9Vf3ok9nDkUGTv/v0ULllSUTAlmzYB4E5MJOYXvyBmyGBiBg8O2PhLSBSKMSYG+Ay421q7\n4Hjrq1BEpKastWzP316xB7Ns1zIy9mUAEOGOoFdSL/q16kf/5P70Tu5NTHhMg2cs3bmTgkVfU7Bo\nIQWLFuHdkwUQsPGXoC4UY8xLwK+B3dbanpWWDwemA27gRWvtg8fZzt+BfGCtCkVEGkr2gWxW7F7B\nsl3LWL57Od/nfI/P+nAZF4mRiSREJvhvUQkV96taHhUWFfBs1lqKN2ygcNGiw8ZfEidPJvmmP9Zq\nm8FeKEPxF8Gcg4VijHEDPwAJqJZSAAARd0lEQVTDgExgCTAaf7k8cMQmJgC9gUQgEshSoYiIUwpK\nC1i5eyUr96xkV+EusouyySnKIedADjlFORSWFVb5uqiwqMPKJjHqUOnER8YTHRaNx+05dHP9/H6E\nO4JwdzhhJqzKAwgOjr+EtWyJp127Wv37qlsoYbXaeh1Zaz83xqQdsXggsNFauwnAGPM6cKG19gH8\nezOHMcacDsQA3YEDxpj3rLX1M5nz+7fCzlX1smkRaRwsFmvBZy0+698T8JU/thZ6WksPW76e/wX+\nj8RRZGPJNT5yjI+9xkduxc3L3sIs9prdbHNZvit/3leLA8uMBQ8Qbg3hFR8N4RY8GM6P7MSEdv8K\n4Ffk5xwplKNIASpPxJwJDDraytbavwIYY8bj30OpskyMMZOByQBt27YNVFYRaSR81lLq81HmtZR6\nfZRWfPTfL/P68NpDZVFVaVhbXhJ1FA60LL/5uThynkMflkIX7HdbSoyl1ECpC8qMpYwj7hsoc1lK\ngbKD62EPu39wnQP19Pd2ZcFUKLVirZ11nOdnADPA/5ZXrT7JucccyhGRBravqJTd+4rIyi8hO7+E\n7ILi8vvFZOeXkJVfTHaB/+P+orIqt+EJc5EU4yExNoIoj5uIMBcRYS48YS48bv/HiDC3//ERz0WE\nu4lwH7G8/Lkwt8FlKt1cVNx3u/xTBbuMwW0MxoDbdXBdcFW6b8o/gv++AYwB/z3//cM+Ysqf//lr\nDj6ub8FUKNuBNpUep5YvE5EQt33vAZZszmFxRg5LNuewYXf+z9YxBuKjPSTGeEiM9dCjdRxJsRHl\njyNIjPWQFOshMcZ/Pzai6jEHqb1gKpQlQGdjTHv8RTIKGBOIDRtjRgAjOjWiGdJEQpW1lo278yvK\nY0lGLtv3HgAgNiKM/u3iuaB3a9omRvsLo7wk4qPDCXMHx1WEQ5VTR3nNA04HkoBdwF3W2pnGmPOA\nafiP7HrJWnt/ID+vjvISCT6lXh9rftpXsQeyNCOH3MJSAJJiIxjYPp4BaQkMSEug24lxuF3aq2ho\nwX6U1+ijLH8PeK+B44hIAyosKWPF1r0s3pzDkowcVmzdy4FSLwDtEqM5s1srBqYlMKB9AmmJ0Xpb\nqhEJpre8RKSJ2r2/iFlfZfDVj9ms2Z5Hmc9iDHQ7IY7fDWhTvgcST3Jc8M+vLkcXEoWiMRQRZ+wr\nKuX5z37kpS8zKPX66Nu2BZOHdmBA+wT6t4snLjJ4puSVuguJQrHWzgfmp6enT3I6i0goKCr1MmdR\nBs98+iN7C0u5oHdrbj67C+0SG/46V9JwQqJQRKRhlHl9vLE8k2kfbWBHXhGndWnJLcO70qN1c6ej\nSQMIiULRW14i9ctaywerd/Lwh+vZtKeAvm1b8NilfRjc0fkZEKXhNOrL19eUDhsWCbyFG7OY+sH3\nrMzMo1NyLH8+pytnd2+lo7OakKA+bFhEGr9VmXk89N/v+WJDFq2bR/LQb07mkn6pOk8khKlQRKRG\nNmcV8MiH63n3ux3ER4dzx/nduOwX7YgMdzsdTRymQhGRatm1r4jpH2/gn0u2ERHm4g9ndOLKoR10\n6K9UCIlC0aC8SO3lFZby3Oc/8vJXm/H6LJcNasuUMzrTslmE09EkyIREoeg8FJGaK/P6mPnlZp7+\n30b2F5dxYe/W3DSsK20To52OJkEqJApFRGpmX1Ep189dwWc/7OH0ri255ZyT6N46zulYEuRUKCJy\nmC3ZBUycvZSMrAL+cVEvxgzSTKdSPSoUEamw6MdsrnltGQBzJg5kSMckhxNJYxISs9EYY0YYY2bk\n5eU5HUUkaM1bvJWxM78hMcbDf649RWUiNRYShWKtnW+tndy8ua4nJHKkMq+Pe+av4bY3VzGkUxJv\nXXcKaUm6iKPUnN7yEglh+4pKmTJ3BZ//sIcrTknjr+d10zS6UmsqFJEQlZFVwMTZS9iSXcgDF/di\n9EANvkvdqFBEQtDCH7O49rXlALwycZCuCiwBoUIRCTGvfbOFu95eQ1pSDDMvT9ekVxIwIVEouvSK\niH/w/b531zFrYQand23JE6P76jpcElAhMfqmo7wk1OUdKGXC7KXMWpjBxF+2Z+blA1QmEnAhsYci\nEsoysgqYMHsJW7MLefDiXozS4LvUExWKSBO2cGMW17y2HJeBV68cxC86aPBd6o8KRaSJOjj43j4p\nhpmXD9BVgqXeqVBEmpgjB9+fHN2XZhovkQagQhFpQvIOlDJl7nK+2JDFlb9sz23nddMc79JgVCgi\nTcTm8jPft+UUMvWSXvxugAbfpWGpUESagMMG3ycOYpAG38UBIXEeii5fL03Zq19vYexLi0luFsHb\n1/1SZSKOCYlC0YmN0hSVeX3c9fZq7vjPaoZ2TuLNa4foSC5xlN7yEmmE8gpLmTLPP/g+6dT23Hqu\nBt/FeSoUkUZm0558rpy9lG25hTx0yclcOqCN05FEABWKSKPy1cYsrnl1GWFuF69d+QsGtk9wOpJI\nBRWKSCPxyqIM7p6/lk4tY3nx8nTaJGi8RIKLCkUkyJV5ffx9wVrmLNrCmSclM21UH535LkFJhSIS\nxPIKS7lu7nK+3JjFVUM7cMvwkzT4LkFLhSISpCoPvj/8m5P5bboG3yW4qVBEgtCXG7K49jX/4Pvc\nSb9gQJoG3yX4qVBEgowG36WxColC0Zzy0hiUen38ff5aXvl6C2d1S2baqL7ERoTEj6g0Ebr0ikgQ\nyCssZfzLi3nl6y1cdVoHnh+brjKRRkffsSIOstbyxYYs7npnDdtzD/DIb3vzm/6pTscSqRUViohD\nVmzN5aEP1rNoUzap8VHMnTSIdA2+SyOmQhFpYBt37+fh/67nv2t2kRTr4Z4LejB6YFs8YSHxDrQ0\nYSoUkQby094DTPvoB/69LJNoTxg3DevCxF+2J0ZjJdJE6DtZpJ7lFpTwzKcbmb1oC1i44pT2XHt6\nRxJjI5yOJhJQKhSRelJYUsZLX27m+c82UVBSxsX9UrnxrM6kxuu8EmmaVCgiAVZS5uOfS7Yy/eON\nZOUXM6x7K/58Tle6tGrmdDSReqVCEQkQn88y/7ufePTDH9iaU8jA9gk8P7Y//dvFOx1NpEGoUETq\nyFrLZz/s4aEP1rN2xz66nRjHy1cM4PQuLTFGVwaW0KFCEamD5Vtzmfr+93yzOYe2CdFMH9WHESe3\nxqVLzEsIUqGI1IC1li3ZhSzOyOHDNTv5aN1ukmIjuPfCHvxugM4lkdCmQhE5Bq/Psm7HPpZm5LAk\nI5fFGTns2V8MQEKMh5uHdWGCziURAVQoIocpKvXyXWYeSzJyWLw5h+VbctlfXAZASosoTumYyID2\nCQxMS6Bjy1i9tSVSiQpFQtq+olKWbcllyeYclmTksHJbHiVeHwCdk2MZ0ac1A9MSGNA+gZQWUQ6n\nFQlujbpQjDGnA/cCa4DXrbWfOhpIgprPZ9m9v9hfIOV7IN/v3IfPQpjL0DOlOZcPaceAtAQGpCUQ\nH+NxOrJIo+JYoRhjXgJ+Dey21vastHw4MB1wAy9aax88xmYskA9EApn1GFeCVGFJGdn5JWTlF5Od\nX0J2QTFZ+SUV9w8+l5VfQm5hCV6fBSAq3E2/di34w5mdGZiWQJ+2LYj2NOq/r0Qc5+RP0CzgKWDO\nwQXGGDfwNDAMf0EsMca8g79cHjji9ROAL6y1nxljWgGPAb9vgNxSCz6fpcTro7jMR0mZj+IyLyVl\nPv+yUv/HysuLyw6tW1TqJaeg5PDCKC+LwhJvlZ8vNiKMxFgPiTEe2iRE07dtPEmxHlo2i+Dk1Bb0\naB1HuFtHZIkEkmOFYq393BiTdsTigcBGa+0mAGPM68CF1toH8O/NHE0uUG9X2lu8OYd5i7cedz1r\nbbW2V9VaR7606nV+vtRWWtlisfbQtioeH/Y5bMV9W2mbtnyZz1r/zVfp/sHlvkP3vT7/tr3l61jr\nPyKq8n1/WXgp8foo9Vbva3M0YS5TXhARJMZ6aJ8UQ2KMh8TYCJJiPSTF+pcnxkaQGOMhMtxdp88n\nIjUXbPv4KcC2So8zgUFHW9kYczFwDtAC/95OVetMBiYDtG3btlahsvP977tXR3VPjK5qtSPPqq5y\nU1UsNJVea8ozmPIVD27yyOeptF7FY8DlMriMwW38y8NcLtwu/32XMbhdBpfxb89tDC6Xf7n/duj1\nLgMRYW48YS48YS4iyj963C4iwt1EuI9YHuaqWP/wdf3L4yLDdNa5SJALtkKpEWvtm8Cbx1lnBjAD\nID09vVZ/Jp/b60TO7XVibV4qIhIygu1N5O1Am0qPU8uXiYhIkAu2QlkCdDbGtDfGeIBRwDt13agx\nZoQxZkZeXl6dA4qISNUcKxRjzDxgEdDVGJNpjJlorS0DpgD/BdYB/2etXVPXz2WtnW+tndy8efO6\nbkpERI7CyaO8Rh9l+XvAew0cR0RE6ijY3vISEZFGKiQKRWMoIiL1LyQKRWMoIiL1LyQKRURE6l+j\nPrGxuowxI4ARQKExZl2lp5oDeUfcP9qycCCrBp+28naO99yxHh8ro3LVPFcekFSDbMpVf7mqk0W5\ngiNXu2q9ylobMjdgxtEeH7x/tGXA0rp8rurmOFqGqjIqV81zlX+sdjblqr9c1cmiXMGX61i3UHvL\na/4xHs+v5rLafq7q5jhahqryKJdyNeZc1cmiXMGX66hMeRvJcRhjllpr053OcSTlqrlgzaZcNaNc\nNdMQuUJtD6UuZjgd4CiUq+aCNZty1Yxy1Uy959IeioiIBIT2UEREJCBUKCIiEhAqFBERCQgVSoAY\nY2KMMUuNMb92OstBxphuxpjnjDH/NsZc43Seg4wxI40xLxhj/mmMOdvpPAcZYzoYY2YaY/4dBFli\njDGzy79Ov3c6z0HB9DWqLIi/p4LyZxDq6XdWbU5eaUo34CVgN7D6iOXDgfXARuDWamzn78AtwK+D\nKVf5a1zAq0GYKx6YGYS5/u309xowFhhRfv+f9ZGnLl+7+voaBSBXwL6nApwrYD+DgcoV6N9Z1loV\nCjAU6Ff5PwLgBn4EOgAeYCXQHegFLDjilgwMwz+75PgAFkqdc5W/5gLgfWBMMOUqf92jQL8gzFVf\nhVKTjLcBfcrXmRssPwP1/TUKQK6AfU8FKlegfwYD9P0V8N9Z1trQuJbXsVhrPzfGpB2xeCCw0Vq7\nCcAY8zpwobX2AeBnu4fGmNOBGPz/oQ4YY96z1vqczlW+nXeAd4wx7wJz65IpULmMMQZ4EHjfWru8\nrpkClau+1SQjkAmkAt9Sz29N1zDX2vrMUttc5dfoC+j3VCByAWsD/TMYoFyxBPh3FoTIxSFrIQXY\nVulxJjDoaCtba/8KYIwZD2QF4j9MIHKVF93FQAT1OwtmjXIB1wNnAc2NMZ2stc8FQy5jTCJwP9DX\nGHNbefHUt6NlfAJ4yhhzPrW8DEZ95HLoa3TcXDTc91SNcjXgz2CNcllrp0Dgf2epUALIWjvL6QyV\nWWs/BT51OMbPWGufwP8LM6hYa7OBq53OAWCtLQCucDrHkYLpa1RZEH9PfUoQ/gweFOjfWTrKq2rb\ngTaVHqeWL3OactVMsOaqLFgzKlfNKBcqlKNZAnQ2xrQ3xnjwD16943AmUK6aCtZclQVrRuWqGeUC\nHeUFzAN2AKX431+cWL78POAH/EdI/FW5lKupZlQu5QrUTReHFBGRgNBbXiIiEhAqFBERCQgVioiI\nBIQKRUREAkKFIiIiAaFCERGRgFChiIhIQKhQREQkIHRxSBEHGWN6ANOBtsAr+OfXmWOtXeJoMJFa\n0JnyIg4xxkQCy4HfApuA74Fl1tqLHQ0mUkvaQxFxzlnACmvtGoDyi/c96mwkkdrTGIqIc/oAKwCM\nMa2BfGvtV85GEqk9FYqIc0rwz6gH8AD+Ob9FGi0Viohz5gJDjTHrgZXAImPMNIczidSaBuVFRCQg\ntIciIiIBoUIREZGAUKGIiEhAqFBERCQgVCgiIhIQKhQREQkIFYqIiASECkVERALi/wNxyGOjZpIL\n+wAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f09ef229f98>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "N = 1e5\n",
    "V = 1e4\n",
    "alpha = np.logspace(-4,4,17)\n",
    "for k_i in [1,10,100,1000]:\n",
    "    plt.loglog(alpha,(k_i+alpha)/(N + V*alpha))\n",
    "plt.ylabel(\"$\\\\hat{\\\\theta}$\")\n",
    "plt.xlabel(\"$\\\\alpha$\")\n",
    "plt.legend([\"$k=1$\",\"$k=10$\",\"$k=100$\",\"$k=1000$\"]);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As the amount of smoothing increases, the probabilities converge to the same value, regardless of the counts. This means that bias increases with $\\alpha$, while variance decreases."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Smoothing: a probabilistic justification\n",
    "\n",
    "Additive smoothing can be justified as a *maximum a posteriori (MAP)* estimate of $\\vec{\\theta}_{y,w}$:\n",
    "\\begin{equation}\n",
    "P(\\vec{\\theta}_y \\mid \\vec{x}^{(TR)}) \\propto P(\\vec{x}_{(TR)} \\mid \\vec{\\theta}_{y}) P(\\vec{\\theta}_{y}),\n",
    "\\end{equation}\n",
    "where $\\vec{x}^{(TR)}$ is a vector of word counts from instances with label $y$ in a *training set*. \n",
    "\n",
    "The *likelihood* $P(\\vec{x}^{(TR)} \\mid \\vec{\\theta}_{y})$ follows a multinomial distribution,\n",
    "\\begin{equation}\n",
    "P(\\vec{x}^{(TR)} \\mid \\vec{\\theta}_{y}) \\propto \\prod_w \\theta_{y,w}^{x^{(TR)}_w}\n",
    "\\end{equation}\n",
    "\n",
    "What about the *prior*, $P(\\vec{\\theta}_{y})$? Since $\\vec{\\theta}_{y}$ is a discrete probability distribution, a natural choice is the [symmetric Dirichlet distribution](https://en.wikipedia.org/wiki/Dirichlet_distribution),\n",
    "\\begin{equation}\n",
    "P(\\vec{\\theta}_{y}; a) \\propto \\prod_w \\theta_{y,w}^{a - 1}\n",
    "\\end{equation}\n",
    "where $a > 0$ is a parameter of the prior. \n",
    "\n",
    "Here is an amazing property of the Dirichlet distribution, which explains why it is the right choice for the prior: if the prior is Dirichlet and the likelihood is multinomial, then the posterior $P(\\vec{\\theta}_y \\mid \\vec{x}^{(TR)})$ is also a Dirichlet distribution, with parameters $\\vec{x}^{(TR)} + a$:\n",
    "\\begin{equation}\n",
    "\\vec{\\theta}_y \\sim \\text{Dirichlet}(\\vec{x}^{(TR)} + a)\n",
    "\\end{equation}\n",
    "\n",
    "Furthermore, the modal value of this distribution --- the *maximum a posteriori* estimate --- is simply:\n",
    "\\begin{equation}\n",
    "\\hat{\\theta}_{y,w} = \\frac{x_w^{(TR)} + a - 1}{V(a - 1) + \\sum_{w'} x^{(TR)}_{w'}}\n",
    "\\end{equation}\n",
    "\n",
    "This is identical to the smoothed estimate, when $\\alpha = a - 1$. Let's check that this is really the maximum a posterior estimate."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 103,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from scipy.stats import dirichlet"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 125,
   "metadata": {},
   "outputs": [],
   "source": [
    "V = 10000\n",
    "a = 1.5 # parameter of the Dirichlet prior\n",
    "N = 1e6\n",
    "\n",
    "# generate the true theta\n",
    "theta_true = np.random.dirichlet(a * np.ones(V))\n",
    "# generate some observed counts\n",
    "x_tr = np.random.multinomial(N,theta_true)\n",
    "\n",
    "log_likelihood = lambda theta_hat : np.dot(np.log(theta_hat),x_tr)\n",
    "log_prior = lambda theta_hat : dirichlet.logpdf(theta_hat,a*np.ones(V))\n",
    "log_posterior = lambda theta_hat : log_likelihood(theta_hat) + log_prior(theta_hat)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 126,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZEAAAEVCAYAAAAhANiZAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAIABJREFUeJzt3Xl4leWd//H3NwlhScISSIAQlkBY\nZAtoDIIbKiqgQMe1tr+OVltbtZ1qO53RWjtja5dR55qfTutUq1PHX692tC4tWhTF4lYXFgUh7JsC\nARIIS1hCSPL9/XGeQEQCh0NOnnOSz+u6zpVz7vMs3wOHfLif537ux9wdERGRWKSEXYCIiCQvhYiI\niMRMISIiIjFTiIiISMwUIiIiEjOFiIiIxKzNh4iZjTGz981skZktMLOSJpa738xKzWy5mT1sEZ3M\n7C9mtiJ47xeNlv+umS0zs4/N7HUz6x9FLV3N7Nlge8vNbHxzflYRkebW5kMEuB+4193HAD8KXn+G\nmU0AzgZGAyOBM4Hzg7cfdPdhwFjgbDObErR/BBS7+2jg2WNt9xgeAl4JtlcELI/5U4mItACFCDjQ\nOXjeBShrYpkOQDrQHmgHbHP3/e4+F8Dda4APgfzg9Vx33x+s/35DO4CZfd/M5ge9lHuDti7AecAT\nDdtz913N+klFRJpZWtgFJIDbgdlm9iCRUJ1w9ALu/p6ZzQW2AAb80t0/00sws67ANCK9iaPdBLwc\nLHcJMBgoCbY108zOA/YAFcBvzawIWAh8x933NcunFBGJgzbREzGzOWa29BiPGcAtwB3u3he4g6An\ncNT6hcBpRHoTfYALzezcRu+nAX8AHnb3dUet+3+AYuCBoOmS4PERkZ7LMCKhkgacDvyXu48F9gF3\nNtsfgohIHFhbnzvLzHYDXd3dzcyA3e7e+ahlvg90cPefBK9/BFS7+/3B6/8G9rr7Pxy13iTgP4Hz\n3b08aPt3YJW7P3rUsr2A9919QPD6XOBOd7+s2T+0iEgzaRM9kRMo48hJ8guB1cdY5lPgfDNLM7N2\nwfLLAczsPiLnUm5vvIKZjQUeBaY3BEhgNnCjmWUGy/Uxs1x33wpsNLOhwXIXAcua4wOKiMSLeiJm\n5xA5j5EGVAO3uvtCMysGvunuXzOzVOARIie+ncgIqu+aWT6wEVgBHAw2+Ut3f9zM5gCjiJxHAfjU\n3acH+/wO8LWgfS/wf9x9rZmNAR4ncgJ/HfBVd98Z1z8AEZFT0OZDREREYqfDWSIiEjOFiIiIxKzV\nXyfSo0cPHzBgQNhliIgklYULF25395wTLdfqQ2TAgAEsWLAg7DJERJKKmX0SzXI6nCUiIjFTiIiI\nSMwUIiIiEjOFiIiIxEwhIiISRzfeeCO5ubmMHDky7FLiQiEiIhJHN9xwA6+88krYZcSNQkREWo3S\n0lImTZrEkCFD+MlPfsK3v/1t5s+ff9LbefbZZznrrLMoKirinHPOoaKiIuaazjvvPLKzs2NeP9G1\n+utERCQJVVXB00/D6tUweDBcey1kZR13lerqaq6++mr++Mc/MnDgQIYNG8YZZ5zBmWeeeXiZc889\nl6qqqs+t++CDDzJp0qTDry+44AKuuuoqAO69916eeeYZbrvttmb6cK2LQkREEss778DUqVBfD/v2\nQUYGfPe7MGsWnHNOk6vNmTOHsWPHMmLECABqamr43ve+95ll3n777ahKePLJJ3n66ac5ePAgW7du\n5Wc/+9ln3p80aRJbt2793Ho//elPmTFjRlT7aC0UIiKSOKqqIgHSuLewL7hD9NSpUFYGmZnHXHXR\nokWMHTsWgLKyMjIzMzn77LM/s0w0PZGnnnqKefPm8de//pXMzEzOO++8w8HUYM6cObF+wlZHISIi\nCaGmtp4V//M8i4dfxMfZ/VnbPZ/LVrzNDQtfItXrIz2Tp5+Gm2465vrp6els3rwZgLvuuouamprP\nLRNNT2TJkiVMmDCBzMxMnnvuOd59911GjRp1ah+uFQvtxLqZjTGz981skZktMLOSoN3M7GEzW2Nm\nH5vZ6Y3Wud7MVgeP68OqXUROTV29s3JrFc8s2Mg9f1rKjF++w8h/mc30TT24Z+LXeL2whP3tOvCT\ni27m777yIEtzB0Z6JGvWNLnNL33pS7z11lsMHTqUoqIixo8fz+23397k8k254YYbeOSRRygpKeGj\njz5i4MCBZGRkxPxZr7vuOsaPH8/KlSvJz8/niSeeiHlbiSi0m1KZ2avAf7j7y2Y2Ffgnd58YPP82\nMBUYBzzk7uPMLBtYABQTubvgQuCME935r7i42DUBo0h43J1Pduxn8aZdLNm0m4837WZp2W7219QB\nkNk+jVF9ujC6bxeKVi5k1L/dQ/7WDQC8NOxc7p10Mzs7dubGxbO446oSOn39xhA/TdthZgvdvfhE\ny4V5OMuBzsHzLkTudQ4wA3jKI+n2vpl1NbPewETgNXevBDCz14DJwB9atGoRaZK7s3VPNYs37mbJ\n5l18HITG7gOHAGiflsKIvM5cU9yX0fldGJ3flYE9MkhJscgGzsmHH916eHvTVrzNees/5BcTv8pv\nxk5j1vb23LeynAuG5obx8eQYwgyR24HZZvYgkcNqE4L2PkTuW95gU9DWVPvnmNnNwM0A/fr1a96q\nReSwyn01fLypISx2sXjTbiqqDgKQlmIM7ZXF1FG9GJ3fldH5XRjSM4t2qcc5ip6VFRmF1Wh0Vpc0\n+Pm7/8Pf3XIlP1hrfPW387l8dG9+NG04uVkdWuiTSlPiGiJmNgfodYy37gYuAu5w9+fM7BrgCWDS\nMZY9ae7+GPAYRA5nNcc2Rdq6qupDLN2853BoLN60i007DwBgBgN7ZHBuYY9ID6NvV4b37kyHdqkn\nv6NzzomMwnr66cg5kMJCuPZaSjIz+UttHb9+Yx2/mruGt1ZVcNfU07i2uO+Rnoy0uDDPiewGurq7\nm5kBu929s5k9Crzh7n8IlltJ5FDWRGCiu38jaP/Mck3RORGRk1d9qI5lW/bw8cYjgbFu+z4afl3k\nd+tIUdC7GJ3flZF9OpPVoV2L1be2Yi8/eH4JH6yv5MwB3fj5FaMozD3+xYhycqI9JxJmiCwHbnH3\nN8zsIuB+dz/DzC4DvsWRE+sPu3tJcGJ9IdAwWutDIifWK4+3H4WIyPEdqqtn1baqI4ekNu5m1bYq\nausjvxtystpTFITF6PwujOrThe6Z7UOuOnL+5Y8LNvHTWcvZX1PLLRMLuXXioNh6P/I5yXBi/evA\nQ2aWBlQTnMMAZhEJkDXAfuCrAO5eaWY/ARomwvnxiQJERD6vvt75aONOZpduY8GGSkrL9nCwth6A\nLh3bMTq/C98YNpBRfbpS1LcLvTp3IHKwILGYGdec2ZcLT8vlvpeW8fDrq3lpcRk//btRjB/UPezy\n2ozQeiItRT0REaitq2fe+kpeXrqV2aVbKa86SHpqCkV9j/QwivK70r97p4QMjGi8uaqCH/5pCRsr\nD3BNcT4/mHoaXTulh11W0kr4w1ktRSEibdXB2jreXbODl5du4bVl29i5/xAd26UycWgOk0f24sJh\nuS16HqMlHKip46HXV/Obt9fRtWM77rl8ODPG5CVtMIZJIRJQiEhbcqCmjjdXlfPy0q38dXk5VQdr\nyWqfxkWn5TJ5ZG/OH5JDx/TWf85gWdke7nphCYs37uLcwT346RdG0a97p7DLSioKkYBCRFq7PdWH\nmLuinJeXbOWNVeVUH6qnW6d2XDK8F5NH9WLCoO60T2v9wXG0unrnd+9/wgOzV1JbX893LhrC184t\nOP51KnKYQiSgEJHWqHJfDXOWbePlpVv425od1NTVk5vVnskjezF5ZC9KBmSTpl+WAGzZfYB/nVnK\n7NJtDOuVxc+vGMXYft3CLivhKUQCChFpLbbtqebV0q28vHQrH6yvpK7eye/WkSkjezF5ZG/G9u2q\ni+6OY3bpVv7lz6Vsq6rm78/qzz9eOrTVnRNqTskwxFdETmBj5X5mB8Hx4ac7cYdBORnccv4gJo/s\nxYi8zjppHKVLR0QO7T04eyVPvf8Js0u3ce+MEVw64liTaki01BMRSTBryvfyytItvFK6laWb9wAw\nvHdnpozsxZRRvXRldjP46NOd3PX8ElZsreKS4T25d8YIenfpGHZZCUWHswIKEUl07s6yLXuYvTTS\n41hdvheAsf26Rg5VjeitkUVxcKiunifeWc//nbOKtJQU/vGSIXxl/ABSdUgQUIgcphCRRFRf7yza\ntOtwcHxauZ8Ug5KCbKaM7M2lI3rRq4tmqG0Jn+7Yz91/WsLbq7dT1Lcrv7hiFKf17nziFVs5hUhA\nISKJoq7embe+ktmlW3ll6Va27qmmXaoxYVAPpozsxcXDeybEnFRtkbszc3EZP35xGbsOHOJr5xZw\n+0VD2sQ1NU1RiAQUIhKmmtp63l27ndmlW3m1dBs79tXQPi2l0VXjPenSUSOEEsWu/TX8bNZynlmw\nib7ZHbnvC6M4f0hO2GWFQiESUIhIS6s+VMebqyp4ZelW5izfRlV1LRnpqVx4Wk+mjOzFxKE5dErX\nwMhE9v66HfzghSWsq9jHjDF53HP5cHq0sV6iQiSgEJGWsPdgLX9dUc4rS7cwd0UFBw7V0aVjOy4e\nHgmOswt7aIryJHOwto5H5q7lkTfW0Ck9jbunnsbVxfltZki1QiSgEJF42bW/hteWbeOVpVt5e812\namrr6ZHZnktH9GTKyN6MG5itKTZagTXlVfzg+aXM21DJuIJsfnbFKAblZIZdVtwpRAIKEWlO5VXV\nvFq6jdmlW3lv7Q5q650+XTty6YjINRyn9+umIaKtUH298/SCjfx81nKqD9Vz2wWFfHPiwFY9J5lC\nJKAQkVPh7qwp38vcleXMWVbO/E8qcYeCHhlMHtmLKSN7MapPlzZziKOtK6+q5icvLefFxWUMzs3k\nmW+Mp1tG67xniaY9EYnRvoO1vLt2B2+sLOeNlRVs3nUAgGG9svjORYOZMrI3Q3pmKjjaoNysDvzn\ndWO5bFRvvvm7hfxp0Wa+enZB2GWFSiEibZ67s277PuauKOfNVRV8sK6Smrp6MtJTObuwB9+6sJDz\nh+SQ11XTYkjE5JG9GNYri5mLyxQiYRcgEoYDNXW8v24Hc1eWM3dlORsrI72NwbmZXD+hPxcMzaV4\nQDbpaToxLsc2fUwe97+yko2V++mb3XanpVGISJuxYfs+5gaHqN5ft4ODtfV0bJfK2YXdufm8QUwc\nktOmfxnIyZk2OhIiMxeXcdsFhWGXE5pQQsTMxgC/BjoAtcCt7j7PzL4M/DNgQBVwi7svDtbZELTV\nAbXRnPCRtq36UB0frK88fJhq/fZ9AAzskcGXx/Vn4tAcSgqydf2GxKRvdidO79eVFxUiobgfuNfd\nXzazqcHricB64Hx332lmU4DHgHGN1rvA3be3eLWSNDZW7ueNleXMXVnBu2u3U32onvZpKYwf1J0b\nJgxg4tAc+nfPCLtMaSWmF+Xxry8uY9W2Kob0bJtT9IcVIg40TJPZBSgDcPd3Gy3zPpDfwnVJkjlY\nW8f89TuDw1TlrK2I9Db6ZXfi2uK+TByWy/iB3dXbkLi4bHQeP35pGTMXlfGPlw4Nu5xQhBUitwOz\nzexBIAWYcIxlbgJebvTagVfNzIFH3f2xpjZuZjcDNwP069ev2YqWxLB514FIb2NFpLexv6aO9NQU\nxg3M5kvj+nPB0BwKemRoCK7EXU5WeyYM6sHMxWV875IhbfI7F7cQMbM5wLHuO3k3cBFwh7s/Z2bX\nAE8AkxqtewGREDmn0XrnuPtmM8sFXjOzFe7+1rH2HQTMYxC52LBZPpCEpqa2ngWfVPLmygrmrixn\n1bbITZv6dO3IFaf34YKhuYwf1F2TGkoophfl8U/PfcziTbsZ07dr2OW0uLj9q3P3SU29Z2ZPAd8J\nXv4ReLzRe6OD11PcfUej7W0Ofpab2QtACXDMEJHkt3V39eGL/d5Zs529B2tpl2qcOSCbq8/oywXD\nchiUowv+JHyXjuzFD/+0lJmLyhQiLagMOB94A7gQWA1gZv2A54GvuPuqhoXNLANIcfeq4PklwI9b\numiJn5raej76dCdvrKpg7opyVmytAqB3lw5MK+rNxKG5nF3Yg8z26m1IYunSsR3nD83hpY/LuPuy\n09rc3Glh/Yv8OvCQmaUB1QTnL4AfAd2BR4L/YTYM5e0JvBC0pQG/d/dXWrxqaTYHaur46NOdfLC+\nknnrK/lo406qD9WTmmIU9+/GnVOGMXFoDkN7Zqm3IQlvelEery3bxgfrdzBhUI+wy2lRoYSIu78D\nnHGM9q8BXztG+zqgqAVKkzjZU32IhRsaQmMHSzbv5lCdYwbDe3fmi2f246yB2Uwo7EHnDrrTnySX\nSaf1pFN6Ki8uLlOIiDSHHXsPMn9D5eGexvIte6h3SEsxRud34aZzBjKuIJvT+3fT7WEl6XVMT+Xi\n4T2ZtWQr904f2aamy1GISLMo23XgM6Gxpjwygqp9Wgqn9+vGty8czLiCbMb066pRVNIqTS/K48+L\nynh7dQUXndYz7HJajP41y0lzdzbs2M+89TsOh8amnZEJDLPap1E8oBtXnN6HcQXZjOrTtU39r0za\nrnMH59ClYztmLi5TiIg0Vl/vrCqvYt76Iz2NiqqDAGRnpFMyIJsbzy6gpCCb03p3bnOjU0QA0tNS\nmDqqF39eVMaBmjo6preNWRIUIvI5tXX1lJbtORwa8zdUsvvAIQB6de7AhEHdKSnIZlxBtq7VEGlk\nWlEef5i3kddXbOPy0Xlhl9MiFCJC9aE6Fm/cdficxsJPdrK/pg6AAd07cemInpQUdGdcQTb53Toq\nNESaMK6gO7lZ7Zm5qEwhIq3X3oO1fPjJTuYFh6YWbdxFTV09ELkF7FVn5HPmgGxKCrLp2blDyNWK\nJI/UFOPy0Xn87v1P2H3gUJsYeagQaQN27a9h/oadzFu/g3nrK1latoe6eic1xRiZ15nrJ/SnpKA7\nxf270S0jPexyRZLatKLe/Pff1jO7dCvXFPcNu5y4U4i0QuV7qpm3ofJwT6NhCpH01BTG9O3KLecP\noiS4RkPTiIg0rzF9u9IvuxMvLi5TiEjic3c27TxwODDmbag8fAe/TumpnNG/G5eN6k1JQTZFfbvq\nvhoicWZmTCvqzX+9sZaKqoPkZLUPu6S4UogkGXdnbcXew0Nt562vZMvuaiAyEdyZA7L5Ukk/Sgqy\nGZ7XmXapukZDpKVNL+rDr+auZdaSLVw/YUDY5cSVQiTB1dU7y7fsORwY8zdUsmNfDRC5IU7DUNuS\ngmyG5GaRoms0REI3tFcWQ3tmMXNxmUJEWlZNbT1LNu8OQmMHCzbspOpgLQD53Tpy/tCcIDS6M6B7\nJw23FUlQ08fk8cDslWzauZ/8bp3CLiduFCIhO1BTx0cbjwy3/fDTyJToAIW5mUwbk0fJgGzOLMim\nT9eOIVcrItGaNjoSIi8u3sItEweFXU7cKERaWOMp0edvqOTjTbs+MyX6dSX9GFeQTfGAbHpktu4T\nciKtWb/unRjTtyszF5cpRCR2mhJdpO2aXpTHj19axpryKgpzs8IuJy4UIs1sy+4Dn5mosKkp0cf2\n69ZmJmgTaasuH92b+/6yjJmLyvjuJUPDLicuFCKnwN35ZMf+I6GxYQcbKz87JfqVp+dTUtBNU6KL\ntEG5nTtw1sDuzFxcxh0XD2mVA2EUIieh8ZToDY/yo6ZE/+oETYkuIkdML8rjzueXsGTzbkbndw27\nnGanEDmOE02JPl5ToovICUwZ2Zt7/ryUmYvKFCLNyczGAL8GOgC1wK3uPs/MJgJ/BtYHiz7v7j8O\n1pkMPASkAo+7+y/iVd8tv1vIW6sq2Kcp0UXkFHTp1I7zh+Tw0sdb+MHU01rdBcFh9kTuB+5195fN\nbGrwemLw3tvufnnjhc0sFfgVcDGwCZhvZjPdfVk8iuuemc6VZ+RTUpBNyYBscjUluojEaFpRHnOW\nlzNvQyVnDewedjnNKswQcaBz8LwLUHaC5UuANe6+DsDM/heYAcQlRO77wqh4bFZE2qCLh/ekY7tU\nZi4ua3UhEuZwoduBB8xsI/AgcFej98ab2WIze9nMRgRtfYCNjZbZFLR9jpndbGYLzGxBRUVFPGoX\nEYlap/Q0Jg3vyctLtnAouAFcaxHXEDGzOWa29BiPGcAtwB3u3he4A3giWO1DoL+7FwH/CfzpZPfr\n7o+5e7G7F+fk5DTXxxERidn0ojx27j/EO2u2h11Ks4rr4Sx3n9TUe2b2FPCd4OUfgceDdfY0Wn+W\nmT1iZj2AzUDjO7zkB20iIgnvvCE96NwhjRcXlXHB0Nywy2k2YR7OKgPOD55fCKwGMLNeFgx7MrMS\nIjXuAOYDg82swMzSgS8CM1u8ahGRGLRPS2XKyN7MLt1K9aG6sMtpNscNETNLNbMVcdr314F/N7PF\nwM+Am4P2q4ClQfvDwBc9ohb4FjAbWA484+6lcapNRKTZTR+Tx76aOv66ojzsUprNcQ9nuXudma00\ns37u/mlz7tjd3wHOOEb7L4FfNrHOLGBWc9YhItJSzhrYnR6Z7Zm5qIypo3qHXU6ziOacSDeg1Mzm\nAfsaGt19etyqEhFphVJTjMtH9+b38z5lT/UhOndI/pm7owmRe+JehYhIGzGtKI8n393Aq6XbuOqM\n/LDLOWUnPLHu7m8CK4Cs4LE8aBMRkZN0er+u5HfryMzFJ7q+OjmcMETM7BpgHnA1cA3wgZldFe/C\nRERaIzNjWlEef1uznR17D4ZdzimLZojv3cCZ7n69u/89kelHdIhLRCRG04vyqKt3Zi3ZEnYppyya\nEElx98bj0XZEuZ6IiBzDsF5ZDM7NbBWHtKIJg1fMbLaZ3WBmNwB/QcNsRURiZmZML8pj/oadbN51\nIOxyTkk0J9a/DzwGjA4ej7n7P8e7MBGR1mxaUR4ALyV5bySqubPc/TnguTjXIiLSZgzokUFRfhdm\nLi7jG+cPCrucmDXZEzGzd4KfVWa2p9Gjysz2NLWeiIhEZ1pRHqVle1hbsTfsUmLWZIi4+znBzyx3\n79zokeXunZtaT0REojOtKA8zmLkoeQ9phTkBo4hIm9azcwfGFWTz4uIy3D3scmJy3BBx9zpgpZn1\na6F6RETalOlFfVi3fR+lZcl5liCaIb4NEzC+bmYzGx7xLkxEpC2YMrIXaSnGi0k6SksTMIqIhKhb\nRjrnDcnhxcVl/PPkYaSkWNglnZRoJ2DcALQLns8nch90ERFpBtOL8ijbXc3CT3eGXcpJi2YCxq8D\nzwKPBk19gD/FsygRkbbk4uE96dAuJSlHaUVzTuQ24GxgD4C7rwZaz13mRURCltE+jYtO68msJVuo\nrasPu5yTEk2IHHT3moYXZpYGJOdYNBGRBDW9KI8d+2r429odYZdyUqIJkTfN7AdARzO7GPgj8GJ8\nyxIRaVsmDs0hq0Na0h3SiiZE7gQqgCXAN4BZ7n73qezUzMaY2ftmtsjMFphZSdD+/aBtkZktNbM6\nM8sO3ttgZksa1jmV/YuIJJr2aalMHtGLV0u3Un2oLuxyohZNiHzb3X/j7le7+1Xu/hsz+84p7vd+\n4F53HwP8KHiNuz/g7mOC9ruAN929stF6FwTvF5/i/kVEEs70MXlUHazljZXlJ144QUQTItcfo+2G\nU9yvAw3zb3UBjtV/uw74wynuR0QkaYwf2J0emelJdbOqJi82NLPrgC8BBUddod4ZqDz2WlG7HZht\nZg8SCbIJR+27EzAZ+FajZgdeNTMHHnX3x45T+83AzQD9+mnGFhFJDmmpKUwd1Zun52+kqvoQWR3a\nhV3SCR3vivV3gS1AD+DfG7VXAR+faMNmNgfodYy37gYuAu5w9+fM7BrgCWBSo2WmAX876lDWOe6+\n2cxygdfMbIW7v3WsfQcB8xhAcXGxRpKJSNKYXpTHU+99wmvLtnHF6flhl3NCTYaIu38CfGJmk4AD\n7l5vZkOAYUROsh+Xu09q6j0zewpoOK/yR+Dxoxb5IkcdynL3zcHPcjN7ASgBjhkiIiLJ6vR+3ejS\nsR0LP9mZFCESzTmRt4AOZtYHeBX4CvDkKe63DDg/eH4hsLrhDTPrErz350ZtGWaW1fAcuARYeoo1\niIgknJQUozA3kzXlyXGjqmgmYDR3329mNwGPuPv9ZrboFPf7deCh4MLFaoLzF4G/A151932N2noC\nL5hZQ82/d/dXTrEGEZGEVJiTyZzl28IuIypRhYiZjQe+DNwUtKWeyk7d/R3gjCbee5Kjejruvg4o\nOpV9iogki8E9M3l6wUYq99WQnZEedjnHFc3hrNuJXLPxgruXmtlAYG58yxIRabsG5WYCJMUhrRP2\nRILp3980s0wzywx6Bf8Q/9JERNqmwY1CpKQgO+Rqji+aqeBHmdlHQCmwzMwWmtmI+JcmItI25XXp\nSMd2qUnRE4nmcNajwHfdvb+79wO+B/wmvmWJiLRdDSO0VpdXhV3KCUUTIhnufvgciLu/AWTErSIR\nEaEwN5O1raQnss7M7jGzAcHjh8C6eBcmItKWFeZmUra7mr0Ha8Mu5biiCZEbgRzg+eCRE7SJiEic\nFAYn1xO9NxLN6KydwD8EV5LXu3viH6QTEUlyDSGyunwvRX27hlxN06IZnXWmmS0BFgNLzGyxmR3z\nQkEREWke/bM70S7VEn6EVjRXrD8B3OrubwOY2TnAb4HR8SxMRKQtS0tNoaBHRsKHSDTnROoaAgQO\nT1mS2Gd6RERagcG5WaxJ8GG+0YTIm2b2qJlNNLPzzewR4A0zO93MTo93gSIibdWg3Ew+rdyf0Pdc\nj+ZwVsPEh/9yVPtYIncbvLBZKxIRESAy/Um9w/rt+zitd+cTrxCCaEZnXdAShYiIyGcVNppDK1FD\nJJrDWSIiEoKCHhmkWGSYb6JSiIiIJKgO7VLpl90poS84VIiIiCSwRJ+I8YTnRMzsimM07waWuHt5\n85ckIiINCnOzeHNVBbV19aSlJt7/+6MZnXUTMJ4jdzOcCCwECszsx+7+/+JUm4hIm1eYm8mhOufT\nyv0MzMkMu5zPiSbW0oDT3P1Kd78SGE5kaO844J/jWZyISFs3uNEcWokomhDp6+7bGr0uD9oqgUOx\n7tjMiszsPTNbYmYvmlnnRu8qyJAzAAAN6UlEQVTdZWZrzGylmV3aqH1y0LbGzO6Mdd8iIski0e+3\nHk2IvGFmL5nZ9WZ2PTAzaMsAdp3Cvh8H7nT3UcALwPcBzGw48EVgBDAZeMTMUs0sFfgVMIVIb+i6\nYFkRkVYrs30aeV06JHWI3EZkwsUxweN/gNvcfd8pXog4BHgreP4acGXwfAbwv+5+0N3XA2uAkuCx\nxt3XuXsN8L/BsiIirdqg3MzkDRF3d+Ad4K/A68BbQdupKuVICFwN9A2e9wE2NlpuU9DWVPvnmNnN\nZrbAzBZUVFQ0Q6kiIuEpDEKkvr45fvU2r2juJ3INMA+4CrgG+MDMropm42Y2x8yWHuMxg8jdEW81\ns4VAFlAT+8f4LHd/zN2L3b04JyenuTYrIhKKwblZHDhUR9nuA2GX8jnRDPG9Gziz4ZoQM8sB5gDP\nnmhFd590gkUuCbY5BLgsaNvMkV4JQH7QxnHaRURarcZ3Oczv1inkaj4rmnMiKUddVLgjyvWOy8xy\ng58pwA+BXwdvzQS+aGbtzawAGEykJzQfGGxmBWaWTuTk+8xTrUNEJNENTuD7rUfTE3nFzGYDfwhe\nXwvMaoZ9X2dmtwXPnydy8h53LzWzZ4BlRG5+dZu71wGY2beA2UAq8N/uXtoMdYiIJLRuGel0z0hP\nyJPrFs05cjO7Ejg7ePm2u78Q16qaUXFxsS9YsCDsMkRETsm1j75Hbb3z3C0TWmR/ZrbQ3YtPtFw0\nPRHc/TnguVOuSkREYlKYm8lLH2/B3TGzsMs5rMkQMbMqItObfO4tIiN/E/MOKSIirdDg3Ex2HzhE\nxd6D5GZ1CLucw5oMEXfPaslCRESkaYW5kV/Ja8r3JlSIJN68wiIi8jmFCTqHlkJERCQJ9Ozcnqz2\naQoRERE5eWbGoNxMVm9TiIiISAwG52aypkIhIiIiMSjMzaSi6iC798d8K6dmpxAREUkSg3sGJ9cr\nqkKu5AiFiIhIkijMOTLMN1EoREREkkSfbh1pn5aSUCfXFSIiIkkiNcUYlJNYJ9cVIiIiSaQwwYb5\nKkRERJLI4NxMNu86wP6a2rBLARQiIiJJpfDwDar2hVxJhEJERCSJJNowX4WIiEgS6d89g7QUS5hh\nvgoREZEk0i41hQE9MhLm5LpCREQkyRQm0DBfhYiISJIpzM3kkx37qamtD7uUcELEzIrM7D0zW2Jm\nL5pZ56D9YjNbGLQvNLMLG63zhpmtNLNFwSM3jNpFRMI2uGcmdfXOhh3hj9AKqyfyOHCnu48CXgC+\nH7RvB6YF7dcD/++o9b7s7mOCR3nLlSsikjgG5URGaCXCeZGwQmQI8Fbw/DXgSgB3/8jdy4L2UqCj\nmbUPoT4RkYQ1KCcTs8SYiDGsECkFZgTPrwb6HmOZK4EP3f1go7bfBoey7jEza2rjZnazmS0wswUV\nFRXNV7WISALomJ5KfreOrC4P/1qRuIWImc0xs6XHeMwAbgRuNbOFQBZQc9S6I4B/A77RqPnLwWGu\nc4PHV5rat7s/5u7F7l6ck5PT3B9NRCR0g3OzEqInkhavDbv7pBMscgmAmQ0BLmtoNLN8IudJ/t7d\n1zba3ubgZ5WZ/R4oAZ5q7rpFRJJBYW4m76zZTl29k5rS5IGZuAtrdFZu8DMF+CHw6+B1V+AvRE66\n/63R8mlm1iN43g64HFja0nWLiCSKwtxMamrr2Vi5P9Q6wjoncp2ZrQJWAGXAb4P2bwGFwI+OGsrb\nHphtZh8Di4DNwG9CqFtEJCE0TMQY9iGtuB3OOh53fwh46Bjt9wH3NbHaGXEtSkQkiTSEyOryvUwa\n3jO0OnTFuohIEurcoR09O7cPvSeiEBERSVKFuZmsCXmYr0JERCRJNQzzdffQalCIiIgkqUG5meyr\nqWPL7urQalCIiIgkqcEJMEJLISIikqQSYZivQkREJEl1z0ina6d2rFaIiIjIyTIzBudmslYhIiIi\nsSjMzQx1Nl+FiIhIEivMzWLn/kPs2HvwxAvHgUJERCSJNZ7+JAwKERGRJBb2MF+FiIhIEuvdpQMZ\n6akKEREROXlmFsyhpRAREZEYDFKIiIhIrApzM9m6p5o91YdafN8KERGRJDc4NwsglIsOFSIiIkku\nzGG+ChERkSTXt1tH0tNS2k5PxMyKzOw9M1tiZi+aWeegfYCZHTCzRcHj143WOSNYfo2ZPWxmFkbt\nIiKJJi01hYE9MtpUT+Rx4E53HwW8AHy/0Xtr3X1M8Phmo/b/Ar4ODA4ek1usWhGRBBfWMN+wQmQI\n8Fbw/DXgyuMtbGa9gc7u/r5H7gP5FPCF+JYoIpI8CnMz2bhzP9WH6lp0v2GFSCkwI3h+NdC30XsF\nZvaRmb1pZucGbX2ATY2W2RS0iYgIkRFa7rC2omV7I3ELETObY2ZLj/GYAdwI3GpmC4EsoCZYbQvQ\nz93HAt8Fft9wvuQk932zmS0wswUVFRXN9ZFERBJWWHc5TIvXht190gkWuQTAzIYAlwXrHAQOBs8X\nmtlaIoe+NgP5jdbND9qa2vdjwGMAxcXFHuNHEBFJGgN6dCLFWj5EwhqdlRv8TAF+CPw6eJ1jZqnB\n84FETqCvc/ctwB4zOysYlfX3wJ/DqF1EJBG1T0tlQPeMthEiwHVmtgpYAZQBvw3azwM+NrNFwLPA\nN929MnjvViKjutYAa4GXW7ZkEZHENig3s8WH+cbtcNbxuPtDwEPHaH8OeK6JdRYAI+NcmohI0hqc\nm8ncFeUcqqunXWrL9BF0xbqISCtRmJtJbb3zyY59LbZPhYiISCvRMBFjS54XUYiIiLQSg3IzAIWI\niIjEoFN6Gn26dmzRk+sKERGRVqSl59BSiIiItCKFuZmsrdhLfX3LXGetEBERaUUG52ZSfaiezbsO\ntMj+FCIiIq3IkbscVrXI/hQiIiKtSEtPxKgQERFpRbp2SqdHZntWb1OIiIhIDAbnZrKmhe4rohAR\nEWllGob5Rm4EG18KERGRVmZY7yxys9qzp7o27vsKZRZfERGJny+P68+Xx/VvkX2pJyIiIjFTiIiI\nSMwUIiIiEjOFiIiIxEwhIiIiMVOIiIhIzBQiIiISM4WIiIjEzFrisvgwmVkF8AnQBdjdxGJNvdcD\n2B6n0prD8T5TImz7ZLdxMstHs2wsf+fHe0/fh5bfRrTrnOr34Xjvt9XvQ393zznhUu7eJh7AYyf7\nHrAg7Lpj/UyJsO2T3cbJLB/NsrH8nev7kDjfh5NZ51S/Dyf4e9f34TiPtnQ468UY30tk8ay7ObZ9\nsts4meWjWTbWv3N9H+Kz7Vi2Ee06p/p9ON77+j4cR6s/nHUqzGyBuxeHXYckBn0fpDF9HyLaUk8k\nFo+FXYAkFH0fpDF9H1BPREREToF6IiIiEjOFiIiIxEwhIiIiMVOIxMjMvmBmvzGzp83skrDrkXCZ\n2UAze8LMng27FgmHmWWY2f8Evxe+HHY9LaVNhoiZ/beZlZvZ0qPaJ5vZSjNbY2Z3Hm8b7v4nd/86\n8E3g2njWK/HVTN+Hde5+U3wrlZZ2kt+NK4Bng98L01u82JC0yRABngQmN24ws1TgV8AUYDhwnZkN\nN7NRZvbSUY/cRqv+MFhPkteTNN/3QVqXJ4nyuwHkAxuDxepasMZQpYVdQBjc/S0zG3BUcwmwxt3X\nAZjZ/wIz3P3nwOVHb8PMDPgF8LK7fxjfiiWemuP7IK3TyXw3gE1EgmQRbeg/6G3mg0ahD0f+FwGR\nL0Sf4yz/bWAScJWZfTOehUkoTur7YGbdzezXwFgzuyvexUmomvpuPA9caWb/RfJOlXLS2mRPpDm4\n+8PAw2HXIYnB3XcQOT8mbZS77wO+GnYdLU09kSM2A30bvc4P2qRt0vdBmqLvRiMKkSPmA4PNrMDM\n0oEvAjNDrknCo++DNEXfjUbaZIiY2R+A94ChZrbJzG5y91rgW8BsYDnwjLuXhlmntAx9H6Qp+m6c\nmCZgFBGRmLXJnoiIiDQPhYiIiMRMISIiIjFTiIiISMwUIiIiEjOFiIiIxEwhIhISMxtgZl9q9PoG\nM/tlE8vOMrOuLVedSHQUIiLhGQB86UQLAbj7VHffFd9yRE6eQkTkKMEd6v5iZovNbKmZXWtmG8zs\n52a2yMwWmNnpZjbbzNY2zOJsEQ8E6ywxs2uP107kVgLnBtu8I2jLM7NXzGy1md3fqKYNZtYj6L0s\nD+6eV2pmr5pZx2CZM83s42B7Dxx9IyWReFCIiHzeZKDM3YvcfSTwStD+qbuPAd4mcrOiq4CzgHuD\n968AxgBFRG4T8ICZ9T5O+53A2+4+xt3/I9jGGCJ3yhwFXGtmjSf6azAY+JW7jwB2AVcG7b8FvhHU\n2GZuiiThUoiIfN4S4GIz+zczO9fddwftMxu9/4G7V7l7BXAwOF9xDvAHd69z923Am8CZx2k/ltfd\nfbe7VwPLgP7HWGa9uy8Kni8EBgT7z3L394L238f86UVOgu4nInIUd19lZqcDU4H7zOz14K2Dwc/6\nRs8bXjfXv6XG261rYrtHL9OxmfYtctLUExE5ipnlAfvd/XfAA8DpUa76NpFDUKlmlgOcB8w7TnsV\nkNUcNQcn3avMbFzQ9MXm2K7IiagnIvJ5o4ict6gHDgG3AM9Gsd4LwHhgMeDAP7n7VjNrqn0HUGdm\ni4mcY9l5inXfBPwmqPtNYPcJlhc5ZZoKXqSVMLNMd98bPL8T6O3u3wm5LGnl1BMRaT0uM7O7iPy7\n/gS4IdxypC1QT0RERGKmE+siIhIzhYiIiMRMISIiIjFTiIiISMwUIiIiEjOFiIiIxOz/AwS0C+qj\nR9+RAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f09eee30630>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "smoothers = [0.01,0.05,0.1,0.25,0.5,1,1.5,2]\n",
    "log_posteriors = [log_posterior((x_tr + smoother)/(x_tr.sum() + V * smoother)) for smoother in smoothers]\n",
    "plt.semilogx(smoothers,log_posteriors)\n",
    "plt.scatter(a-1,log_posterior((x_tr + a-1)/(x_tr.sum()+V*(a-1))),50,color='r')\n",
    "plt.text(a-1+.1,log_posterior((x_tr + a-1)/(x_tr.sum()+V*(a-1))),\"$\\\\alpha = a-1$\")\n",
    "plt.ylabel(\"log posterior\")\n",
    "plt.xlabel(\"smoothing\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As predicted, setting the smoothing equal to $a - 1$ maximizes the posterior log-probability. \n",
    "\n",
    "In real situations, the parameter $a$ is not known, but [it can be estimated from data](https://tminka.github.io/papers/dirichlet/minka-dirichlet.pdf)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## All your Bayes: from smoothing to marginalization\n",
    "\n",
    "Some people feel that smoothing has an unsatisfying \"two step\" feel to it:\n",
    "\n",
    "- Use the training data and the smoothing hyperparameter to estimate $\\hat{\\theta}$\n",
    "- Use the estimate $\\hat{\\theta}$ to compute the probability of the observed data $P(\\vec{x} \\mid y)$\n",
    "\n",
    "Rather than using the *point estimate* $\\hat{\\theta}$, it would be nicer to consider all possible values for $\\theta$. This is the [Bayesian](https://en.wikipedia.org/wiki/Bayesian_inference) approach. \n",
    "\n",
    "Specifically, we would like to compute the *posterior predictive distribution*,\n",
    "\\begin{align}\n",
    "P(\\vec{x} \\mid y, \\vec{x}^{(TR)}; a) = &{} \\int P(\\vec{x} | y, \\vec{\\theta}_{y}) P(\\vec{\\theta}_{y} \\mid \\vec{x}^{(TR)}; a) d\\vec{\\theta}_{y}\n",
    "\\end{align}\n",
    "If we could solve this integral, we would never need to choose any specific value for $\\vec{\\theta}_{y}$. Instead, we marginalize over all possible values, with each value counting in proportion to its posterior probability $P(\\vec{\\theta}_{y} \\mid x^{(TR)}; a)$. \n",
    "\n",
    "Integrating over the vector $\\vec{\\theta}$ may not seem promising, but miraculously, this integral has a closed form solution! \n",
    "\n",
    "\\begin{align}\n",
    "P(\\vec{x} \\mid y, \\vec{x}^{(TR)}; a) = &{} \\int P(\\vec{x} | y, \\vec{\\theta}_{y}) P(\\vec{\\theta}_{y} \\mid \\vec{x}^{(TR)}; a) d\\vec{\\theta}_{y}\\\\\n",
    "\\propto &{} \\prod_w \\frac{\\Gamma(x_w + x^{(TR)}_w + a_w)}{x_w! \\Gamma(x^{(TR)}_w + a_w)}\n",
    "\\end{align}\n",
    "\n",
    "The solution is expressed in terms of the [**Gamma function**](https://en.wikipedia.org/wiki/Gamma_function), $\\Gamma(n)$. It is a generalization of the factorial function: for non-negative integer inputs $\\Gamma(n) = (n-1)!$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 135,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from scipy.special import gamma, factorial"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 171,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAEACAYAAAC9Gb03AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAIABJREFUeJzt3Xl8ldW97/HPLyEhhCFABiAJIUFm\nERJMcaooVSmoAeS0VUtPW/WY2tPh3N7b2+GFbbWnaXvP6eiprYdWj7avtGptQaCorbaKVqwoIYAy\nkxDIQAYChAxkWvePnWBAAjvJTp69d77v1ysvutd+9vP84HH/srqetX7LnHOIiEj4ivA6ABER6V9K\n9CIiYU6JXkQkzCnRi4iEOSV6EZEwp0QvIhLmlOhFRMKcEr2ISJhTohcRCXNK9CIiYW6I1wEAJCQk\nuPT0dK/DEBEJKW+//Xa1cy7xYsd5mujNLAfImTJlCm+99ZaXoYiIhBwzO+TPcZ4O3Tjn1jvncuPi\n4rwMQ0QkrHma6M0sx8xWnzhxwsswRETCmnr0IiJhLmjG6M/V0tLCkSNHaGpqGvjAQkBMTAypqalE\nRUV5HYqIBDkLho1HsrOz3bkPY4uKihg5ciTx8fGYmUeRBSfnHDU1NdTV1ZGRkeF1OCLiETN72zmX\nfbHjgnYefVNTk5J8N8yM+Ph4/b8dEfFLUD+MVZLvnv5tRMRfehgrIjJQ8vMhPR0iInx/5ucPyGWD\nYmWsiEjYy8+H3FxoaPC9PnTI9xpg5cp+vXRQD90Eg8jISDIzM8/8PPLIIzQ2NnLdddfR1tbW7eea\nm5tZsGABra2tAxitiAStVatoa2wif+5iisYk+9oaGmDVqn6/tIZuLmLYsGFs27btzM99993HY489\nxooVK4iMjOz2c9HR0dxwww089dRTAxitiAStkhLKR8azavHneSPtsrPa+1vQzroJZvn5+SxbtuzM\n64ULF/KXv/wFgPvvv58vfOELACxfvpz8ARqDE5Egl5ZGcUdPflJt2Vnt/S0kxugfXP8O75adDOg5\nZyWP4ls5l/b4c83NzRw8eJCu1TYffPBBvvnNb1JZWUlBQQHr1q0DYPbs2WzZsiVQIYtIKMvLo+ih\n3wOQ0ZnoY2MhL6/fLx20K2ODVXV1NaNHjz6rbcGCBTjn+NGPfsTLL798ZkgnMjKS6Oho6urqGDly\npBfhikiwWLmSQ+UxxFQ0M+5ULUya5Evy/fwgFjxO9M659cD67Ozsey90XG963v1l2LBh71uotGPH\nDsrLy4mPj39fQj99+jQxMTEDGaKIBKnixDTShzQS0d79RI7+oDH6HhozZgxtbW1nkn15eTkrV67k\n2WefZcSIETz//PNnjq2pqSEhIUH1aEQEgKLqeibFxw74dZXoe2HRokW89tprNDQ0sGLFCn74wx8y\nc+ZMvvGNb/Dggw+eOe5vf/sbt9xyi4eRikiwaGt3HD7WSHrC8AG/thL9RZw6dep9bZ/73Od44okn\niI2NZfPmzdx0002Ab6x+8+bNZ4777W9/y2c+85kBi1VEglfZ8Uaa29rJiFeiDwnz5s1j4cKFF10w\ntXz5cqZNmzaAkYlIsCquqQdg0mBL9KGwMrY7d99990UXTH3yk58cwIhEJJgV1/hKH2QMtqGbUFgZ\nKyISCMXV9cRERTBu1NABv7aGbkREBkBxdT3p8cM9KTGuRC8iMgCKa3yJ3gtK9CIi/czLqZWgRC8i\n0u86p1ame7BYCpToRUT6XefUyrDp0ZvZ9Wb2qpk9YmbXB/r8IiKhprjal+i9mFoJfiZ6M3vMzCrN\nbOc57YvNbI+Z7Tezr3U0O+AUEAMcCWy4A8+f3aR64+677yYpKYnZs2efaTvfrlQPPPAAP/jBDwJ6\nbREZWEXVDQyLiiRp5MBPrQT/e/SPA4u7NphZJPAwsASYBdxpZrOAV51zS4CvAg8yUPpp011/dpPq\njU9/+tNnFUAD7UolEq4O1fiKmXkxtRL8TPTOuU3AsXOa5wP7nXMHnXPNwJPAMudce8f7tcDA/Prq\n3HT30CFw7r1NdwOQ7LvuJtXdTlK9sWDBAsaOHfu+du1KJRJ+imrqPRu2gb7Vo08BDnd5fQS4wsxW\nAB8GRgM/6+7DZpYL5AKk9XUrrVWr3ttZvVPnprt9KOp/7m5S3e0kFUjalUokvLS2tXP4WAOLZo33\nLIaAbzzinPsj8Ec/jlttZuVATnR09OV9umh3m+v2cdPdc3eT6m4nqU433ngjFRUV7ztPXl7eWXvM\nXoh2pRIJL+Unmmhpc2QkeDO1EvqW6EuBiV1ep3a0+c3fHaYuKi3NN1xzvvY+OHc3qQvtJAXw4osv\n9ul6nbruSvXAAw8E5Jwi4o2ijhk3Xq2Khb5Nr9wCTDWzDDOLBu4AejSWEbDqlXl5vk12uwrAprtd\nd5O60E5SgaRdqUTCi9dz6MH/6ZW/AzYD083siJnd45xrBT4PvADsAp52zr3Tk4sHrHrlypWwerVv\ns10z35+rVwdk091Fixbx5z//+YI7SfXGnXfeyVVXXcWePXtITU3l0UcfBd6/K9UjjzzCr3/96z5d\nS0S8U1zdQGy0d1MrAcw5593FzXKAnClTpty7b9++s97btWsXM2fO9CawLrZu3cqPf/xjfvOb3wzI\n9VasWMH3v/99vzYsCZZ/IxHp3t2Pb6H8RBPP/du1AT+3mb3tnMu+2HGqR38R/uwmFSjalUok/PjK\nE3v3IBa0w5RfLrabVKBoVyqR8NLa1s7h2gZPx+dBPXoRkX5TdrxjaqWHM25A1StFRPpN0ZkNwTV0\n0+3QjZcPioOd/m1Eglx+Pof+9f8AkHHzwoDV3+qNoB26iYmJoaamRgntPJxz1NTUnFlUJSJBpqP+\nVhExxDY3krhnR8Dqb/VGwEsgBEpqaipHjhyhqqrK61CCUkxMDKmpqV6HISLn01F/q3hMMpNqyzEI\nSP2t3vI00XeZR/++96KiosjIyBj4oERE+qqjztah0ROYUVX8vvaBFrRDNyIiISstjVaLoGT0eNJr\ny85q94Jm3YiIBFpeHqXjJtEaOYT02nJfWwDqb/VW0I7Ri4iErJUrKa4bAsWQfrzMV38rL8+T8XlQ\nohcR6RfFc66A4ndIf+ctGOntDLmgnkcvIhKqiqrrGR4dSeII76pWdtLDWBGRfuDbEHy4ZxuCd6WH\nsSIi/aC4psHTDcG7UqIXEQmwzg3Bva5x00mJXkQkwEqPN9La7jwvT9xJiV5EJMAOVJ0CYLISvWbd\niEh42l1RB8C08SM9jsRHs25ERAJsd3kdKaOHMSomyutQAA3diIgE3J6KOqYHSW8elOhFRAKqubWd\nA1WnmKFELyISng5UnaK13alHLyISrvZ0PIidMX6Ux5G8p18SvZkNN7O3zOzW/ji/iEiw2l1RR1Sk\nMTkxOKZWgp+J3sweM7NKM9t5TvtiM9tjZvvN7Gtd3voq8HQgAxURCQW7K05ySeIIoiKDZ8DE30ge\nBxZ3bTCzSOBhYAkwC7jTzGaZ2U3Au0BlAOMUEQkJeyrqgupBLPhZj945t8nM0s9png/sd84dBDCz\nJ4FlwAhgOL7k32hmG51z7QGLWEQkSJ1oaKH8RBMzJgTP+Dz0beORFOBwl9dHgCucc58HMLNPA9Xd\nJXkzywVyAdI82kdRRCSQdlecBAiqGTfQj7NunHOPO+c2XOD91c65bOdcdmJiYn+FISIyYPYc9c24\nmRlEM26gb4m+FJjY5XVqR5vfVOtGRMLJ7oo64oZFMW6U97tKddWXRL8FmGpmGWYWDdwBrOvJCVTr\nRkTCye7yk0wfPzIodpXqyt/plb8DNgPTzeyImd3jnGsFPg+8AOwCnnbOvdOTi6tHLyLhor3dsffo\nKWYG2fg8+D/r5s5u2jcCG3t7cefcemB9dnb2vb09h4hIMCg93sip061MD7LxeVA9ehGRgOisQR9s\nM25A9ehFRPouP589X/4mANNvuBLy8z0O6Gzq0YuI9EV+PuTmsitqNBOPVzDiwF7IzQ2qZK8evYhI\nX6xaBQ0N7ElMZ3pVsa+tocHXHiSCp+qOiEgoKimhKTKKorEpzOxM9B3twaIvJRD6zMxygJwpU6Z4\nGYaISO+lpbG/MZK2iMj3evQd7cFCQzciIn2Rl8ee1GkAzOhM9LGxkJfnXUzn8LRHLyIS8lauZE95\nDNFHW0ivLYdJk3xJfuVKryM7Q0M3IiJ9tGtMKlOjmxnS1up1KOeloRsRkT7aU1EXlAulOmnWjYhI\nHxyrb6ay7nTQlSbuSoleRKQPgnWzka60MlZEpA/2dNS4CbZ9YrvSGL2ISB/sLq9j7PBoEkcG12Yj\nXWnoRkSkD3YfrWP6uODbbKQrJXoRkV5qb3fsDfIZN6BELyLSayXHGmhsaWPmBCV6EZGw9N5mI8E7\ntRI060ZEpNd2lB5nSIQxfZx69N3SrBsRCWUFJceZOWEUw6IjvQ7lgjR0IyLSC23tjsLDx8lKG+11\nKBelRC8i0gt7j9ZR39ymRC8iEq4KSo4DMC9tjMeRXJwSvYhIL2wtqWXs8GjSxsZ6HcpFKdGLiPRC\nQUktWRNHB/WK2E4BT/RmNtPMHjGzZ8zss4E+v4iI1040tHCgqp55k4J/2Ab8TPRm9piZVZrZznPa\nF5vZHjPbb2ZfA3DO7XLO3Qd8DLgm8CGLiHir4HAtAFkTg/9BLPjfo38cWNy1wcwigYeBJcAs4E4z\nm9Xx3lLgT8DGgEUqIhIkCkqOE2EwJ5wSvXNuE3DsnOb5wH7n3EHnXDPwJLCs4/h1zrklQPDsjisi\nEiBbS2qZNm4kI4Z6uu223/oSZQpwuMvrI8AVZnY9sAIYygV69GaWC+QCpKWl9SEMEZGB097u2Hb4\nOLfOSfY6FL8F/NeRc+5l4GU/jlttZuVATnR09OWBjkNEpD8crD5FXVMr80JgoVSnvsy6KQUmdnmd\n2tHmN9W6EZGQkp/P1ttzAci66yOQn+9xQP7pS6LfAkw1swwziwbuANb15ASqXikiISM/H3JzKRia\nwKimU0zeuQVyc0Mi2fs7vfJ3wGZgupkdMbN7nHOtwOeBF4BdwNPOuXd6cnH16EUkZKxaBQ0NFCRP\nJ6tsNxE4aGjwtQc5v8bonXN3dtO+kT5MoTSzHCBnypQpvT2FiMjAKCmhLnoYexInsXjv62e1BzvV\noxcR8UdaGtsnTMNZBFlle85qD3baYUpExB95eRRMmg1AZmeij42FvDwPg/KPevQiIv5YuZKtiz7C\nlBPlxDU3wKRJsHo1rAz+daGhsaxLRMRjzjkK2mK56YYPwC/avQ6nRzR0IyLih+KaBmobWsgKgY1G\nzqWhGxERPxSUdFSsDKEVsZ208YiIiB8KSo4zYugQpiaN9DqUHtPQjYiIH7aW1DJ3YhyREcG/o9S5\nNHQjInIRDc2t7K6oI2ti6I3Pg4ZuREQuaseRE7S1O+ZNCr3xeVCiFxG5qK0lxwHIVI9eRCQ8vX6g\nmilJIxg7PNrrUHpFD2NFRC6gqaWNfxQdY8HURK9D6TU9jBURuYB/FB2jubWdBdMSvA6l1zR0IyJy\nAZv2VhE9JIIrMuK9DqXXlOhFRC7g1X1VzE8fy7DoSK9D6TUlehGRbpSfaGTv0VMhPWwDSvQiIt16\ndW81AAumhe6DWNCsGxGRbm3aV0XSyKFMHxd69W260qwbEZHzaGt3vLa/mmunJmIWevVtutLQjYjI\neewoPcHxhpaQH58HJXoRkfN6dW8VZvDBKUr0IiLhJT8f0tPZ9NgaZtccIv7ZZ7yOqM+U6EVEOuXn\nQ24uJ8ur2JoygwW734DcXF97COuXRG9my83sl2b2lJkt6o9riIgE3KpV0NDA5klzaIuI5NriAmho\n8LWHML8TvZk9ZmaVZrbznPbFZrbHzPab2dcAnHNrnXP3AvcBtwc2ZBGRflJSAsCmjHkMP93AvNLd\nZ7WHqp706B8HFndtMLNI4GFgCTALuNPMZnU55P6O90VEgl9aGg5for+qZDvR7a1n2kOZ34neObcJ\nOHZO83xgv3PuoHOuGXgSWGY+/w94zjm3NXDhnq2gpJZfvXqQ8hON/XUJERlM8vI4NCGDw6PHs6Co\nwNcWGwt5ed7G1UdD+vj5FOBwl9dHgCuALwA3AnFmNsU598i5HzSzXCAXIK2Xvy1f3lPFT1/aR97G\nXXwgfSw5c5O5efZ44kcM7dX5RGSQW7mSTTXRUAYLigtg0iRfkl+50uvI+qSvif68nHMPAQ9d5JjV\nwGqA7Oxs15vrfOmmaSzLTGbD9nLWFZbxjbU7eWDdO1wzJYGcORP48OzxjIqJ6s2pRWSQ2hSXTlrT\nKdKPlXodSsD0NdGXAhO7vE7taPOLmeUAOVOmTOl1AJMTR/DFG6byhQ9NYXdFHesKy1hfWMb/fWY7\nq9buZOH0RHLmJnPDjHEhXWZURPpfc2s7mw/UsDwrxetQAqqviX4LMNXMMvAl+DuAj/v7YefcemB9\ndnb2vX2MAzNj5oRRzJwwiq98eDoFh4+zvrCMDdvLeeGdo8RGR3LTrHEsnZvMtVMTiR6iJQQicrat\nJbXUN7eFfLXKc/md6M3sd8D1QIKZHQG+5Zx71Mw+D7wARAKPOefe6cE5+9yj7+a8zEsbw7y0Mdx/\nyyz+UVTD+sIyNu6o4NltZcQNi2LJ7PEsnZvMFZPjiYwI7YJFIhIYm/ZWERlhXH1J6O4mdT7mXK+G\nxwMqOzvbvfXWW/1+nebWdl7bX8W6bWX8+d2jNDS3kThyKLdcNoGcucnMSxsd8lXqRKR3nHMs+vEm\nxg6P5qnPXOV1OH4xs7edc9kXO65fHsb6q7969N2JHhLBh2aM40MzxtHY3MZfd1eyvrCM375ZwuOv\nF5M6Zhg5c5PJmZPMzAkjlfRFBpFd5XXsqzzFvy+f7XUoATeoevTdOdnUwp/fOcr6wjJe219NW7tj\nStIIcuYkszQzmYyE4Z7FJiID4/vP7eZXrx7kzVU3MnZ4tNfh+MXfHr2nib5Lj/7effv2eRZHVzWn\nTrNxZwXrC8t4s8i3PuyylDhy5k7g1jnJJI8e5nGEIhJo7e2Oa//jb0wdN4LH75rvdTh+C4lE38nr\nHn13yk80sqGwnPXby9h+xLfd4fz0seRkamGWSDh5q/gYH3lkMz/62FxWzEv1Ohy/KdEHWFF1PesL\ny1hXWMb+ylNERhgfnJLA0rnJLLp0HCO1MEskZH3z2Z08teUwb3/jJkYM9fTRZY+ERKIPxqGbi3HO\nnbUw60htI0OHRPChGUksnZvMwhlJxERpYZZIqGhta+eK777ElZPjeXjlPK/D6ZGQSPSdQqFHfz7O\nObaWvLcwq/rUaUYMHcKiS30Ls66ZkkBUpBZmiQSzV/ZW8anH3uSRT1zO4tnjvQ6nR0JiemWoMzMu\nnzSGyyeN4f5bZvLGwWOsKyzl+Z0V/HFrKWOHR3PzZeNZOjeF7EljiNDCLJGgs25bGSOHDuH66eG1\nGrYrDd30g9Otbbyyp4p1hWW8uOsoTS3tJMfFcOvcZJbOTebS5FGaoy8SBJpa2sj+zossmT2e//zo\nXK/D6TEN3QSJ+tOtvLjrKOu2lfHK3ipa2x2TE4eztCPpT04c4XWIIoPWczvK+Wz+Vn5zz3yunRp6\nPXoN3QSJ4UOHsCwzhWWZKRxvaOa5nRU8u62Un760j5+8uI/ZKaNYNjeFW+dOYEKc5uiLDKR1hWUk\njIjmqsnhVdvmXEr0A2h0bDR3zk/jzvlpVJxoYsN238ydvI27+O5zu5ifPpalmcncPHsCY0JkZZ5I\nqKprauGl3ZXc+YGJDAnzSRMaow8CnXP0n91WyoGqeoZEGAumJbJ0bjI3zRrH8BCa1ysSKp55+whf\n/n0hf/js1Vw+aYzX4fSKxuhDkHOOd8tP+ubobyuj7EQTw6IiubGjjv5101RHX6TP8vNh1So+Of9u\nDiZN4tUrh2CfCM2tAjVGH4LMjEuT47g0OY6vfngGbx2qZV1hKRt3+GrvnKmjn5nMFRmqoy/SY/n5\nkJtLNVH8PT2T3H/8Efuf34MR8vvCXoh69CGgpa2d1/ZXs25bGS+8U0FDcxvjRg0lZ04yyzJTmJ2i\n6ZoifklPh0OH+HXWLXxz0Wd5/tHPMaP6kG8T8OJir6PrMfXow0hUZAQLpyexcHoSjc1tvumahWU8\nsbmYX71WREaCb7rmskxN1xS5oJISANZeej3Tqg75knyX9nClRB9ihkVH+jZHmZvMiYYWnttZzrPb\nynjor/v46Uv7uCwljmWZydw6J5nxcTFehysSXNLS2H46iq0pM7n/r786qz2cadZNmOicrrmu0FdS\n2QyuzIhnWWYyS2ZPIC5W1TVFyM/ni8/s5K/p83j9559mVHMDxMbC6tUhOUavWTeD2MGqUzy7zZf0\ni6rriYo0rp+exPLMFG6YqeqaMniVHW/k2u+/xF27/8r9637i68nn5YVkkgclesE3XXNH6QnWFpSx\nYXsZlXXvVddcnpnC1ZfEh/1CEZGuvrdxF7989SCbvrKQ1DGxXofTZ3oYK5gZc1JHMyd1NKtumckb\nB2t4dlspz+3wVddMGBHNrXN8D3EzJ47WzB0Ja6dOt/LbN0tYctmEsEjyPaFEP0hERhjXTEngmikJ\nfHvZbF7eU8nagjJ++2YJj79eTHp8LEszU1iumTsSpn7/1mHqmlr5lw9meB3KgNPQzSB3sqmF5zsK\nrb1+oAbnODNzZ+ncZJJGaeaOhL62dsf1P/gbSSNj+MNnr/Y6nIDR0I34ZVRMFB/LnsjHsidy9GRT\nR82dMr7zp118d+Murr4kgWWZySyePV774krI+vM7FRw+1siqm2d6HYonAt6jN7PJwCogzjn3EX8+\nox598NlfeYp120pZu62MkmMNDB0SwY0zx7EsM5nrpyep5o6ElH/6xetU1jXx8pcXhlXpkID26M3s\nMeBWoNI5N7tL+2Lgp0Ak8Cvn3PedcweBe8zsmd6FLsFgStII/vei6XzppmlsLTnOs9tK2bC9nD/t\nKCduWBS3zJnA8kxtkSjBb2tJLW8fquVbObPCKsn3hL9DN48DPwN+3dlgZpHAw8BNwBFgi5mtc869\nG+ggxTtd98X9xq2zeG1fNWu3lbJmaym//UcJKaOHsSwzmeVZKUwbN9LrcEXe59HXihgZM4SPZk/0\nOhTP+JXonXObzCz9nOb5wP6OHjxm9iSwDFCiD1NRkREsnJHEwhlJ1J9u5c/vVrC2oIz/3nSQn798\ngFkTRrE8K5mlc1NUfkGCwuFjDTy3o5x7r53MiEG8r0Nf/uYpwOEur48AV5hZPJAHZJnZ151z3zvf\nh80sF8gFSAvzOhPhaPjQIdyWlcptWalU1Z1mw/Yy1m4r47sbd/O953Zz9SXxLMtMYYke4oqHnni9\nmAgzPn1NuteheCrgv+KcczXAfX4ct9rMyoGc6OjoywMdhwycxJFDueuaDO66JoODVadYu62MtQWl\nfOWZ7Xxj7U5unDWO2zJTWKCNU2QAVZ86zZNbDnPzZdqPuS+JvhToOuiV2tHmN+fcemB9dnb2vX2I\nQ4LI5MQR/O+bpvGlG6dScPg4a7aWsmF7GX/aXs6Y2ChuneMbz5+XppW40r9++Oc9NLW08W83TvU6\nFM/1JdFvAaaaWQa+BH8H8PGenKBL9co+hCHByMyYlzaGeWlj+GbOLDbtrWJNQSlPv3WY37xxiLSx\nsSzPSuG2rBQyEoZ7Ha6EmZ2lJ3hyy2HuuSaDS7TS27959Gb2O+B6IAE4CnzLOfeomd0M/ATf9MrH\nnHN5vQlC8+gHj7qOlbhrCkrZfNC3Ejdz4mhuy0rh1jkTiB8x1OsQJZTl5+NWreL2az7LgcRJ/HVe\nO3Gf7FH/M6SERPVK1aMf3CpONPHstlLWFJSyu6KOIRHGddMSWZ6Vwk2zxqmcsvRMx36wG9Lm8fll\nX+N7z/8Xd+57NWRrzfsjJBJ9J/XoZVf5SdYWlLJ2WylHT/rKKS+ZPZ7b5qVwZUa8FmXJxaWn01ha\nwY3/8gvimupY/8SXiHTtIbsfrD9CotaNxuil08wJo5g5YRRfWTyDNw7WsKaglI07yvn920eYEBfD\nsswUVszToiy5gJISVl91O6VxSfxoww99Sb6jfbBTj16CVmNzG3/ZdZQ1W4+waV81be2OS5NHcVtW\nCkszk0kaqUVZ8p6yWZl8aMk3uGH/mzy87j/ee0M9eiV6CQ3Vp06zvrCMNQWlbD9yggiDa6cmsmJe\nCotmjWdYtMbzB7svfm8NL9TAS7/8DKknq3yNIbwfrD80dCNhJWHEe4uy9lfWsaaglLUFZfzbk9sY\nHh3J4tkTWDEvhSsnxw/awlWD2ZbiY6w7Ec0Xk5t8u0fVWcjvBxtI6tFLyGpvd7xZfIw1W33j+XWn\nWxk/KoZlWcmsyEpl+niN5w8Gbe2O5Q//naq60/z1y9cRGz14atpo6EYGlaaWNl7cdZQ1W0t5ZW8V\nre2OWRNGsWJeinbKCnM/eXEvP3lxHw/dmcXSuclehzOgQiLRax699IeaLuP5hR3j+R+cmsiKrBQW\nXTpuUPX4wt2mvVV86n/eZEVWKj/46JxBV1YjJBJ9J/Xopb/srzzF2gLfoqzS443ERkeyePZ4VmSl\nctUlGs8PZWXHG7nloVcZNyqGNf96zaB8IK9EL9JFe7tjS/Ex1hSU8qcd5dQ1tTJu1FCWZaawPDOF\nWcmjvA5ReqClrZ3b/3szeyrqWPeFDw7aejZK9CLdaGpp46VdlawpKOXlPZW0tjtmjB/J8qwUlmUm\nD/qStqHg3ze8y6OvFfGzj2dx65zBNS7flRK9iB+O1TezYbtvPL+g5DhmcNXkeJZnadOUYPXcjnI+\nm7+VT1+dzgNLL/U6HE+FRKLXw1gJJkXV9Wfq7RyqaWDokAhumjWO27J8m6ZERWrTFK8VVdeT81+v\nMSVpBE9/5qpBv5FNSCT6TurRSzBxzlFw+DhrC0pZX1hGbUMLY2KjuGXOBJZnpnD5pDGDbnZHMGhs\nbuO2n/+dipNN/OmL15IyWkNsSvQiAdDS1s6mvVWs3VbGX96toKmlndQxw1iWmczyzBSmqsjagKg/\n3co9T2zhH0XHeOxTH2DhjCSvQwoKSvQiAXbqdCsv7Kxg7bZS/r6/mnYHsyaMYllmMjlzk0lWDzOw\n8vNh1SpOVlTz6U98j8KEDH7dQuVjAAAJPElEQVR0RxbLMlO8jixoKNGL9KPKuiY2FJbzbGEZhYeP\nAzA/fSw5mcncctkExg6P9jjCENexicjx9gg++bFvsyspg4deeIglX71HtWu6UKIXGSDF1fWsLyzj\n2cIy9leeYkiE8cGpCdw6J5lFl45jlGbu9Fx6OjWVtXzi9u9wIH4iv1j7XW44sCWsSw73Rkgkes26\nkXDinGNXeR3rCstYX1hG6fFGoiMjWDAtkZy5E7hh5jhGDFX5BX9Ujozn47d/hyNxSfzyj9/h2uJt\nvjfMoL3d2+CCSEgk+k7q0Uu4cc6x7fBx1heWs3FHORUnmxg6JIKF05O4ec4EPjQjSUm/GyU1DXzy\nm09RFT2CR//wba48vPO9N9WjP0tI1KMXCVdmRlbaGLLSxnD/LTN561AtG7aXsXFHBc+/U0H0kAiu\nm5bIktnjuWHmOOKGaXjHOcfTbx3m2+vfJSIukV8/dT+Xd03ysbG++vLSY+rRiwygtnbH24dq2bij\nnOd3VlBxsomoSOODUxJYdOl4bpiZNCi3SKysa+Jrf9jBX3dXcvUl8fznR+eS8qc/wqpVvj1ftYnI\neWnoRiTItbc7th05znM7ynluZwVHahsxg8yJo7lp1jgWzRrHJYkjwn5x1sYd5axas4OG5ja+tmQG\nn7oqnQhVFfWLEr1ICHHOsbuijr+8e5S/vHuUHaUnAMhIGM7C6UlcPz2R+RljiYkKn1K8tfXNPLD+\nHZ7dVsbc1Dh++LFMpiQNziqUveVZojez4cDPgWbgZedc/sU+o0QvcrbyE428uKuSv7x7lDcO1tDc\n2k5MVARXTY7nummJXD89ifSE4V6H2Sulxxt59NUintxSQnNrO1+8YSr/ev0lDFEtoR4LaKI3s8eA\nW4FK59zsLu2LgZ8CkcCvnHPfN7N/Bo4759ab2VPOudsvdn4lepHuNTa38UZRDa/sqeKVvVUUVdcD\nkDpmGFdOjufKyfFckTGWiWNjz/5gx8rSAR/j7ua675adZPWmA6zfXo4BS+cmc9/1lzBNZSR6LdCz\nbh4Hfgb8ussFIoGHgZuAI8AWM1sHpAI7Og5r60HMInIew6IjWTg9iYXTffVdDtXU88reKl7fX8NL\nu47yzNtHAEgZ7Uv88zPGMKfw70z9X/cxpP6U7ySHDkFuru9/92ey71jRSkMDAPVlR3n9e//Nb4qG\ns+lUFMOjI7nr6nTu/mCGSkYMIL+HbswsHdjQ2aM3s6uAB5xzH+54/fWOQ48Atc65DWb2pHPujoud\nWz16kd5pb3fsrazjjQM1vHHwGP8oqqG2oQWAmJYmLj16kDkV+5hbvo/ZR/eTNiqa6IMH+i+e9Ax2\nNRqvZFzOpox5vJ06k5bIKBIaT3LX8g/wiSsmERerqaSBEvAx+vMk+o8Ai51z/9Lx+p+BK4Cv4uv9\nNwGvdTdGb2a5QC5AWlra5YcOHfIrDhHpXnu7o7imnu0LcygcP5Ud46ewc9wlNEX5pmxGtrcxMXEk\nGQnDmZw4goyE4aTHDyd+RDRjYqMZHRvl1wPflrZ2yo83cbi2gZJjvp9DNfW8uXkX1SPGADCjsojr\nirayoGgr2aXvMrSluV//7oORZwumnHP1wF1+HLcaWA2+Hn2g4xAZjCIijMmJI5h8qojlL70MQKtF\nsC8hjV1JGRRlzOJg5t0UVdXzxsFjNLa8f3R1eHQkozuSPvjm/re2u44/22lpdVSdOk1b+3tf26hI\nI2X0MK6qOcB1r2xiQVEBSfW175100qR+/XvLhfUl0ZcCE7u8Tu1o81uXWjd9CENE3icv78xY+RDX\nzsyqYmbWV8LX74GPzwN8vf+jdU2U1DRQ29DMsfoWahuaqa1v5lhDMycaWjCDyAhjSEREx59GZIQx\nblQMaWNjmTg2lrT4WMaPiiEywiC/DNb9x5kxekArWoNAXxL9FmCqmWXgS/B3AB/vyQmcc+uB9dnZ\n2ff2IQ4ROVfnA9cLzLqJiDAmxA0L7GboflxXBp6/0yt/B1wPJABHgW855x41s5uBn+CbXvmYc65H\nv7ZVvVJEpPe0MlZEJMz5m+g9XYpmZjlmtvrEiRNehiEiEtY8TfTOufXOudy4uDgvwxARCWvq0YuI\nhDn16EVEwpzKxYmIhDkN3YiIhDkN3YiIhLmgmEdvZieA862YigPO7e6f25YAVPdTaBdzvvgG4hz+\nfOZix1zo/e7e8+d+gHf3JBD3ozfn8ff4QN8Tf+/TYLsf/n4mHL4jk5xziRc9yjnn+Q+w2t/2c9uA\nt4It7v4+hz+fudgxF3q/L/fDy3sSiPvRm/P4e3yg74m/92mw3Q9/PzOYviPB8jB2fQ/auzvWC4GI\npTfn8OczFzvmQu8P5vvRm/P4e3yg70lP7pMXvLof/n5m0HxHgmLopi/M7C3nxxJgGTi6J8FF9yP4\nDPQ9CZYefV+s9joAeR/dk+Ci+xF8BvSehHyPXkRELiwcevQiInIBSvQiImFOiV5EJMyFdaI3s+Vm\n9ksze8rMFnkdz2BnZpPN7FEze8brWAYzMxtuZk90fDe0x5/HBuJ7EbSJ3sweM7NKM9t5TvtiM9tj\nZvvN7GsXOodzbq1z7l7gPuD2/ow33AXofhx0zt3Tv5EOTj28PyuAZzq+G0sHPNhBoCf3YyC+F0Gb\n6IHHgcVdG8wsEngYWALMAu40s1lmdpmZbTjnJ6nLR+/v+Jz03uME7n5I4D2On/cHSAUOdxzWNoAx\nDiaP4//96HdDBuIiveGc22Rm6ec0zwf2O+cOApjZk8Ay59z3gFvPPYeZGfB94Dnn3Nb+jTi8BeJ+\nSP/pyf0BjuBL9tsI7s5eyOrh/Xi3v+MJtZucwns9EfD9B5tygeO/ANwIfMTM7uvPwAapHt0PM4s3\ns0eALDP7en8HJ93enz8C/2RmvyAIlucPIue9HwPxvQjaHn0gOOceAh7yOg7xcc7V4HteIh5yztUD\nd3kdh/gMxPci1Hr0pcDELq9TO9rEG7ofwU33J7h4dj9CLdFvAaaaWYaZRQN3AOs8jmkw0/0Ibro/\nwcWz+xG0id7MfgdsBqab2REzu8c51wp8HngB2AU87Zx7x8s4Bwvdj+Cm+xNcgu1+qKiZiEiYC9oe\nvYiIBIYSvYhImFOiFxEJc0r0IiJhToleRCTMKdGLiIQ5JXoRkTCnRC8iEuaU6EVEwtz/B7XUFurP\nnmInAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f09ed16f4a8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "x = np.logspace(-2,1)\n",
    "plt.loglog(x,gamma(x))\n",
    "plt.scatter(np.linspace(1,10,10),factorial(np.linspace(1,10,10)-1),color='r')\n",
    "plt.legend(['$\\\\Gamma(x)$','$(x-1)!$']);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that the Gamma function accepts inputs that are less than one. This scenario is relevant for natural language processing: in many applications, word distributions are sparse, assigning high probability to a few words, and very low probability to everything else. To obtain such sparsity, we can choose a prior $\\vec{\\theta} \\sim \\text{Dirichlet}(a)$, where $a < 1$. \n",
    "\n",
    "A particularly common choice is $a = \\frac{1}{2}$. In this case, the value $\\Gamma(\\frac{1}{2})$ will appear in the denominator of the marginal likelihood for all words that are unseen in the training data.\n",
    "\n",
    "Given that $\\Gamma(a) = (a-1)!$ when $a$ is a positive integer, it may be hard to see what the output should look like when $a \\in (0,1)$. [**Euler's reflection formula**](https://en.wikipedia.org/wiki/Reflection_formula) gives an answer:\n",
    "\n",
    "\\begin{align*}\n",
    "\\Gamma(a) \\Gamma(1-a) = & \\frac{\\pi}{\\sin a \\pi}\\\\\n",
    "\\Gamma\\left(\\frac{1}{2}\\right) = & \\sqrt{\\frac{\\pi}{\\sin \\frac{\\pi}{2}}}\\\\\n",
    "= & \\sqrt{\\pi}\\\\\n",
    "\\end{align*}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 169,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEKCAYAAAAIO8L1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAIABJREFUeJzt3Xd8lfX5//HXlU0GMwNIQsLeECAg\noBZEoDjAgQMVR61ita6ftrbafm2rtdbRqm1d1PoV98YqLlDBxZCwJeyREEYCAcIImef6/ZHAN4QD\nnIScc59xPR+P8+Ak903O++aEXPmM+/MRVcUYY4ypL8zpAMYYY/yTFQhjjDFuWYEwxhjjlhUIY4wx\nblmBMMYY45YVCGOMMW5ZgTDGGOOWFQhjjDFuWYEwxhjjVoTTAU5FYmKiZmZmOh3DGGMCyqJFi3ap\natLJzgvoApGZmUlOTo7TMYwxJqCISJ4n51kXkzHGGLesQBhjjHHLCoQxxhi3rEAYY4xxywqEMcYY\nt6xAGGOMccsKhDHGGLdCskAszt/DI5+tdjqGMcY0mMulPPRxLj9uLfH6a4VkgVi5tYRn52xgfdEB\np6MYY0yDLN9awr+/3cS6ov1ef62QLBCje6UAMCu30OEkxhjTMLNydxAeJozqnuL11wrJAtGuRTP6\npbVgZu4Op6MYY0yDzFxZyGkdW9MiNtLrr+XTAiEi4SKyRERmnOCciSKiIpLtzSxje6WwJH8vRfvK\nvPkyxhjTZDbtOsi6ogOM7eX91gP4vgVxB7DqeAdFJKH2nAXeDjKmV1sAvlhV5O2XMsaYJjGrttdj\ndLAVCBFJA84DXjjBaQ8CjwBe/7W+W0o8GW1irZvJGBMwZq4spHf75qS1ivXJ6/myBfEkcA/gcndQ\nRAYC6ar6sS/CiAhjeqYwd30xB8qrfPGSxhjTaLsOlLMofw9jfNR6AB8VCBE5HyhS1UXHOR4G/B24\n24OvNUVEckQkZ+fOnaeUa2zvtlRUu/h6zal9HWOM8bYvVxWiCmNru8d9wVctiNOBCSKyGXgTGCUi\nr9Y5ngD0AebUnjMU+NDdQLWqTlXVbFXNTko66YZIJzQooxWt46KO9OsZY4y/mpVbSGrLZvRsl+Cz\n1/RJgVDVe1U1TVUzgUnAV6o6uc7xElVNVNXM2nPmAxNU1avbxYWHCWf3SObL1UVUVrvt+TLGGMcd\nLK/im3W7GNs7BRHx2es6eh+EiDwgIhOczDC2d1v2l1WxYONuJ2MYY8xxfbtuJxVVLp92L4EDe1Kr\n6hxgTu3z+49zzkhf5TmjSyIxkWHMyt3BGV0TffWyxhjjsZm5hbSMjWRwZiufvm5I3kldV7OocH7S\nNYlZuYWoqtNxjDHmKFXVLr5aXcSoHslEhPv2R3bIFwiAMb1S2FZSxspt+5yOYowxR1m4eQ97Syt9\ndvd0XVYggLN7phAmMHOlzWYyxviXmbk7iI4I4yfdTm3WZmNYgQBax0WRndmamba6qzHGj6gqM1cW\nckaXRGKjfD5kbAXisLG9Uli9Yz/5xaVORzHGGABWbd/P1r2HGNvb991LYAXiiMPTx2xtJmOMv5iZ\nuwMRGNXDCoSjOrSJpUfbBNtEyBjjN2blFjKoQyuSEqIdeX0rEHWM7ZXCws272X2wwukoxpgQV7Cn\nlJXb9jnWvQRWII4ypldbXApfrbY9IowxzvqitjdjjI/vnq7LCkQdfVKb065FjE13NcY4bmZuIV2T\n4+mYGOdYBisQdYgIY3ql8M26nRyqqHY6jjEmRO0trWDBpt0+3fvBHSsQ9Yzt1ZayShffrd/ldBRj\nTIiavaaIapcytrdz3UtgBeIYp3VqTUJMhHUzGWMcM3NlIckJ0fRLbeFoDisQ9USGhzGqdo+Iapct\n3meM8a2yymq+XruTMb1SCAvz3d4P7liBcGNMrxR2H6xgUd4ep6MYY0LM3A27KK2odnz8AaxAuDWi\nWxJR4WG2Fakxxudm5RYSHx3BsM5tnI5iBcKdhJhIhndpw0zbI8IY40PVLmVWbiEjuycRHRHudBwr\nEMczplcKecWlrC084HQUY0yIWLplD7sOVPhF9xJYgTiuMT1r3iDrZjLG+MrM3EIiw4WzeiQ7HQWw\nAnFcyc1jyEpvaXtEGGN8ZtbKQoZ2akPzmEinowBWIE5obO8UlheUsL3kkNNRjDFBbn3RATbuOujI\n1qLH49MCISLhIrJERGa4OXaXiOSKyHIR+VJEMnyZzZ3Db5QtAW6M8bbDe9GMDtUCAdwBrDrOsSVA\ntqr2A94FHvVZquPonBRP1+R4Ply6zekoxpgg9+HSbWSlt6Rdi2ZORznCZwVCRNKA84AX3B1X1dmq\neni/z/lAmq+yHY+IcOGAVHLy9thWpMYYr8ndto/VO/Zz8cBUp6McxZctiCeBewCXB+f+HPjUu3E8\nc+GAmjds+pKtDicxxgSr6UsKiAgTzu/X3ukoR/FJgRCR84EiVV3kwbmTgWzgseMcnyIiOSKSs3Pn\nziZOeqzUls0Y2qk105cU2E1zxpgmV+1S/rt0GyO7J9M6LsrpOEfxVQvidGCCiGwG3gRGicir9U8S\nkdHA74AJqlru7gup6lRVzVbV7KSkJG9mPuLigWlsLi5lyZa9Pnk9Y0zo+H79Lor2lzPRz7qXwEcF\nQlXvVdU0Vc0EJgFfqerkuueIyADgeWqKg1/t+XlOn7ZER4QxfbF1Mxljmtb0JVtpHhPBqJ7+cXNc\nXY7eByEiD4jIhNoPHwPigXdEZKmIfOhgtKMkxEQytndbPlq+jYoqT4ZQjDHm5A6WV/HZjzs4r197\nv1h7qb4IX7+gqs4B5tQ+v7/O50f7OktDXDwglY+WbWP2miJ+6vAuT8aY4PDZjzs4VFntd7OXDrM7\nqT10ZtdEEuOjrJvJGNNkpi/ZSnrrZmRntHI6iltWIDwUER7G+P7t+Wp1ESWllU7HMcYEuB0lZXy/\nYRcXZaUi4uzOccdjBaIBJg5Mo6LaxYwVdme1MebU/HfpVlThooGO3xN8XFYgGqB3++Z0TY63biZj\nzCmbvmQrAzq0pGNinNNRjssKRAOICBcNtKU3jDGn5sjSGgP8c3D6MCsQDXRhVioitvSGMabxpi8p\nIDLc/5bWqM8KRAO1b9mMoR3b8L4tvWGMaYSqahcf1C6t0crPltaozwpEI1w0MJW84lIW59vSG8aY\nhvl+QzE795f7ffcSWIFolHP6tCUmMozpSwqcjmKMCTDTFxf47dIa9VmBaISEmEjG9mrLjOXbbekN\nY4zHDpZX8fnKQs7v759La9RnBaKRLhqYyt7SSmav8at1BY0xfuzI0hoB0L0EViAa7cwuiSTGR9s9\nEcYYj01fspUOrWMZ5KdLa9RnBaKRIsLDmFC79Mbe0gqn4xhj/NzhpTUuHOC/S2vUZwXiFFw8MLVm\n6Y3l252OYozxcx8cXlojQLqXwArEKendvjndUuLtpjljzAmpKtMXb2Wgny+tUZ8ViFMgIlw0II1F\neXvIKz7odBxjjJ/K3b6PNYX7/XphPnesQJyiCwe0t6U3jDEnNH3x1pqlNfq2czpKg1iBOEXtWjRj\nWKc2TF+y1ZbeMMYco6raxX+XbeOsAFhaoz4rEE3gogG29IYxxr0jS2v46baiJ2IFogmc07cdMZFh\nvL/Ylt4wxhzt/cUFtGgWyVk9/H9pjfqsQDSB+OiII0tvlFVWOx3HGOMn9pdV8vnKHZzXr11ALK1R\nnxWIJjJpcDolhyrtnghjzBHvLSqgrNLF5dnpTkdpFJ8WCBEJF5ElIjLDzbFoEXlLRNaLyAIRyfRl\ntlM1rHMbuibHM23uZhusNsbgcikvz8sjK70l/dNbOh2nUXzdgrgDWHWcYz8H9qhqF+AJ4BGfpWoC\nIsI1wzNZsbXEBquNMXy7fhcbdx3kuuGZTkdpNJ8VCBFJA84DXjjOKRcA02qfvwucLYGyYEmtiwek\nkhAdwcvzNjsdxRjjsJfnbiYxPppzA+zeh7p82YJ4ErgHON4GCqnAFgBVrQJKgDa+idY04qIjuCQ7\njU9WbKdof5nTcYwxDskvLuWrNUVcOSSdqIjAHer1SXIROR8oUtVFTfC1pohIjojk7Ny5swnSNa1r\nhmVSWa28viDf6SjGGIe8PG8z4SJcNTTD6SinxFel7XRggohsBt4ERonIq/XO2QqkA4hIBNACKK7/\nhVR1qqpmq2p2UlKSd1M3QsfEOEZ2T+K1Bfm225wxIai0ooq3c7Ywrk9bUprHOB3nlPikQKjqvaqa\npqqZwCTgK1WdXO+0D4Fra59fUntOQE4HunZ4Jjv3l/PZyh1ORzHG+NgHS7axr6wqoAenD3O0c0xE\nHhCRCbUf/gdoIyLrgbuA3zqX7NSM6JpEZptYps3d7HQUY4wPqSrT5m6mV7vmAbNr3Ik0uECISJyI\nNPqWQFWdo6rn1z6/X1U/rH1epqqXqmoXVR2iqhsb+xpOCwsTrh6WyaK8Pfy4tcTpOMYYH5m/cTdr\nCvdz3fDMgNk17kROWiBEJExErhSRj0WkCFgNbBeRXBF5TES6eD9m4Lk0O43YqHBeslaEMSFj2tzN\ntIqNZEJWe6ejNAlPWhCzgc7AvUBbVU1X1WTgDGA+8IiI1B9PCHnNYyK5eGAqHy7bxu6Dtme1McFu\n695DzMzdweWDOxATGXjrLrnjSYEYraoPqupyVT0yLUdVd6vqe6o6EXjLexED1zXDMqmocvHmQpvy\nakywe21+HgCTh3ZwOEnTOWmBUNVKABFpIyI3i8jPRGSIiDSrf445WreUBIZ3bsOr8/KoqrYpr8YE\nq7LKat5cuIXRPVNIaxXrdJwm05BB6ulAEvAX4DGgRERWeyVVELl2eCbbSsr4YlWh01GMMV7yUW1X\ncjBMba2rIQUiQVUfAApVdQRwBfC2d2IFj9E9U0ht2Yxpc/OcjmKM8QJVZdq8zXRNjmdY54BaHeik\nGlIgDi8uVC4izVT1PWCsFzIFlfAwYfLQDOZtLGbNjv1OxzHGNLHF+Xv5ces+rgmSqa11NaRAPC4i\nrakZkH5RRG4DAnORcx+bNDid6Igwps3b7HQUY0wTmzZ3MwkxEVw8IPD2nD4ZjwtE7Yyl3ar6d+AT\natZNuthryYJIq7goLshqz/TFWykptfF8Y4JF0b4yPlmxnUsHpRMXHeF0nCbXqKU2VPUVVb1HVXOb\nOlCwumZYJocqq3ln0RanoxhjmsjrP+RT5VKuGRbYq7YeT+AuVB5g+qS2IDujFa/Mz8PlCsg1CI0x\ndVRUuXhtQT4juyeRmRjndByvaMxaTOO9ESQUXDs8k7ziUuasLXI6ijHmFH3643Z27i/n2iCb2lpX\nY1oQDzV5ihBRsz58NC/ZlFdjAt60uZvpmBjHiK7+ty9NU2lMgQiueVw+FBkexlWnZfDN2p1s3HnA\n6TjGmEZaUVDC4vy9XD00g7Cw4P2R2JgCYR3op+CKIR2IDBf+890mp6MYYxrpP99tJDYqnEuy05yO\n4lU2SO1jSQnRXDIonXdyCthecsjpOMaYBtq48wAfLtvG5KEZNI+JdDqOV1mBcMAtIzvjUuW5ORuc\njmKMaaB/zV5PVEQYN57ZyekoXteYAmGrzp2i9NaxXDIojTcWbqFwX9nJ/4Ixxi9s3nWQ/y7dxuTT\nMkhKiHY6jtc1uECo6hhvBAk1vzyrC9Uu5VlrRRgTMP41ez0RYcKUEcHfegDrYnJMeutYLh6Qyhs/\n5FNkrQhj/F5+cSnTl2zlytM6kJwQ43Qcn7AC4aBbR3WhyqU8/81Gp6MYY07i6dnrCQ8TfjGis9NR\nfMYnBUJEYkTkBxFZJiIrReRPbs7pICKzRWSJiCwXkXN9kc1JGW3iuDArldcW5LFzf7nTcYwxx7Fl\ndynvLS7gyiEdSGkeGq0HOMUCISI/8/DUcmCUqvYHsoBxIjK03jm/B95W1QHAJOCZU8kWKG4d1YWK\nKhdTv7GxCGP81TNz1hMmodV6gFNvQRzTEnBHaxy+dTiy9lH/hjsFmtc+bwFsO8VsAaFjYhwXZKXy\nyvw8dh2wVoQx/qZgTynvLirg8sHptG0ROq0H8KBA1Hb3uHusAFI8fSERCReRpUARMEtVF9Q75Y/A\nZBEpoGa/ids8v4zAdrgV8W8bizDG7xyeaXjzyNBqPYBnLYgU4BpgvJtHsacvpKrVqpoFpAFDRKRP\nvVOuAF5S1TTgXOAVETkmn4hMEZEcEcnZuXOnpy/v1zonxTO+f3tenpdHsbUijPEb2/Ye4u2cLVyW\nnU77ls2cjuNznhSIGUC8qubVe2wG5jT0BVV1LzAbGFfv0M+Bt2vPmQfEAIlu/v5UVc1W1eykpOBZ\nRfG2UV0oq6rmBVujyRi/cbj1cMtZXRxO4gxPCsQNqvqduwOqeiWAnGSnbhFJEpGWtc+bAWOA1fVO\nywfOrj2nJzUFIjiaCB7okpzAeX3b8fLczew5WOF0HGNC3o6SMt5auIVLBqWRGoKtB/CsQMwWkdtE\npEPdT4pIlIiMEpFpwLUn+Rrtar/OcmAhNWMQM0TkARGZUHvO3cCNIrIMeAO4TlVDauXY28/uSmll\nNS98Z2MRxjjtua834FLllpGh2XoA8GSX7XHA9cAbItIR2EvNb/fhwEzgSVVdcqIvoKrLgQFuPn9/\nnee5wOmeRw8+3VISOLdPO6bNzePGMzvRMjbK6UjGhKTCfWW8/kM+Ewemkd461uk4jjlpC0JVy1T1\nGVU9HcigphtooKpmqOqNJysOpmFuO7sLB8qreNHGIoxxzHNfb6DapfwyRMceDvNkmuuUw89VtVJV\nt9cONBsv6NG2OeN6t+V/v99MSWml03GMCTlF+8t4fUE+Fw1IpUOb0G09gGdjEL/wegpzlNvP7sr+\n8ipe/N5aEcb42tSvN1LlUm4N8dYDeFYggnfDVT/Vq31zxvZK4cXvN1FyyFoRxvjKzv3lvLogjwuy\n2pOZGOd0HMd5UiBCaiaRv7j97K7sL6ti2tzNTkcxJmT8+9uNVFS5rPVQy5MC0V9E9rl57BeRfV5P\nGKL6pLZgdM9kXvh2I3tL7b4IY7ytaF8Zr8zLY0L/9nRKinc6jl/wZBZTuKo2d/NIUNXmJ/v7pvHu\nHtudA+VVPPnFOqejGBP0Hv18DdUu5c7R3ZyO4jcavJqriMSJSLg3wpij9WzXnCtP68Ar8/NYW7jf\n6TjGBK1lW/by7qICrj+jo4091OHJNNcwEblSRD4WkSJgDbBDRHJF5DERsc46L7prTHfiosJ5cEYu\nIXZjuTE+oar88aOVJMZHc+so+3FWl0dLbQCdgXuBtqqapqpJwBnAfOAREZnsxYwhrXVcFHeO7sa3\n63bx5aoip+MYE3T+u3QbS/L3cs+47sRHe7K4ROjw5F9jtKoeM9dSVXcD7wHviUhkkyczR1w9LIPX\nFuTx549zObNbItER1sNnTFMorajir5+upm9qCy4ZmOZ0HL/jySB1JYCItBGRm0XkZyIypHZV1qPO\nMd4RGR7G/eN7s7m4lJe+3+x0HGOCxnNzNrBjXxl/nNCLsDC75au+hgxSTweSgL8AjwElIlJ/yW7j\nJSO6JXF2j2T++dV6ivaXOR3HmIC3ZXcpz3+zkQuy2jMoo7XTcfxSQwpEgqo+ABSq6ghqdoB72zux\njDu/O68n5VXVPP75GqejGBPw/vrpasJE+O05PZyO4rcaUiAO/9paLiLNVPU9YKwXMpnj6JQUz89O\n78g7iwpYXmDrJRrTWPM3FvPxiu3cPLIz7VqE5mZAnmhIgXhcRFoDbwEvishtQEvvxDLHc+uoLrSJ\ni+KBj2zaqzGNUe1S/vRRLqktmzHlJ52cjuPXPC4Qqvqequ5W1b8DnwDpwMVeS2bcah4Tya/Gdicn\nbw8fLd/udBxjAs5bC7ewavs+7j23BzGRNiPwRDy5Ue6YoX1VfUVV76ndBe6ke1KbpnVpdjq92zfn\n4U9Wcaii2uk4xgSMkkOVPD5zDUMyW3Ne33ZOx/F7vtqT2jSh8DDhD+N7s72kjOe+3uB0HGMCxj++\nXMee0gruH98L+7325DwpEOOAamr2pN5Wu8TGJmAdNTOZnlTVl7yY0bgxpGNrzu/Xjue+3sDWvYec\njmOM31tfdIBpczczaXA6fVJbOB0nIHhSIK5xsyf1ANuT2nn3ntsTgIc/WeVwEmP8358/zqVZZDh3\nj+3udJSA0aAtR21Pav+S2rIZvxjRmRnLt/PDpt1OxzHGb81eXcScNTu5Y3RXEuOjnY4TMHyy5aiI\nxIjIDyKyTERWisifjnPeZbVdWCtF5PVTfd1Q8IsRnWnXIoY/fbSSapdNezWmvooqFw/OyKVTYhzX\nDMt0Ok5A8dWWo+XAKFXtD2QB40RkaN0TRKQrNSvGnq6qvYE7m+B1g16zqHDuPbcnK7ft482F+U7H\nMcbvvDR3Ext3HeR/zu9FVESDt8AJaT7ZclRrHKj9MLL2Ub/w3Ag8rap7av+OrW3tofH92jGsUxse\n/mQ1BXtKnY5jjN/YsPMAf5u5ltE9kzmrR7LTcQKOz7YcFZFwEVkKFAGzVHVBvVO6Ad1E5HsRmS8i\n447zdaaISI6I5OzcudPTlw9qIsKjl/RDVbnn3eW4rKvJGKqqXdz99jKaRYXzl4v6Oh0nIPmsvaWq\n1aqaBaQBQ0SkT71TIoCuwEhqps/+W0SOWcpDVaeqaraqZiclJXk7dsBIbx3L787rxdwNxby6IM/p\nOMY47vlvNrJ0y14euKAPyc1jnI4TkHzeIVc7A2o2NfdX1FUAfFg7U2oTsJaagmE8dMWQdH7SLYmH\nP1nN5l0HnY5jjGNW79jHk1+s5dy+bRnfz+6YbiyfFAgRSTrcGqjdaGgMUH8viQ+oaT0gIonUdDlt\n9EW+YCEiPDKxLxHhwq/fXWazmkxIqqztWmoeE8mDF/SxO6ZPga9aEO2oWbJjObCQmjGIGSLygIhM\nqD3nc6BYRHKpaWH8WlWLfZQvaLRr0Yw/jO/Nws17+N/vNzkdxxif+9dX61m5bR8PXdSXNnbPwynx\nyQ7dqrocGODm8/fXea7AXbUPcwomDkzlsx+38+jnaxjZPZkuyfFORzLGJ37cWsLTs9dzYVZ7xvVp\n63ScgGeTgoOQiPCXi/sSGxXO3e8so6ra5XQkY7yuvKqau95eSuu4KP40of4cGNMYViCCVHJCDA9e\n0IdlW/by/Dc2lGOC35NfrGNt4QEemdiPFrGRTscJClYggtj4/u05r287nvxiLat3eHRPozEBaXH+\nHp7/egOXZ6fbDXFNyApEkHvwwj60aBbJ3W8vo9K6mkwQKqus5lfvLKNdi2b8/vyeTscJKlYgglzr\nuCgeuqgvK7ft419frXc6jjFN7rHP17Bx50EevaQfCTHWtdSUrECEgJ/2bstFA1L51+z1rCgocTqO\nMU1mwcZiXvx+E1cPzeD0LolOxwk6ViBCxB/H9yYxPoq731lKeZXtY20C38HyKn717jLSW8Xy23N6\nOB0nKFmBCBEtYiP568R+rC08wBOz1jkdx5hT9vCnqyjYc4jHL+1PXLRPbukKOVYgQshZ3ZOZNDid\n57/ZwNdrbSVcE7g+XbGdV+fnc/3pHRnSsbXTcYKWFYgQc//4XnRPSeD2N5aQX2x7R5jAs7ZwP3e/\ns4wBHVpyzzjbX9qbrECEmNioCKZenQ3AlFdyKK2ocjiRMZ4rOVTJTa8sIjYqgmevGkR0RLjTkYKa\nFYgQ1KFNLP+4YgBrCvfzm/dWULMMljH+zeVS/t9bS9myu5RnJw+kbQvb48HbrECEqBHdkvj1T7vz\n0bJtvPCtrfpq/N+TX67jq9VF3D++F4MzbdzBF6xAhLCbR3TmnD5tefjTVcxdv8vpOMYc16zcQv7x\n5TouGZTG1UMznI4TMqxAhDAR4bFL+9M5KZ5fvr6Ygj02aG38z/qiA/y/t5bSL60Ff77QNgDyJSsQ\nIS4+OoKp12RT5VJuemURZZV2E53xH/vLKrnplRyiIsJ4dvIgYiJtUNqXrEAYOibG8dSkLHK37+O+\n6TZobfyDy6Xc/fYyNheX8q8rB5DaspnTkUKOFQgDwKgeKdx5djfeX7yVaXM3Ox3HGJ6Zs56ZuYXc\nd25Phne2dZacYAXCHHHbqC6M7pnCgx+vYv5G2w7cOGf26iL+NmstF2S15/rTM52OE7KsQJgjwsKE\nv1/en4w2sdz6+mK2lxxyOpIJQZt3HeT2N5fQo21z/npxPxuUdpAVCHOU5jGRTL06m7JKF794dbEN\nWhufOlhexU2vLCI8TJh69SCaRdmgtJN8UiBEJEZEfhCRZSKyUkT+dIJzJ4qIiki2L7KZY3VJjudv\nl/Vn2Za93P3OMqpdNmhtvK+y2sVtbyxhXdF+/nnFANJbxzodKeT5qgVRDoxS1f5AFjBORIbWP0lE\nEoA7gAU+ymWO46e923LvOT34ePl27nvfZjYZ76quXUbjq9VFPHBBH87smuR0JAP4ZBF1rfnpcqD2\nw8jah7ufOA8CjwC/9kUuc2I3jejMgfIq/vnVeuKiI/if83taf7BpcqrKfe+vYMby7dx7Tg8m253S\nfsNnYxAiEi4iS4EiYJaqLqh3fCCQrqof+yqTObm7xnTjZ6dn8uL3m3jyC9toyDQtVeXBGat4K2cL\nt43qwk0jOjsdydThs22YVLUayBKRlsB0Eemjqj8CiEgY8HfgupN9HRGZAkwB6NChg/cCG6BmOY7/\nOa8XB8ureOrLdcRHR3DjTzo5HcsEiSe+WMeL32/iuuGZ3DWmm9NxTD0+n8WkqnuB2cC4Op9OAPoA\nc0RkMzAU+NDdQLWqTlXVbFXNTkqyfkpfCAsTHr64H+f1bcdDn6zi9QX5TkcyQWDqNxv4x5fruHRQ\nGvef38u6L/2QT1oQIpIEVKrqXhFpBoyhZqwBAFUtARLrnD8H+JWq5vginzm58DDhicuzKK2o4ncf\nrCAuOpwLslKdjmUC1OsL8vnLJ6s5r287/jqxH2FhVhz8ka9aEO2A2SKyHFhIzRjEDBF5QEQm+CiD\nOUWHF0wbktmau95exsyVO5yOZALQf5du5XcfrOCs7kk8cXkW4VYc/JYE8vTF7OxszcmxRoavHSiv\n4qoXFrBq2z5evG4wZ3S1dXKMZ2au3MHNry0mO6MV064fYquzOkREFqnqSe81szupTYPFR0cw7WeD\n6ZQUx40v57Aob7fTkUwA+G5ymfLcAAAP50lEQVTdLm59fQl9Ulvwn+sGW3EIAFYgTKO0jI3i5Z8P\noW2LGK7734X8uLXE6UjGjy3K282NL+fQKSmOaT8bTHy0zyZQmlNgBcI0WnJCDK/ecBoJ0RFc/Z8F\nLM7f43Qk44fmbtjFdS8upG2LGF7++RBaxkY5Hcl4yAqEOSWpLZvx+o1Dad4skiumzuezH23g2vyf\nD5Zs5doXf6Btixheu+E0khNinI5kGsAKhDllmYlxvHfzcHq2a87Nry3ipe83OR3JOExVeXr2eu58\naymDMlrx7i+G0952hAs4ViBMk0iMj+aNG4cyumcKf/wolz/PyMVlq8CGpKpqF/dN/5HHPl/DBVnt\nmXb9EFrERjodyzSCFQjTZJpFhfPc5EFcOyyDF77bxK1v2H4SoeZgeRU3vpzDGz/kc8vIzjxxWRbR\nETZbKVDZVALTpMLDhD9O6E1aq1ge+mQVRfsW8O9rsmkVZwOTwa5ofxnXv7SQ3G37eOiiPlx1mq3K\nGuisBWGanIhw40868fSVA1m+tYSJz84lv7jU6VjGi9YX7eeip+eyoeggL1ybbcUhSFiBMF5zXr92\nvHbDaewureCiZ75n6Za9TkcyXrBgYzEXPzOX8ioXb900lFE9UpyOZJqIFQjjVYMzW/PezcOJjQ5n\n0tR5zMotdDqSaUIfLtvG1f/5gaSEaKbfMpx+aS2djmSakBUI43Wdk+J5/+bT6Z6SwE2v5PDMnPW2\nz3WAq6p28feZa7j9jSVkpbfkvZuH2x7SQcgKhPGJpIRo3pgylHP6tOPRz9Zw1Qvz2bb3kNOxTCPk\nF5dy6fPz+MdX65k4MM3ujg5iViCMz8RGRfCvKwfw6CX9WF5Qwrgnv+Hj5dudjmU8pKq8t6iAc576\nhvVFB/jHFQP422X9bdG9IGbTXI1PiQiXZaczJLM1d7y1lF++vpjZa9L444TetoCbHyspreS+D1bw\n8fLtDOnYmicuzyLV7owOetaCMI7ITIzj3V8M4/ZRXXh/cQHnPvWtLfbnp+ZtKGbcU9/w+Y87uGdc\nd964cWjQFoekpCRExOPHFVdc4XRkr7ICYRwTGR7GXWO789ZNw6h2KZc+N4+nvlhHVbXL6WgGqKhy\n8chnq7nyhfnERIbz/i3DuWVkl6DdAS4vL4/w8HBiY2sG22NiYkhJSTnqERcXB0B0dDQpKSmcddZZ\nTkb2OisQxnGDM1vz6Z1nMr5fO574Yi2XT53Plt12Y52TNuw8wMRn5/LsnA1MGpzOx7efEfRTWDMy\nMtixYwcTJtTsgvzLX/6SHTt2HPW45JJLALjhhhvYsWMHU6ZMcTKy11mBMH6heUwkT04awFOTsli7\nYz/nPPUt7y4qIJC3xA1ELpfy2oI8zv/HdxTsKeW5yYN4+OJ+xEaFzvjQihUrAOjbt+8xx5YtW3bc\nY8EodN51ExAuyEplYIdW3PX2Un71zjJenZ/H78/rSXZma6ejBb2563fx0CerWLltH2d0SeRvl/Un\npXlo7d9QUVHBmjVrgGOLQFVVFatWrXJ7LFhZgTB+J711LG9OGcZ7iwv428w1XPLcPM7p05bfjOtB\nZmKc0/GCzvqi/Tz8yWq+XF1EastmPDUpi/H92hMWpGMNJ7Jq1SqqqqoIDw+nV69eRx1bvXo15eXl\nAPTp08eJeD7nkwIhIjHAN0B07Wu+q6p/qHfOXcANQBWwE7heVfN8kc/4n/Cwmumw5/drxwvfbuK5\nrzfwxapCJg/N4PZRXW112Cawc385T36xljcXbiE2MpzfntOD64ZnhvR9DYe7l7p27UpMzNGtp8Pd\nSxkZGTRv3tzn2ZzgqxZEOTBKVQ+ISCTwnYh8qqrz65yzBMhW1VIRuRl4FLjcR/mMn4qNiuD2s7sy\naUg6T8xay7S5m3l3UQG3jerCtcMzba+BRjhUUc1/vtvIs3M2UF7l4uqhGdx+dldaW9Fl+fLlgI0/\nHOaTAqE1I40Haj+MrH1ovXNm1/lwPjDZF9lMYEhOiOHhi/tx3fCOPPzpKv7yyWpenpfHb8b14Px+\n7RAJve6QhnK5lOlLtvL4zDVsLyljbK8UfntODzolxTsdzW8cbkH069fvmGMnKh7BymdjECISDiwC\nugBPq+qCE5z+c+BTnwQzAaV72wRe+tkQvl23k4c+XsVtbyzhP99t4vazuzCiW3LQztE/FVXVLr5Y\nVcQ/vlxH7vZ99EtrwZOXZ3FapzZOR/M7NoPpaOLraYQi0hKYDtymqj+6OT4ZuBUYoarlbo5PAaYA\ndOjQYVBeng1ThKpql/Le4gL+PnMtO/aVkdqyGVee1oHLstNJSoh2Op7jdpSU8ebCfN78YcuRf597\nxnUP2QHok9m9ezdt2tQUzQ0bNtCpU6cjx4qLi0lMTARqWhKBXiREZJGqZp/0PCfmmYvI/UCpqj5e\n7/OjgX9SUxyKTvZ1srOzNScnx0spTaCoqHIxM3cHr83PZ97GYiLDhZ/2bsvkoRmc1rF1SHU/uVzK\n9xt28dr8fGatKqTapZzZNZHJQzM4u0cyEeF269PxfP3114wcOZL4+Hj27dt31PdNTk4OgwcPBqCw\nsJDk5GSnYjYJTwuEr2YxJQGVqrpXRJoBY4BH6p0zAHgeGOdJcTDmsKiIMM7v157z+7VnfdEBXluQ\nx3uLCpixfDtdk+O56rQOXDwojeYxkU5H9Zo9Byt4d1EBr/+Qz6ZdB2kVG8kNZ3TkytM6kNHGpgZ7\n4nD3Up8+fY75pSIs7P8K6x133MG4ceO49tprfZrPCb4ag2gHTKsdhwgD3lbVGSLyAJCjqh8CjwHx\nwDu1b06+qk7wUT4TJLokx/OH8b2556c9+Gj5Nl6bn8cfP8rlkc/WcEFWey4fnE7/tJZB0cVS7VIW\n5+/hjR/ymbF8OxVVLgZltOL2s7twTp92IT1dtTFONP7Qr18/xo8fz2effcaMGTPo2bOnr+M5wpEu\npqZiXUzGEysKSnh1fh7/XbaVskoXrWIjGd4lkTO7JHJG10TSWgXOTmh5xQf5dt0uvlu3i7kbdrGv\nrIq4qHAuGpjKVadl0LNdaMzPN6fGr8cgmooVCNMQJYcqmb26qOYH7PqdFO6rmQPRMTGOM2qLxbDO\nbfyqK2pvaQVzNxQfybxld80ufKktmx3JfFaPZNtLwzSIFQhjTkBVWV90oPYH7y7mbyymtKKa8DCh\nf1oLzuiaRK92zenQOpYObWJ98gN4X1kl+cWlbNldyo/bSvhu3S5WbC3BpRAfHcHQTm34SbdEzuiS\nSMfEuJAafDdNywqEMQ1QUeViSf4evlu/i2/W7WJFwV5cdf5rtI6LqikWdR9tav5MTogmzIMf1i5V\nduwrI393TRHIKy498jx/dyl7SiuPnBseJmSlt+SMLomc2TWR/uktibQZSKaJWIE4mZEjj/3cZZfB\nLbdAaSmce+6xx6+7ruaxaxfUrgt/lJtvhssvhy1b4Oqrjz1+990wfjysWQM33XTs8d//HkaPhqVL\n4c47jz3+l7/A8OEwdy7cd9+xx598ErKy4Isv4M9/Pvb4889D9+7w0Ufwt78de/yVVyA9Hd56C559\n9tjj774LiYnw0ks1j/o++QRiY+GZZ+Dtt489PmdOzZ+PPw4zZhx9rFkz+LT23sgHH4Qvvzz6eJs2\n8N57Nc/vvRfmzTv6eFoavPpqzfM776z5N6yrWzeYOrXm+ZQpsHbt0cezsmr+/QAmT6YqfwvlldWU\nVbkoq6xmXae+vHDujeTvLuW+F39Py9J9R/317zP688/Ta3YXe+ntPxBTdfQtPF92HsK/T7sYgDdf\n/y1Qs/1qVEQYMRFhrBtxDtuvup7MZnDGHdcSExl+9E1/9r0XMt97FBQcfXzYMHj44ZrnEydCcfHR\n19QIfjXN1ZhAExEmRERHEFd7v11q9yRG3nAaAK75KVQUhlNeVU1ZpYvKahfDOycSMbobAJlfxRJR\nfvQMop90SyJhTDeSE6LpObs50ZHhRIeHcbjh0TMrFUZ0rvnlxMYTjJ8I3RaEMcaEKE9bENapaYwx\nxi0rEMYYY9yyAmGMMcYtKxDGGGPcsgJhjDHGLSsQxhhj3LICYYwxxi0rEMYYY9wK6BvlRGQnEIh7\njiYCu5wO4WOhds2hdr1g1xxIMlQ16WQnBXSBCFQikuPJXYzBJNSuOdSuF+yag5F1MRljjHHLCoQx\nxhi3rEA4Y6rTARwQatccatcLds1Bx8YgjDHGuGUtCGOMMW5ZgfASERknImtEZL2I/NbN8V+IyAoR\nWSoi34lILydyNqWTXXOd8yaKiIpIwM/+8OB9vk5Edta+z0tF5AYncjYlT95nEblMRHJFZKWIvO7r\njE3Ng/f5iTrv8VoR2etEzianqvZo4gcQDmwAOgFRwDKgV71zmtd5PgH4zOnc3r7m2vMSgG+A+UC2\n07l98D5fB/zL6aw+vuauwBKgVe3HyU7n9vY11zv/NuBFp3M3xcNaEN4xBFivqhtVtQJ4E7ig7gmq\nWndT4zgg0AeDTnrNtR4EHgHKfBnOSzy95mDiyTXfCDytqnsAVLXIxxmbWkPf5yuAN3ySzMusQHhH\nKrClzscFtZ87ioj8UkQ2AI8Ct/som7ec9JpFZCCQrqof+zKYF3n0PgMTRWS5iLwrIum+ieY1nlxz\nN6CbiHwvIvNFZJzP0nmHp+8zIpIBdAS+8kEur7MC4SBVfVpVOwO/AX7vdB5vEpEw4O/A3U5n8bGP\ngExV7QfMAqY5nMcXIqjpZhpJzW/T/xaRlo4m8p1JwLuqWu10kKZgBcI7tgJ1f1NMq/3c8bwJXOjV\nRN53smtOAPoAc0RkMzAU+DDAB6pP+j6rarGqltd++AIwyEfZvMWT7+0C4ENVrVTVTcBaagpGoGrI\n/+dJBEn3EliB8JaFQFcR6SgiUdR803xY9wQRqfsf5jxgnQ/zecMJr1lVS1Q1UVUzVTWTmkHqCaqa\n40zcJuHJ+9yuzocTgFU+zOcNJ71m4ANqWg+ISCI1XU4bfRmyiXlyzYhID6AVMM/H+bwmwukAwUhV\nq0TkVuBzamZAvKiqK0XkASBHVT8EbhWR0UAlsAe41rnEp87Daw4qHl7z7SIyAagCdlMzqylgeXjN\nnwNjRSQXqAZ+rarFzqU+NQ343p4EvKm1U5mCgd1JbYwxxi3rYjLGGOOWFQhjjDFuWYEwxhjjlhUI\nY4wxblmBMMYY45YVCGOMMW5ZgTDGGOOWFQhjmpiIXFK7SN2y2r0+kpzOZExj2I1yxjQxEWlz+M5h\nEfkDsEtVn3Y4ljENZi0IY5redSLyg4gsA24hOPa+MCHI1mIypgmJyDXUbDAzSlUPiMg3wEqHYxnT\nKNaCMKZp9QXm1haHicBwYIXDmYxpFBuDMKYJiUhv4H2gBJgJXKaq3ZxNZUzjWIEwxhjjlnUxGWOM\nccsKhDHGGLesQBhjjHHLCoQxxhi3rEAYY4xxywqEMcYYt6xAGGOMccsKhDHGGLf+P11k+HmY8gjM\nAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f09ed350a90>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "min_a = .25\n",
    "max_a = .75\n",
    "a_vals = np.linspace(min_a,max_a,30)\n",
    "gamma_a = [gamma(a)*gamma(1-a) for a in a_vals]\n",
    "plt.plot(a_vals,gamma_a)\n",
    "plt.plot([min_a,max_a],[np.pi,np.pi],'r--')\n",
    "plt.text(0.9*max_a,1.01*np.pi,'$\\\\pi$',fontsize=24)\n",
    "plt.xlabel('$a$')\n",
    "plt.ylabel('$\\\\Gamma(a)\\\\Gamma(1-a)$');"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}