olgabot · March 4, 2014 19:19
diff --git a/botvinnik_bnfo285_hw3.ipynb b/botvinnik_bnfo285_hw3.ipynb
 {
 "metadata": {
  "name": ""
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "# Problem 1\n",
      "\n",
      "## Problem 1 (a)"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import prettyplotlib as ppl\n",
      "\n",
      "p = 0.3\n",
      "n = 100000\n",
      "data = np.random.geometric(p=p, size=n)\n",
      "ppl.hist(data, bins=20)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 4,
       "text": [
        "<matplotlib.axes._subplots.AxesSubplot at 0x1092861d0>"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAFBtJREFUeJzt3XGs3tV93/H3B1wTkjJTK5MxYIj/uEh1RUTird6WZkmE\ni9xqMkxD4FRFVudFkdyl0aROxf2j2V9OyLQRpApULU5tvJRiBcW4KqK4LO3SP4gTZhoXx8NMuZN9\nF19ndmbCojE7fPfHc+54Yl37Xt9703ue6/dLurrn+f7O+fl3dCQ+nN/ze56bqkKSpN5cs9gXIEnS\ndAwoSVKXDChJUpcMKElSlwwoSVKXDChJUpdmFVBJbkzylSTfSXI0yYYkK5McTPJakheS3DjUf0eS\n40mOJblnqL4+yZF27LGh+nVJnm71l5LcvrDTlCSNmtnuoB4DnquqnwfeDxwDHgYOVtUdwIvtNUnW\nAQ8C64BNwONJ0s7zBLCtqsaAsSSbWn0bcKbVHwUemffMJEkjbcaASrIC+HBVfQmgqi5U1TlgM7Cn\nddsD3Nfa9wJPVdX5qhoHXgc2JFkN3FBVh1q/J4fGDJ/rGeDuec1KkjTyZrODWgt8P8kfJvkvSf5D\nkvcAq6pqsvWZBFa19s3AyaHxJ4FbpqlPtDrt9wkYBCBwLsnKuUxIkrQ0zCaglgEfBB6vqg8C/5t2\nO29KDb4v6af6nUm7d++e+jf88ccff/wZzZ8rMpuAOgmcrKpvttdfYRBYp5LcBNBu351uxyeANUPj\nb23nmGjti+tTY25r51oGrKiqs8MXMT4+PrsZSZKWhBkDqqpOASeS3NFKG4FXgT8BtrbaVmB/ax8A\ntiRZnmQtMAYcaud5oz0BGOAh4NmhMVPnup/BQxeSpKvYsln2+xTw5STLgf8G/AZwLbAvyTZgHHgA\noKqOJtkHHAUuANvrna9M3w7sBq5n8FTg862+C9ib5DhwBtgyz3lJkkbcrAKqqv4a+PvTHNp4if47\ngZ3T1F8G7pym/hYt4CRJAr9JQpLUKQNKktQlA0qS1CUDSpLUJQNKktQlA0qS1CUDSpLUJQNKktQl\nA0qS1CUDSpLUJQNKktQlA0qS1CUDSpLUJQNKktQlA0qS1CUDSpLUJQNKktQlA0qS1CUDSpLUJQNK\nktQlA0qS1CUDSpLUJQNKktQlA0qS1CUDSpLUJQNKktQlA0qS1CUDSpLUJQNKktQlA0qS1CUDSpLU\nJQNKktSlWQVUkvEk305yOMmhVluZ5GCS15K8kOTGof47khxPcizJPUP19UmOtGOPDdWvS/J0q7+U\n5Pa5TujHb7/N//3xhTn/nH/7x3P9pyVJC2jZLPsV8NGqOjtUexg4WFWfT/I77fXDSdYBDwLrgFuA\nP08yVlUFPAFsq6pDSZ5Lsqmqnge2AWeqaizJg8AjwJa5TmrHN5+d61D++R3/kF9YefOcx0uSFsZs\nAwogF73eDHyktfcAf8EgpO4Fnqqq88B4kteBDUn+O3BDVR1qY54E7gOeb+f6TKs/A/z+Fc7jJ7x5\n/q05j71Qb8/nn5YkLZDZvgdVDHZC30ryiVZbVVWTrT0JrGrtm4GTQ2NPMthJXVyfaHXa7xMAVXUB\nOJdk5ZVMRJK0tMx2B/Whqvpekr8LHExybPhgVVWSWvjLkyRdrWa1g6qq77Xf3we+CvwiMJnkJoAk\nq4HTrfsEsGZo+K0Mdk4TrX1xfWrMbe1cy4AVF73fJUm6yswYUEneneSG1n4PcA9wBDgAbG3dtgL7\nW/sAsCXJ8iRrgTHgUFWdAt5IsiFJgIeAZ4fGTJ3rfuDFec9MkjTSZnOLbxXw1UGmsAz4clW9kORb\nwL4k24Bx4AGAqjqaZB9wFLgAbG9P8AFsB3YD1wPPtSf4AHYBe5McB84wjyf4JElLw4wBVVXfBe6a\npn4W2HiJMTuBndPUXwbunKb+Fi3gJEkCv0lCktQpA0qS1CUDSpLUJQNKktQlA0qS1CUDSpLUJQNK\nktQlA0qS1CUDSpLUJQNKktQlA0qS1CUDSpLUJQNKktQlA0qS1CUDSpLUJQNKktQlA0qS1CUDSpLU\nJQNKktQlA0qS1CUDSpLUJQNKktQlA0qS1CUDSpLUJQNKktQlA0qS1CUDSpLUJQNKktQlA0qS1CUD\nSpLUJQNKktQlA0qS1CUDSpLUpVkFVJJrkxxO8ift9cokB5O8luSFJDcO9d2R5HiSY0nuGaqvT3Kk\nHXtsqH5dkqdb/aUkty/kBCVJo2m2O6hPA0eBaq8fBg5W1R3Ai+01SdYBDwLrgE3A40nSxjwBbKuq\nMWAsyaZW3wacafVHgUfmNyVJ0lIwY0AluRX4VeCLwFTYbAb2tPYe4L7Wvhd4qqrOV9U48DqwIclq\n4IaqOtT6PTk0ZvhczwB3z3k2kqQlYzY7qEeBfw28PVRbVVWTrT0JrGrtm4GTQ/1OArdMU59oddrv\nEwBVdQE4l2TlFcxBkrQEXTagkvwT4HRVHead3dNPqKrinVt/kiQtiGUzHP9HwOYkvwq8C/g7SfYC\nk0luqqpT7fbd6dZ/AlgzNP5WBjunida+uD415jbgfyRZBqyoqrPzmZQkafRddgdVVb9bVWuqai2w\nBfhPVfUQcADY2rptBfa39gFgS5LlSdYCY8ChqjoFvJFkQ3to4iHg2aExU+e6n8FDF5Kkq9xMO6iL\nTd3K+xywL8k2YBx4AKCqjibZx+CJvwvA9nYLEGA7sBu4Hniuqp5v9V3A3iTHgTMMglCSdJWbdUBV\n1V8Cf9naZ4GNl+i3E9g5Tf1l4M5p6m/RAk6SpCl+k4QkqUsGlCSpSwaUJKlLBpQkqUsGlCSpSwaU\nJKlLBpQkqUsGlCSpSwaUJKlLBpQkqUsGlCSpSwaUJKlLBpQkqUsGlCSpSwaUJKlLBpQkqUsGlCSp\nSwaUJKlLBpQkqUsGlCSpSwaUJKlLBpQkqUsGlCSpSwaUJKlLBpQkqUsGlCSpSwaUJKlLBpQkqUsG\nlCSpSwaUJKlLBpQkqUsGlCSpS5cNqCTvSvKNJK8kOZrks62+MsnBJK8leSHJjUNjdiQ5nuRYknuG\n6uuTHGnHHhuqX5fk6VZ/KcntP42JSpJGy2UDqqr+D/CxqroLeD/wsSS/BDwMHKyqO4AX22uSrAMe\nBNYBm4DHk6Sd7glgW1WNAWNJNrX6NuBMqz8KPLKQE5QkjaYZb/FV1Y9aczlwLfADYDOwp9X3APe1\n9r3AU1V1vqrGgdeBDUlWAzdU1aHW78mhMcPnega4e86zkSQtGTMGVJJrkrwCTAJfq6pXgVVVNdm6\nTAKrWvtm4OTQ8JPALdPUJ1qd9vsEQFVdAM4lWTm36UiSloplM3WoqreBu5KsAP4syccuOl5J6qd1\ngZKkq9Osn+KrqnPAnwLrgckkNwG023enW7cJYM3QsFsZ7JwmWvvi+tSY29q5lgErqursFc9EkrSk\nzPQU33unntBLcj3wy8Bh4ACwtXXbCuxv7QPAliTLk6wFxoBDVXUKeCPJhvbQxEPAs0Njps51P4OH\nLiRJV7mZbvGtBvYkuYZBmO2tqheTHAb2JdkGjAMPAFTV0ST7gKPABWB7VU3d/tsO7AauB56rqudb\nfRewN8lx4AywZaEmJ0kaXZcNqKo6AnxwmvpZYOMlxuwEdk5Tfxm4c5r6W7SAkyRpit8kIUnqkgEl\nSeqSASVJ6pIBJUnqkgElSeqSASVJ6pIBJUnqkgElSeqSASVJ6pIBJUnqkgElSeqSASVJ6pIBJUnq\nkgElSeqSASVJ6pIBJUnqkgElSerSTH/y/ar1ya//0ZzH/sGHf20Br0SSrk7uoCRJXTKgJEldMqAk\nSV0yoCRJXTKgJEldMqAkSV0yoCRJXTKgJEldMqAkSV0yoCRJXTKgJEldMqAkSV0yoCRJXTKgJEld\nmjGgkqxJ8rUkryb5myS/1eorkxxM8lqSF5LcODRmR5LjSY4luWeovj7JkXbssaH6dUmebvWXkty+\n0BOVJI2W2eygzgP/qqp+AfgHwG8m+XngYeBgVd0BvNhek2Qd8CCwDtgEPJ4k7VxPANuqagwYS7Kp\n1bcBZ1r9UeCRBZmdJGlkzRhQVXWqql5p7TeB7wC3AJuBPa3bHuC+1r4XeKqqzlfVOPA6sCHJauCG\nqjrU+j05NGb4XM8Ad89nUpKk0XdF70EleR/wAeAbwKqqmmyHJoFVrX0zcHJo2EkGgXZxfaLVab9P\nAFTVBeBckpVXcm2SpKVl1gGV5GcZ7G4+XVU/HD5WVQXUAl+bJOkqNquASvIzDMJpb1Xtb+XJJDe1\n46uB060+AawZGn4rg53TRGtfXJ8ac1s71zJgRVWdveLZSJKWjNk8xRdgF3C0qr4wdOgAsLW1twL7\nh+pbkixPshYYAw5V1SngjSQb2jkfAp6d5lz3M3joQpJ0FVs2iz4fAn4d+HaSw622A/gcsC/JNmAc\neACgqo4m2QccBS4A29stQIDtwG7geuC5qnq+1XcBe5McB84AW+Y5L0nSiJsxoKrqr7j0TmvjJcbs\nBHZOU38ZuHOa+lu0gJMkCfwmCUlSpwwoSVKXDChJUpcMKElSlwwoSVKXDChJUpcMKElSlwwoSVKX\nDChJUpcMKElSlwwoSVKXDChJUpcMKElSlwwoSVKXDChJUpcMKElSlwwoSVKXDChJUpcMKElSlwwo\nSVKXDChJUpcMKElSlwwoSVKXDChJUpcMKElSlwwoSVKXDChJUpcMKElSlwwoSVKXDChJUpcMKElS\nlwwoSVKXZgyoJF9KMpnkyFBtZZKDSV5L8kKSG4eO7UhyPMmxJPcM1dcnOdKOPTZUvy7J063+UpLb\nF3KCkqTRNJsd1B8Cmy6qPQwcrKo7gBfba5KsAx4E1rUxjydJG/MEsK2qxoCxJFPn3AacafVHgUfm\nMR9J0hIxY0BV1deBH1xU3gzsae09wH2tfS/wVFWdr6px4HVgQ5LVwA1Vdaj1e3JozPC5ngHunsM8\nJElLzFzfg1pVVZOtPQmsau2bgZND/U4Ct0xTn2h12u8TAFV1ATiXZOUcr0uStETM+yGJqiqgFuBa\nJEn6/+YaUJNJbgJot+9Ot/oEsGao360Mdk4TrX1xfWrMbe1cy4AVVXV2jtclSVoi5hpQB4Ctrb0V\n2D9U35JkeZK1wBhwqKpOAW8k2dAemngIeHaac93P4KELSdJVbtlMHZI8BXwEeG+SE8DvAZ8D9iXZ\nBowDDwBU1dEk+4CjwAVge7sFCLAd2A1cDzxXVc+3+i5gb5LjwBlgy8JMTZI0ymYMqKr6+CUObbxE\n/53AzmnqLwN3TlN/ixZwkiRN8ZskJEldMqAkSV2a8RafFs8nv/5Hcx77Bx/+tQW8Ekn62+cOSpLU\nJQNKktQlA0qS1CUDSpLUJQNKktQln+L7KZnPE3jgU3iS5A5KktQlA0qS1CUDSpLUJQNKktQlA0qS\n1CUDSpLUJQNKktQlA0qS1CUDSpLUJQNKktQlA0qS1CUDSpLUJQNKktQlA0qS1CUDSpLUJQNKktQl\nA0qS1CX/oq4uyb8KLGkxuYOSJHXJHdQSNp8dkLsfSYvNHZQkqUsGlCSpSwaUJKlLBpQkqUvdPCSR\nZBPwBeBa4ItV9cgiX5I64KPu0tWri4BKci3w+8BGYAL4ZpIDVfWdxb0yzZdPEkqaq15u8f0i8HpV\njVfVeeCPgXsX+ZokSYuoix0UcAtwYuj1SWDDXE50TcJvv3/jnC9k9btXzHmsliZ3gdLiSFUt9jWQ\n5J8Bm6rqE+31rwMbqupTQ32+yCC4JEmjabyqds+2cy87qAlgzdDrNVwURlX1L/5Wr0iStKh6eQ/q\nW8BYkvclWQ48CBxY5GuSJC2iLnZQVXUhyb8E/ozBY+a7fIJPkq5uXbwHJUnSxXq5xXdJSTYlOZbk\neJLfWezrWQhJxpN8O8nhJIcW+3rmIsmXkkwmOTJUW5nkYJLXkryQ5MbFvMYrdYk5/ZskJ9taHW4f\nKB8ZSdYk+VqSV5P8TZLfavWRXavLzGlk1yrJu5J8I8krSY4m+Wyrj+Q6XWY+V7RGXe+g2gd4/ytD\nH+AFPj7qt/+SfBdYX1VnF/ta5irJh4E3gSer6s5W+zzwP6vq8+1/Jn6uqh5ezOu8EpeY02eAH1bV\nv1/Ui5ujJDcBN1XVK0l+FngZuA/4DUZ0rS4zpwcY7bV6d1X9KMky4K+A3wY2M7rrNN187uYK1qj3\nHdRS/gBvFvsC5qOqvg784KLyZmBPa+9h8B+NkXGJOcEIr1VVnaqqV1r7TeA7DD53OLJrdZk5wWiv\n1Y9aczmD9+J/wGiv03TzgStYo94DaroP8N5yib6jpIA/T/KtJJ9Y7ItZQKuqarK1J4FVi3kxC+hT\nSf46ya5RucUynSTvAz4AfIMlslZDc3qplUZ2rZJck+QVBuvxtap6lRFep0vMB65gjXoPqH7vP87P\nh6rqA8CvAL/Zbi0tKTW4d7wU1u8JYC1wF/A94N8t7uXMTbsV9gzw6ar64fCxUV2rNqevMJjTm4z4\nWlXV21V1F3Ar8I+TfOyi4yO1TtPM56Nc4Rr1HlAzfoB3FFXV99rv7wNfZXArcymYbO8PkGQ1cHqR\nr2fequp0NcAXGcG1SvIzDMJpb1Xtb+WRXquhOf3HqTkthbUCqKpzwJ8C6xnxdYKfmM/fu9I16j2g\nltwHeJO8O8kNrf0e4B7gyOVHjYwDwNbW3grsv0zfkdD+ozDlnzJia5UkwC7gaFV9YejQyK7VpeY0\nymuV5L1Tt7uSXA/8MnCYEV2nS81nKmybGdeo66f4AJL8Cu/8nahdVfXZRb6keUmylsGuCQYflP7y\nKM4pyVPAR4D3MrjH/HvAs8A+4DZgHHigqv7XYl3jlZpmTp8BPsrgdkQB3wU+OfSeQPeS/BLwn4Fv\n887toR3AIUZ0rS4xp98FPs6IrlWSOxk8BHFN+9lbVf82yUpGcJ0uM58nuYI16j6gJElXp95v8UmS\nrlIGlCSpSwaUJKlLBpQkqUsGlCSpSwaUJKlLBpQkqUv/D816spbHcMs9AAAAAElFTkSuQmCC\n",
       "text": [
        "<matplotlib.figure.Figure at 0x1077213d0>"
       ]
      }
     ],
     "prompt_number": 4
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Problem 1 (b)"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from sklearn.cross_validation import Bootstrap\n",
      "\n",
      "Ns = [100, 1000, 10000]\n",
      "B = 1000\n",
      "\n",
      "true_p = 0.3\n",
      "\n",
      "fig1, axes1 = plt.subplots(ncols=3, sharex=True, figsize=(12, 4))\n",
      "fig2, axes2 = plt.subplots(ncols=3, sharex=True, figsize=(12, 4))\n",
      "bins = np.arange(0.2, 0.4, 0.005)\n",
      "\n",
      "for n, ax1, ax2 in zip(Ns, axes1.flat, axes2.flat):\n",
      "    # Generate random data\n",
      "    data = np.random.geometric(p=true_p, size=n)\n",
      "    \n",
      "    # Do parametric bootstrap\n",
      "    p_mle = 1/data.mean()\n",
      "    p_bootstrap_parametric = np.zeros(B)\n",
      "    data_bootstrap = np.zeros(n)\n",
      "    for i in range(B):\n",
      "        data_bootstrap = np.random.geometric(p=p_mle, size=n)\n",
      "        p_bootstrap_parametric[i] = 1/data_bootstrap.mean()\n",
      "        \n",
      "    # Plot parametric bootstrap\n",
      "    ppl.hist(p_bootstrap_parametric, ax=ax1, bins=bins)\n",
      "    ymin, ymax = ax1.get_ylim()\n",
      "\n",
      "    # vertical line indicating MLE\n",
      "    vline = ax1.vlines(p_bootstrap_parametric.mean(), 0, ymax)\n",
      "#     print p_mle, ymax\n",
      "#     std_error = np.square(np.abs(p_bootstrap - p_bootstrap.mean())).sum()\n",
      "    std_error = p_bootstrap_parametric.std()\n",
      "    print 'Parametric: \\hat{{p}} = {:.5f}\\tstderr = {:.5f}\\t1/\\sqrt{{N}} = {:.5f}, stderr=(1/sqrt(N))*{:.5f}'.format(p_bootstrap_parametric.mean(), \n",
      "                                                                                  std_error, 1/np.sqrt(n), std_error/(1/np.sqrt(n)))\n",
      "    ax1.annotate('$\\hat{{p}} = {:.4f}$'.format(p_mle), (p_mle, ymax), fontsize=14, \n",
      "                 xytext=(10, -15), textcoords='offset points')\n",
      "    ax1.set_title('Parametric, $N = {}$'.format(n))\n",
      "    ax1.set_xlabel('$p^*$')\n",
      "        \n",
      "        \n",
      "    # Do the non-parametric bootstrap\n",
      "    p_bootstrap_non_parametric = np.zeros(B)\n",
      "    bs = Bootstrap(n=n, n_iter=B)\n",
      "    data_bootstrap = np.zeros(n)\n",
      "    for i, (train_index, test_index) in enumerate(bs):\n",
      "        data_bootstrap[:(n/2)] = data[train_index]\n",
      "        data_bootstrap[(n/2):] = data[test_index]\n",
      "        p_bootstrap_non_parametric[i] = 1/data_bootstrap.mean()\n",
      "        \n",
      "    # Plot non-parametric bootstrap\n",
      "    ppl.hist(p_bootstrap_non_parametric, ax=ax2, bins=bins)\n",
      "    ymin, ymax = ax2.get_ylim()\n",
      "\n",
      "    vline = ax2.vlines(p_mle, 0, ymax)\n",
      "#     print p_mle, ymax\n",
      "#     std_error = np.square(np.abs(p_bootstrap - p_bootstrap.mean())).sum()\n",
      "    std_error = p_bootstrap_non_parametric.std()\n",
      "    print 'Non-parametric: \\hat{{p}} = {:.5f}\\tstderr = {:.5f}\\t1/\\sqrt{{N}} = {:.5f}, stderr=(1/sqrt(N))*{:.5f}'.format(p_bootstrap_non_parametric.mean(), \n",
      "                                                                                  std_error, 1/np.sqrt(n), std_error/(1/np.sqrt(n)))\n",
      "    ax2.annotate('$\\hat{{p}} = {:.4f}$'.format(p_mle), (p_mle, ymax), fontsize=14, \n",
      "                 xytext=(10, -15), textcoords='offset points')\n",
      "    ax2.set_title('Non-parametric $N = {}$'.format(n))\n",
      "    ax2.set_xlabel('$p^*$')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Parametric: \\hat{p} = 0.35897\tstderr = 0.02880\t1/\\sqrt{N} = 0.10000, stderr=(1/sqrt(N))*0.28800\n",
        "Non-parametric: \\hat{p} = 0.35866\tstderr = 0.02493\t1/\\sqrt{N} = 0.10000, stderr=(1/sqrt(N))*0.24932"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "Parametric: \\hat{p} = 0.31611\tstderr = 0.00813\t1/\\sqrt{N} = 0.03162, stderr=(1/sqrt(N))*0.25704"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "Non-parametric: \\hat{p} = 0.31568\tstderr = 0.00811\t1/\\sqrt{N} = 0.03162, stderr=(1/sqrt(N))*0.25649"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "Parametric: \\hat{p} = 0.30440\tstderr = 0.00252\t1/\\sqrt{N} = 0.01000, stderr=(1/sqrt(N))*0.25212"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "Non-parametric: \\hat{p} = 0.30439\tstderr = 0.00254\t1/\\sqrt{N} = 0.01000, stderr=(1/sqrt(N))*0.25436"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAA1gAAAEYCAYAAABBWFftAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3X2YVeV56P/vrQMawQhIRUREzU8jiEpQSashaYx6mVgR\nrcGXaLyU5NgkVdMkJpqeRpIeqTHqSdLEVo1JKFUqMWLUxp8iVXw5BWt8F60YMydA5K2oRaIG9D5/\n7DW4GWaGYWYxe8/s7+e69uVaz3rW2s+azb5d936e9azITCRJkiRJ3bddrRsgSZIkSX2FCZYkSZIk\nlcQES5IkSZJKYoIlSZIkSSUxwZIkSZKkkphgSZIkSVJJTLAkSZIkqSQmWJIkSZJUEhMslS4inomI\nD9e6HZL6JmOMpG3NOKPuiMysdRsaUkQ0A7sBbwPrgLuAv8zMdbVsV0eKNp+bmf9W67Z0V0SMA87M\nzK9UlU0GxgDvAMsyc2ZH5VI9M8bUVhkxxtijemecqa1tGWeMP92Umb5q8AJ+AxxVLO8BPA383VYe\no6kGbf5YT7YH+CAwB1jacnxgGPAvwJ3AEV045peAW4GfVJXtAvyqav3fgV3bKR9a638/vnxt6WWM\n6fR71mOMMfb46hUv40yn37M3xRnjTwkvhwjWgcz8HfD/A2MBIuLiiHgxIv47Ip4tfkWg2NYcEV+N\niKeAtRGxfSfqfyUinoqItRFxQ0QMi4i7IuK1iJgbEYOq6u8RET+PiJUR8VJEnF+UzwT2Au4ojvOV\nDtrTHBEfK7aPjIhbi+Otjoi/38q/zcLib/MC8OdF2QoqAemTmfl/uvD3vhr4RaviDwOLqtafBI5q\np/yjW/ueUi0ZYzr829RjjDH2qNcxznT4t+lNccb4U4KmWjegwQVUvrjAx4GfF+UvAh/KzOURMQX4\n54h4X/FlBDitqL86M9+OiI7qJ3Ay8DGgH/A48AHgHOB54JfABcC3ImI74A4qv7KcCowE7o2I/8zM\nsyLiQ8DU3LxbvXV7EsiI2J5K8LgX+BSVbubDtuoPVGnTeuD7wEXAzcWmAZn5RlFnX+CzHRxmQWa2\nDkLRan1P4NWq9VeB/YA17ZRLvYExZkt/oPqNMcYe9RbGmS39gXpfnDH+dJMJVu0EcFtEbABeo/Ll\nnQ6Qmbe0VMrM2RFxCTCBSsBI4PuZuayqTkf1Af4+M1cBRMSDwIrMfLJYn0MlYAEcTqUL+H8V67+J\niB9RCTr3tHMem7WnygRgOHBRZr5TlD3c8Z9lM+OBR4FngKsjYnxmPla8b8s5vwRcspXHbX3z4SDg\nzar1PwADi3ptlUv1zhjTOfUaY4w96g2MM53T2+KM8aebHCJYOwmcmJmDM3PvzPzLzHwLICI+HRGP\nR8QrEfEKle72oVX7Lqk+UCfqr6hafqPV+pu8+6UZBezRcpziWJdQuYG1I0vaKR8J/N+qgNQVB2fm\nU8UxrgHOj4j3A//ZjWPC5r/6rG1V9h4qv+C0Vy7VO2NM59RrjDH2qDcwznROb4szxp9usgerzkTE\nKOA6KmNg/z0zMyIeZ9N/6LmV9Td7m3bKlwC/ycz929ne3pST7ZX/FtgrIrbPzLc7aE9Hqn8E+BGV\nIQeLgO+1FHaxW711m3/Npl3+uwKPUekWry4fWpRLvZIxZjP1FmOMPer1jDOb6S1xxvhTEhOs+jOA\nyhdmNbBdRHya4obRkup35BEqN3d+Ffh7Kl3Co4EdM/NRKr8WvQ/o7NSmjwAvA5dHxKVUxi2Pb7mZ\nMyJ+CmRmntPWzhHRr2gDVCq+GhG3AB/NzO9UlXelW711YH4AuKJq/VDgYirTzlaXjwe+tpXvJdUT\nY0yhTmOMsUd9gXGm0MvijPGnJA4RrDOZuQi4isqUmMupBJiHyqrfslur5SyO9TbwZ8A44CVgFZVf\nlN5b1P074H8WXe5f6sS5vAOcAPx/VH4BWgJMqaqyZ3ttjYjDqdwEemxEjKja9H3gwS29d0ci4i+B\nc4E/jYhLI+K9WXlmxxUR8T8j4hvAFZm5sr3y7ry/VEvGmIp6jzHGHvVmxpmK3hpnjD/d1+GDhiPi\nx8DxwMrMPKgoG0LlH8sooBmYkpmvFtsuofJhvw1ckJnt3UyoBhcR/anMAnRwN7rc1YdEZQaqf6Iy\nTj6B6zLz+xExDfgMlf9JAnw9M+8q9jHmqE3GGG1JcQ/Mv1QV7Qv8DfDPeJ2jTjDOqD1bSrAmAq8D\n/1SVYF1BZQrLKyLia8DgzLw4IsYAN1GZvWUEleks9+/mTYGSGkRE7A7snplPRMRA4FfAZCq/FK7N\nyjM/qusbcySVIirTaC+jMmPc+XidI6kbOhwimJkPAq+0Kp4EzCiWZ1C5AAI4EZiVmeszs5nKDXwT\nymuqpL4sM5dn5hPF8uvAc1QuYqDtm5mNOZLKcjTwYmYuwescSd3UlXuwhuW7D4lbAQwrlvcAllbV\nW8q7F0eS1GkRsTeVh0guKIrOj4gnI+KGiBhUlBlzJJXlNGBWsex1jqRu6dYkF1kZX9j+GMM2tv30\npz9t2ceXr4Z7TZs2reZtKOm1zRTDA28BLix6sv4B2IfKDcsvU7kRuj2btc2Y42trXn3oO9qXXttU\ncR/NCcDPWm/zOseXr4Z8dVtXEqwVxb0SRMRwoGVWkWVUHsbWYs+ibBPNzc1deEupb/jmN79Z6ybU\ntWI6258D/5yZtwEUMxplcaHzI94dkmPMUen8jjakjwO/ysyWiXS8zpHULV1JsG4Hzi6WzwZuqyo/\nLSL6R8Q+wH5Unh0gSVsUEQHcACzKzO9WlQ+vqnYS8HSxbMyRVIbTeXd4IHidI6mbOnzQcETMAj4C\nDI2IJcA3gMuB2RExlWL6UoDMXBQRs6k8mXoD8PnsaIpCSdrUkcCZwFMR8XhR9nXg9IgYR6Xb/jfA\neWDMkdR9ETGAygQXn60q9jpHUrd0mGBl5untbDq6nfrTgendbZSkxpOZD9F2r/pdHexjzJHUZcUD\nVYe2KluD1zmSuqFbk1xI2rZWrVq15UqSas7vqiSphQmWVKd+9rOf8ZWvfIVbb721lOM9+uijXHjh\nhcycOZO/+Iu/4Ne//nW7de+//35mzZrFDTfcwJlnnsm8efM2bjv33HOJiPUR8UZEPBIR41u2RcTS\niHil1euHVdt3johbImIkUh9Rr99VgLVr13LKKaewZMmSTcrPPfdc+vXrx3ve8x4mTJjAY489tnHb\nnnvuyeDBgzd5feELXyjl3CSpEXQ4RFBSbbz00kvsvPPOzJgxg3/913+lubmZvffeu8vHe+uttzjl\nlFNYuHAhw4YNY/To0Zx++uk88kjb92d/8pOf5KqrrmLq1KkMGjSISZMmsXLlSgYMGMCoUaOg8uyX\nqHpWDBExDLiCyo3gLfclfAX4m2L7VCqzbp0MfKnLJyNtY+c9eNMm69dOPKPduvX8Xb3hhhtYunQp\nt956K1dfffUm+40aNYply5aRmQwbNmxj+YoVK/jqV7/KpEmTqMw7A1deeSV/+7d/2+VzkqRGY4Il\n1aF9992XfffdF4Djjz++28d74IEHGDhw4MYLqUMPPZTnnnuu3YvB+fPnbyx/55132LBhwybbM3Pl\nZjtVzMzMVwAiYgpwU2a+WuxzQ1F+abdPSKoT9fxdnTp1KtD+1PO77bZbm+VnnXUWgwcPBmD27Nmc\nccYZDBo0qM26kqTNmWBJdWbdunVcc801LFiwgHPPPZc1a9bw2GOPccIJJ3DUUUcB8Morr/Cd73yH\njiawampq4tJLL6WpqYnm5mZ23XXXjdsigsGDB/Pss8+2edE2ZsyYjctz5sxh2rRpDBgwAIA33niD\niPgfwH8DHwOuzsznWvVmjQDGZObsbv0xpDpW79/Vjrzxxhtcd911vPe972XevHl86UtfYvTo0Zv0\nZi1btoxFixYxZcqUzvw5JEkFEyypzsyZM4fPfe5z/OIXv+DNN9/krLPOYsKECUyZMoUnn3wSgMGD\nBzN9eucnslq9ejU77bTTJmU77rgja9eubXefJ554gnnz5jFw4EC++MUvbiw/6KCDAH6SmesjYiVw\nW0Qc0Gq64suAaZ1uoNQL1ft3tSMHHXQQU6ZMoV+/fuy2225MnjyZ559/fuOwQIC//uu/Ztq0aZ1u\nuySpwkkupDpzwgkn0K9fPxYvXsxJJ50EwG9/+1tWr17d5WMOGjRos1/QX3/9dYYOHdrOHjBu3Di+\n/OUvc/jhhzNx4kTWrVsHwGmnnUZmri+qvUTlYZsHtewXEbsBH83M5i43WOoF6v272pHTTjuNfv36\nAZVhjosXL+bpp5/euH3lypXcd9993bqfTJIalT1YUo20vpEeKjfT77LLLjzwwANMmDCB7bar/AZy\n9913c+yxx26st2bNGq688soOhx1tv/32TJs2jaamJg444ACuvfbajds2bNjAmjVrWias2MSCBQuY\nPHkyCxcuZNSoUUycOJHzzjuPu+++mz322INjjjmG119/fYfMfAvYudjtD1WH+DjwX1vzt5B6o3r+\nrp588sntvt+CBQs45phjWL16NTvssMPG3rH+/ftvrHPXXXdtMlRRktR5JlhSHbr//vs5+OCDgcrz\nde644w7mzp27cfuQIUO2atjRxIkTWbVqFUuWLGHkyJHMnz+fAw88kP322w+AefPmMXToUA455BCa\nmpoYO3Ysw4cPByqzpPXv359x48axww47cNFFF/GNb3zjreLQRwIPZ+bzVW83Fvj9FpoUW9gu9Qr1\n+l1trTrBGzlyJBdddBE77LADAA8//DBHHnkkBxxwwMY6zzzzzGZDFSVJnWOCJdWh+fPnc8QRRzBr\n1iwWLlzILbfcwl577dXl4zU1NTFz5kwuu+wyjjjiCO677z5uvvnmjdt/+MMfMn78eA455BAOO+ww\nzjnnHH7wgx+w3Xbb8dBDD3HnnXdunClt/PjxRMSXge2pDA88qdXbvQb8Z+s2RMQZwIeoTOF+eUQ8\nlJk/bF1P6k3q+bt600038dBDDxERXHzxxXzoQx/iC1/4AiNGjGD8+PFcddVVvP322yxevJg5c+Zs\n0o5ddtmF97///V0+D0lqZNHRsIVtYdq0aelNs2pUEbHxl+T2hgj+4Q9/YOTIkSxfvnyTG87rTN02\nrDVjjrZGRPA/Hrhxk7L2noPVS76rfUGv+uMac6Rer9sxx0kupDqzcOFCxo4d6wWbVOf8rkqS2mKC\nJdWRZ555hm9961usXr16k/s4JNUXv6uSpPZ4D5ZUR8aOHevFmtQL+F2VJLXHHixJkiRJKokJliRJ\nkiSVxARLkiRJkkpigiVJkiRJJTHBkiRJkqSSmGBJkiRJUklMsCRJkiSpJCZYkiRJklQSEyxJktSQ\nImJQRNwSEc9FxKKI+GBEDImIuRHxQkTcExGDqupfEhGLI+L5iDi2lm2XVL9MsCRJUqP6HvDLzBwN\nHAw8D1wMzM3M/YF5xToRMQY4FRgDHAdcExFeR0najIFBkiQ1nIjYBZiYmT8GyMwNmfkaMAmYUVSb\nAUwulk8EZmXm+sxsBl4EJvRsqyX1BiZYkiSpEe0DrIqIn0TEYxFxfUQMAIZl5oqizgpgWLG8B7C0\nav+lwIiea66k3sIES5IkNaImYDxwTWaOB9ZRDAdskZkJZAfH6GibGsi0adNq3QTVERMsSZLUiJYC\nSzPzP4r1W6gkXMsjYneAiBgOrCy2LwNGVu2/Z1Em8c1vfrPWTVAdaap1AyRJ6mnnPXjTZmXXTjyj\nBi1RrWTm8ohYEhH7Z+YLwNHAs8XrbODbxX9vK3a5HbgpIq6mMjRwP+CRnm+5pHpngiVJkhrV+cCN\nEdEf+DVwDrA9MDsipgLNwBSAzFwUEbOBRcAG4PPFEEJJ2oQJliRJakiZ+SRweBubjm6n/nRg+jZt\nlKRez3uwJEmSpBpYtWpVrZugbcAeLEmSJKmH/exnP+POO+/kxBNP5OSTT+728R599FFmzpzJYYcd\nxsMPP8xFF13E+973vjbr3n///bz88sv8/ve/57777uOcc87hYx/72FYd58tf/jInn3wyRx555BaP\n2WjswZIkSZJ60EsvvcTOO+/MjBkz2GGHHWhubu7W8d566y1OOeUUvv71r3PWWWfxmc98htNPP73d\n+p/85CdZv349U6dO5aSTTmLSpEmsW7eu08eZP38+N954I2+//fYWj9mI7MGSJEmSetC+++7Lvvvu\nC8Dxxx/f7eM98MADDBw4kGHDKs/FPvTQQ3nuuedobm5m77333qz+/PnzN5a/8847bNiwodPHee21\n13jqqacYPXp0p47ZiLqcYEXEJcCZwDvA01Rm3hkA3AyMoph5JzNf7X4zJfV1ETES+CdgNyoP77wu\nM78fEUNoJ64Ucehc4G3ggsy8pxZtlySps9atW8c111zDggULOPfcc1mzZg2PPfYYJ5xwAkcddRQA\nr7zyCt/5znfoaKLKpqYmLr30UpqammhubmbXXXfduC0iGDx4MM8++2ybCdaYMWM2Ls+ZM4dp06Yx\nYMCATh3n+uuv54ILLuDWW2/t1DEbUZcSrIjYG/gsMDoz34qIm4HTgAOBuZl5RUR8jcoT0S9u90CS\n9K71wF9l5hMRMRD4VUTMpfLjzWZxJSLGAKcCY6g8k+be4nk279TqBCRJ2pI5c+bwuc99jl/84he8\n+eabnHXWWUyYMIEpU6bw5JNPAjB48GCmT+/8hJWrV69mp5122qRsxx13ZO3ate3u88QTTzBv3jwG\nDhzIF7/4xU4d54477uATn/gE/fv37/QxG1FXe7D+m8rF0E4R8TawE/A74BLgI0WdGcD9mGBJ6oTM\nXA4sL5Zfj4jnqCROk2g7rpwIzMrM9UBzRLwITAAW9HDTJUkNpK0HlW+NE044gX79+rF48WJOOukk\nAH7729+yevXqLh9z0KBBm/V2vf766wwdOrTdfcaNG8e4ceO4/vrrmThxIvPnz2eXXXZp9zi/+93v\nePXVVzfpqWpdt61jNmIvVpcSrMxcExFXAb8F3gDuzsy5ETEsM1cU1VYAw0pqp6QGUvSSfwBYCLQX\nV/Zg02RqKZWETJKkurXLLrvwwAMPMGHCBLbbrjLf3N13382xxx67sc6aNWu48sorOxwiuP322zNt\n2jSampo44IADuPbaazdu27BhA2vWrGHUqFGb7bdgwQImT57MwoULGTVqFBMnTuS8887j7rvvZvTo\n0Vx33XVtHueee+5h+fLlfPvb3wbghRdeYNasWaxbt44hQ4a0e8wyZkjsbbo6RPB9wBeBvYHXgJ9F\nxJnVdTIzI8InnEvaKsXwwJ8DF2bm2ojYuK0TccWYI0mqe/fffz8HH3wwUHkW1h133MHcuXM3bh8y\nZMhWDRGcOHEiq1atYsmSJYwcOZL58+dz4IEHst9++wEwb948hg4dyiGHHEJTUxNjx45l+PDhQGVG\nw/79+zNu3Dj22muvdo/TcqwW//iP/8gZZ5zBhz/8YR599NF2j9mIujpE8DDg/2TmfwFExK3AnwDL\nI2L3zFweEcOBlSW1U1IDiIh+VJKrmZl5W1G8op24sgwYWbX7nkWZJEl1bf78+RxxxBHMmjWLhQsX\ncsstt7DXXnt1+XhNTU3MnDmTyy67jCOOOIL77ruPm2++eeP2H/7wh4wfP55DDjmEww47jHPOOYcf\n/OAHbLfddjz00EPceeedG2c17Og4AMuWLeN73/sey5cv56qrrmLt2rUcf/zxHR6z0XQ1wXoe+JuI\neA/wJnA08AiwDjgb+Hbx39vaPYIkVYlKV9UNwKLM/G7VpttpO67cDtwUEVdTGRq4H5U4JElS3frD\nH/7AM888w7333ktEdPi8qq1x1FFHbZyF8NOf/vQm21rP+PepT31q43LrySg6Og7AiBEjuOKKK7ji\niis6fcxG09V7sJ6MiH8CHqUyTftjwHXAzsDsiJhKMZ1ySe2UelRbN7BeO/GMGrSkoRxJ5dEPT0XE\n40XZJcDltBFXMnNRRMwGFgEbgM9nR4PVJUmqAwsXLmTs2LFUD4FX39Ll52Bl5hXAFa2K11DpzZKk\nrZKZDwHbtbO5zbiSmdOBzg9SlySphp555hm+9a1vsXr1aubOncsxxxxT6yZpG+hygiVJkiSp88aO\nHbvJZBbqm9r7tViSJEmStJVMsCRJkiSpJCZYkiRJklQSEyxJkiRJKokJliRJkiSVxARLkiRJkkpi\ngiVJkiRJJTHBkiRJkqSSmGBJkqSGFBHNEfFURDweEY8UZUMiYm5EvBAR90TEoKr6l0TE4oh4PiKO\nrV3LJdUzEyxJktSoEvjTzPxAZk4oyi4G5mbm/sC8Yp2IGAOcCowBjgOuiQivoyRtxsAgSZIaWbRa\nnwTMKJZnAJOL5ROBWZm5PjObgReBCUhSKyZYkiSpUSVwb0Q8GhGfLcqGZeaKYnkFMKxY3gNYWrXv\nUmBEzzRTUm/SVOsGSL3JeQ/etFnZtRPPqEFLJEklODIzX46IPwLmRsTz1RszMyMiO9i/o22SGpQ9\nWJIkqSFl5svFf1cBc6gM+VsREbsDRMRwYGVRfRkwsmr3PYsySdqECZYkSWo4EbFTROxcLA8AjgWe\nBm4Hzi6qnQ3cVizfDpwWEf0jYh9gP+CRnm21pN7AIYKSJKkRDQPmRARUroduzMx7IuJRYHZETAWa\ngSkAmbkoImYDi4ANwOcz0yGCkjZjgiVJkhpOZv4GGNdG+Rrg6Hb2mQ5M38ZNk9TLOURQkiRJkkpi\ngiVJkiRJJXGIoCRJW+AjGiRJnWUPliRJkiSVxARLkiRJkkpigiVJkiRJJTHBkiRJkqSSmGBJkiRJ\nUklMsCRJkiSpJCZYkiRJklQSn4MllaCtZ+SAz8mRJElqNPZgSZIkSVJJTLAkSZIkqSQmWJIkSZJU\nEhMsSZIkSSqJCZYkSZIklaTLCVZEDIqIWyLiuYhYFBEfjIghETE3Il6IiHsiYlCZjZXUt0XEjyNi\nRUQ8XVU2LSKWRsTjxevjVdsuiYjFEfF8RBxbm1ZLkiS9qzs9WN8DfpmZo4GDgeeBi4G5mbk/MK9Y\nl6TO+glwXKuyBK7OzA8Ur7sAImIMcCowptjnmoiwV16SJNVUly5GImIXYGJm/hggMzdk5mvAJGBG\nUW0GMLmUVkpqCJn5IPBKG5uijbITgVmZuT4zm4EXgQnbsHmSJElb1NUHDe8DrIqInwCHAL8CvggM\ny8wVRZ0VwLDuN1GSOD8iPg08Cnw5M18F9gAWVNVZCoyoReNUv9p6CLgPAJckbUtdHU7TBIwHrsnM\n8cA6Wg0HzMykMrRHkrrjH6j8qDMOeBm4qoO6xhxJklRTXU2wlgJLM/M/ivVbqCRcyyNid4CIGA6s\n7H4TJTWyzFyZBeBHvDsMcBkwsqrqnkWZJElSzXQpwcrM5cCSiNi/KDoaeBa4Azi7KDsbuK3bLZTU\n0Iofa1qcBLTMMHg7cFpE9I+IfYD9gEd6un2SereI2L6YofSOYr3dGZGduVRSZ3T1HiyA84EbI6I/\n8GvgHGB7YHZETAWagSndbqGkhhERs4CPAEMjYglwKfCnETGOyvC/3wDnAWTmooiYDSwCNgCfL3q5\nJGlrXEgljuxcrLfMiHxFRHytWL+41cylI4B7I2L/zHynFo2WVL+6nGBl5pPA4W1sOrrrzZHUyDLz\n9DaKf9xB/enA9G3XIkl9WUTsCXwCuAz4UlE8icoPPVCZEfl+KknWxplLgeaIaJm5tHqyHUnqVg+W\n1Os5w5gkNbT/DVwEvLeqrL0ZkZ25VFKn+FBOSZLUcCLiz4CVmfk4bT9rrzMzIjssWdJm7MGSJEmN\n6AhgUkR8AtgReG9EzARWRMTumbm81YzIzlwqqVNMsNQQajkUsPV7n/fgTQ5DlKQay8yvA18HiIiP\nAF/JzLMi4goqMyF/m01nRL4duCkirqYyNNCZSyW1yQRLkiTp3eF+l9PGjMjOXCqps0ywJElSQ8vM\n+cD8YnkN7cyI7MylkjrDSS4kSZIkqSQmWJIkSZJUEocISpIkSd3kszXVwgRL6iUM3JIkSfXPIYKS\nJEmSVBITLEmSJEkqiQmWJEmSJJXEBEuSJEmSSmKCJUmSJEklMcGSJEmSpJKYYEmSJElSSUywJEmS\nJKkkJliSJEmSVBITLEmSJEkqiQmWJEmSJJXEBEuSJEmSSmKCJUmSJEklMcGSJEmSpJKYYEmSJElS\nSUywJEmSJKkkTbVugKRNnffgTZuVXTvxjJodR5IkSZ1nD5YkSZIklcQES5IkNZyI2DEiFkbEExGx\nKCL+rigfEhFzI+KFiLgnIgZV7XNJRCyOiOcj4tjatV5SPTPBkiRJDScz3wQ+mpnjgIOBj0bEh4CL\ngbmZuT8wr1gnIsYApwJjgOOAayLC6yhJmzEwSJKkhpSZvy8W+wPbA68Ak4AZRfkMYHKxfCIwKzPX\nZ2Yz8CIwoedaK6m3MMGSJEkNKSK2i4gngBXAfZn5LDAsM1cUVVYAw4rlPYClVbsvBUb0WGMl9RrO\nIihJkhpSZr4DjIuIXYC7I+KjrbZnRGRHh9imDZTUK9mDJUmSGlpmvgb8K3AosCIidgeIiOHAyqLa\nMmBk1W57FmWStIluJVgRsX1EPB4RdxTr7c68I0lbEhE/jogVEfF0VZkzekkqXUQMbYknEfEe4Bjg\nceB24Oyi2tnAbcXy7cBpEdE/IvYB9gMe6dlWS+oNutuDdSGwiHe7yNuceUeSOuknVGbnquaMXpK2\nheHAvxX3YC0E7sjMecDlwDER8QJwVLFOZi4CZlO57rkL+HxmOkRQ0ma6fA9WROwJfAK4DPhSUTwJ\n+EixPAO4H5MsSZ2UmQ9GxN6tituLKxtn9AKaI6JlRq8FPdJYSb1aZj4NjG+jfA1wdDv7TAemb+Om\nSerluvNr7/8GLgLeqSprb+YdSeoqZ/SSJEm9RpcSrIj4M2BlZj4ORFt1im5zu84llaYTccWYI0mS\naqqrQwSPACZFxCeAHYH3RsRMipl3MnN5q5l3JNWJ8x68abOyayee0W55HWgvrjijlyRJqjtd6sHK\nzK9n5sjM3Ac4Dfi3zDyL9mfekaSuckYvSZLUa5T1oOGWYTmXA7MjYirQDEwp6fiSGkBEzKIyocXQ\niFgCfIOJEvagAAAMhklEQVR24kpmLoqIlhm9NuCMXpIkqQ50O8HKzPnA/GK53Zl3JGlLMvP0djY5\no5ckSeoVfGaMJEmSJJXEBEuSJEmSSmKCJUmSJEklMcGSJEmSpJKYYEmSJElSScqapl3qUXX8UFxJ\nkiQ1MHuwJEmSJKkk9mBJktRF9qZLklqzB0uSJEmSSmKCJUmSJEklMcGSJEmSpJKYYEmSJElSSUyw\nJEmSJKkkJliSJEmSVBKnaZd6OaeJliRJqh/2YEmSJElSSUywJElSw4mIkRFxX0Q8GxHPRMQFRfmQ\niJgbES9ExD0RMahqn0siYnFEPB8Rx9au9ZLqmQmWJElqROuBv8rMA4E/Br4QEaOBi4G5mbk/MK9Y\nJyLGAKcCY4DjgGsiwusoSZvxHiz1Kd6PJEnqjMxcDiwvll+PiOeAEcAk4CNFtRnA/VSSrBOBWZm5\nHmiOiBeBCcCCHm66pDrnLy+SJKmhRcTewAeAhcCwzFxRbFoBDCuW9wCWVu22lEpCJkmbMMGSJEkN\nKyIGAj8HLszMtdXbMjOB7GD3jrZJalAmWJIkqSFFRD8qydXMzLytKF4REbsX24cDK4vyZcDIqt33\nLMokaRMmWJIkqeFERAA3AIsy87tVm24Hzi6WzwZuqyo/LSL6R8Q+wH7AIz3VXkm9h5NcSJKkRnQk\ncCbwVEQ8XpRdAlwOzI6IqUAzMAUgMxdFxGxgEbAB+HwxhFCSNmGCJUmSGk5mPkT7I3mObmef6cD0\nbdYoSX2CQwQlSZIkqSQmWJIkSZJUEhMsSZIkSSqJCZYkSZIklcRJLiRJvd55D960Wdm1E8+oQUsk\nSY3OHixJkiRJKokJliRJkiSVxARLkiRJkkpigiVJkiRJJTHBkiRJkqSSdCnBioiREXFfRDwbEc9E\nxAVF+ZCImBsRL0TEPRExqNzmSmpUEdEcEU9FxOMR8UhRZsyRJEl1pas9WOuBv8rMA4E/Br4QEaOB\ni4G5mbk/MK9Yl6QyJPCnmfmBzJxQlBlzJElSXelSgpWZyzPziWL5deA5YAQwCZhRVJsBTC6jkZJU\niFbrxhxJklRXuv2g4YjYG/gAsBAYlpkrik0rgGHdPb4amw8PVZUE7o2It4FrM/N6jDmSJKnOdCvB\nioiBwM+BCzNzbcS7Py5nZkZEdrN9ktTiyMx8OSL+CJgbEc9XbzTmSJKketDlWQQjoh+V5GpmZt5W\nFK+IiN2L7cOBld1voiRBZr5c/HcVMAeYgDFHkiTVma7OIhjADcCizPxu1abbgbOL5bOB21rvK0lb\nKyJ2ioidi+UBwLHA0xhzJElSnenqEMEjgTOBpyLi8aLsEuByYHZETAWagSndbqEkVe6tmlMMQ24C\nbszMeyLiUYw5kiSpjnQpwcrMh2i/9+vorjdHkjaXmb8BxrVRvgZjjiRJqiNdvgdLkiRJkrQpEyxJ\nkiRJKkm3n4MllcHnXUmSJKkvsAdLkiQ1nIj4cUSsiIinq8qGRMTciHghIu6JiEFV2y6JiMUR8XxE\nHFubVkvqDUywJElSI/oJcFyrsouBuZm5PzCvWCcixgCnAmOKfa6JCK+hJLXJ4CBJkhpOZj4IvNKq\neBIwo1ieAUwulk8EZmXm+sxsBl6k8rBzSdqMCZYkSVLFsMxcUSyvoPIMPoA9gKVV9ZYCI3qyYZJ6\nDye5UI9yMgtJUm+QmRkR2VGVHmuMpF7FBEuSpJL5Y1KvtSIids/M5RExHFhZlC8DRlbV27Mok6TN\nOERQkiSp4nbg7GL5bOC2qvLTIqJ/ROwD7Ac8UoP2SeoF7MGSJEkNJyJmAR8BhkbEEuAbwOXA7IiY\nCjQDUwAyc1FEzAYWARuAz2emQwQltckES5IkNZzMPL2dTUe3U386MH3btUhSX+EQQUmSJEkqiQmW\nJEmSJJXEBEuSJEmSSmKCJUmSJEklMcGSJEmSpJKYYEmSJElSSUywJEmSJKkkJliSJEmSVBITLEmS\nJEkqiQmWJEmSJJWkqdYNUN903oM3bVZ27cQzatASSZIkqefYgyVJkiRJJbEHS51ij5QkSZK0ZfZg\nSZIkSVJJ7MFSt9izJaknGXMkSfXOHixJkiRJKok9WJIkSVIr9pirq0ywJEnqIV6wSVLf5xBBSZIk\nSSqJCZYkSZIklcQES5IkSZJK4j1YkqS6471KkqTeqvQEKyKOA74LbA/8KDO/XfZ7SFILY46knmTM\n0dZq6wcj8EejvqzUBCsitgd+ABwNLAP+IyJuz8znynwfbTv+aqzexJgjqScZcyR1Rtn3YE0AXszM\n5sxcD/wLcGLJ7yFJLYw5knqSMUfSFpWdYI0AllStLy3KJGlbMOZI6knGHElbFJlZ3sEi/hw4LjM/\nW6yfCXwwM8+vqvMjKgFJUu/VnJk/rXUjjDlSQ6iLeAPGHKlBdDvmlD3JxTJgZNX6SFoFmcz8TMnv\nKalxGXMk9SRjjqQtKnuI4KPAfhGxd0T0B04Fbi/5PSSphTFHUk8y5kjaolJ7sDJzQ0T8JXA3lelL\nb3BmHUnbijFHUk8y5kjqjFLvwZIkSZKkRlbqEMGIOC4ino+IxRHxtTa2fyoinoyIpyLi4Yg4uLP7\n1qNunm9zUf54RDzSsy3vuk6c84nFOT8eEb+KiKM6u2896ub59rrPuLOfUUQcHhEbihu+t2rfnmyv\nMad3x5xGizdgzOmgXs1jTqPFGzDmtLG9T8WcRos30IMxJzNLeVHpKn8R2BvoBzwBjG5V50+AXYrl\n44AFnd233l7dOd9i/TfAkFqfxzY45wFVywdReV5IX/6M2zzf3vgZd/YzKur9G3An8Oe1+nyNOX07\n5jRavOnuOffVz7iqXk1jTqPFm+6ec1/999iXYk6jxZut+ZzKiDll9mBt8eF7mfnvmflasboQ2LOz\n+9ah7pxvi9j2zSxVZ855XdXqQGB1Z/etQ9053xa96TPu7Gd0PnALsKoL+/Zoe405vTrmNFq8AWNO\nPcecRos3YMzp6zGn0eIN9GDMKTPB2tqH700FftnFfetBd84XIIF7I+LRiPjsNmjfttCpc46IyRHx\nHHAXcMHW7FtnunO+0Ps+4y2eb0SMoBJQ/qEoarmJsxafrzGnb8ecRos3YMyp55jTaPEGjDl9PeY0\nWryBHow5Zc4i2OnZMiLio8C5wJFbu28d6c75AhyZmS9HxB8BcyPi+cx8sOxGlqxT55yZtwG3RcRE\nYGZEHLBtm7XNdOl8gfcXm3rbZ9yZ8/0ucHFmZkQE7/56VYvvsDGnHX0k5jRavAFjTlvqJeY0WrwB\nY07blfpOzGm0eAM9GHPKTLC2+PA9gOIGyOupPAn9la3Zt85053zJzJeL/66KiDlUuh7r/R/mVn1O\nmflgRDQBQ4p6ffIzbtFyvhGxa2b+Vy/8jDtzvocC/1KJOQwFPh4R6zu5b9mMOX075jRavAFjTj3H\nnEaLN2DM6esxp9HiDfRkzGnv5qytfVFJ1n5N5eav/rR9s9xeVG4Q++Ot3bfeXt08352AnYvlAcDD\nwLG1PqeSzvl9vDv9/3jg1338M27vfHvdZ7y1nxHwE+DkWn2+xpy+HXMaLd6UcM598jNuVb9mMafR\n4k0J59wn/z32pZjTaPGmK59Td2JOaT1Y2c7D9yLivGL7tcA3gMHAPxSZ4frMnNDevmW1bVvozvkC\nuwO3FmVNwI2ZeU8NTmOrdPKc/xz4dJHtvw6c1tG+tTiPzurO+dILP+NOnu9W7VsH7TXm9NKY02jx\nBow51HHMabR4A8Yc+njMabR4Az0bc3zQsCRJkiSVpNQHDUuSJElSIzPBkiRJkqSSmGBJkiRJUklM\nsCRJkiSpJCZYkiRJklQSEyxJkiRJKokJliRJkiSVxARLkiRJkkpigqVtJiIm17oNkhqHMUdSTzLm\nqD0mWNomIuIQ4AMR8alat0VS32fMkdSTjDnqSFOtG6DeLSLGAycCS4DlwPsz8yoqyft/A2/VsHmS\n+hhjjqSeZMxRV0Rm1roN6sUi4kjgT4BFmfnLiJiXmR+rdbsk9U3GHEk9yZijrnCIoLolMx8GPgg8\nEBEB7F7jJknqw4w5knqSMUddYYKlMuyama8DRwG/qHVjJPV5xhxJPcmYo63iPVjqloh4H9AUEScA\nhwOX1rhJkvowY46knmTMUVd4D5a6JSLOAjIz/7nWbZHU9xlzJPUkY466wiGC6rKIGA58BhhR67ZI\n6vuMOZJ6kjFHXWUPliRJkiSVxB4sSZIkSSqJCZYkSZIklcQES5IkSZJKYoIlSZIkSSUxwZIkSZKk\nkphgSZIkSVJJTLAkSZIkqSQmWJIkSZJUkv8HEllubyWrwrEAAAAASUVORK5CYII=\n",
       "text": [
        "<matplotlib.figure.Figure at 0x109b8cb10>"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAA1gAAAEYCAYAAABBWFftAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3X2cVdV96P/PFwc0ghGQBFERNddEkShBJS2GtjHqL40B\n0SpRI/EKySUPjaZRU03aSGy1xqg3SdUWE5NSGqhoxKckP0Qq+HALlvqIaEXNXIEKQtEUSVTQ7/3j\nbMbDMDMMM5t5Op/363Ves/faa++z9pzhy/nutfbakZlIkiRJktqvV2c3QJIkSZJ6ChMsSZIkSSqJ\nCZYkSZIklcQES5IkSZJKYoIlSZIkSSUxwZIkSZKkkphgSZIkSVJJTLAkSZIkqSQmWOpyImJZRPxB\nZ7dDUs9lnJG0qxlnaldkZme3oduLiHrgPcDBmfnbouzzwGcz8+Od2baupPg9Tc7Mf+nstuyMiBgJ\nnJOZF1WVTQCGA+8AqzNzZkvlUnsZZ1qnluOM8UftZZxpHeOM8WdH6jq7AT1IL+AC4G86uyGtFRF1\nmbmlA98ygdhV7YmIjwKXAMcCB2XmlogYDPwA6AdcmZn/ZyeP+XXgY8Bvqsr2Bv4yM48u1v81In4J\nbGmi/FeZub6t5yQ1YpzZsVqMM8Yflck4s2PGmR2X13T8cYhgORK4Brio+GPdRkQcHhELI+LVort4\nXKPt9RFxYUQ8ERGvRcQ/R8Tuzb1ZUf+SiHg6IjZExE+21i/Kn4+I/y62T2hi329ExJPAxojYraV9\nivoXRcSTEbExIm6OiMER8auI+E1EzI+I/lX194uIn0fEKxHxYkR8tSifCRwI3F0c56IW2lMfEZ8o\ntg+NiNuL462PiL9t9kPIXAL8/8BzwJ8UZWuBe4AzdjYYFftfB9zZqPgPgOVV608AxzdT7hU/lcU4\n825940zryo0/2lnGmXfrG2faV17T8ccerPIsBRYCFwF/ubUwIuqAu4EfAycAY4E7I+KYzHyuqJbA\nGcD/B7wJPAz8T2B6C+93NnAS8Nvi+H9RvO/zwMcyc01ETAT+KSL+R2auqdr3TOCPgfWZ+XZENLXP\nB4p/zAmcBnwC6A08BnwEOA94FvglcD5weUT0KtoyF/gMMBS4LyL+IzMnRcTHgClNdKk3bk8CGRG7\nUQkm9wGfpdLtfExzv5Di/TcDPwQuBm4pNvXNzN9V1TsE+EILv9vFmVkdhBpfpToAeK1q/TXgUGBD\nM+VSWYwzxpnm4ozxR2UxzhhndjbOGH8aMcEqTwLfBh6OiB9Ulf8elX8QVxXr90fEPcBZwHeq6v1w\na9CIiLuBkTt4r+szc3VR/wrgb6l0z97WUClzTkRcCowG7qra94db9y3qNbfP3UXx32bmuuK9HgTW\nZuYTxfpcKsEKKl3ZgzLzr4v1X0fEj6kEnHtbOJdt2lNlNDAEuDgz3ynKHm7h9zKKyn8My4DrImJU\nZj5avMe7b5j5InBpC8dpqo3V+gNvVK2/RaXLPpspl8pinDHONBdnjD8qi3HGOLOzccb404hDBEuU\nmU9TuUJxCe/+Ee8HrGxU9f8C+zcqq74i81ugX0ScXXQ/b4yIXzSqX33Ml4r3ISI+FxGPRaX7/lVg\nBLBPC/s2t8+gqiprq5Z/12j9Dd79RzQM2G/rcYpjXQq8n5Y1/v1sNRT4v1XBaEeOzMwni/o3Al+N\niA8B/9HK/ZvT+IrPxkZl76Fy9aa5cqk0xhnjzE6WSzvNOGOcKam8ZtmDVb7LgEeBa4v1/wSGRkRk\nNkzZOIxKd3SLMnMWMKuZzQc2Wl4dEQcCP6Iy7vVfMzMj4jG2/wfVcAUjIoYBN1EZQ9vSPtWa27YS\n+HVmfrC5U9rJ8peAAyNit8x8u4X2bFV9weDHVIYXLKdyU2iDNnSpN27fC2zbtb8Plc/8tUblg4py\nqWzGGePMjsqNP2ov44xxpq3lNR9/7MEqWWa+QGWs7AVU/pCXULmC842I6B0RfwR8GvjnFg7TUjDY\nuv3LEbF/RAwEvlW8Z18q43rXA70i4jwqV29a0rdo587s05xHqNzY+Y2IeE9UbvAcERFb/9GtBT6w\nk8d7GbgqIvaMiD0iYkxTFSOiN5UuaQAy8zXgNuDjmflWdd3MfDEzL23h1fgm0MafxwPA0VXrRwML\nmigfVZRLpTLOGGdaUW78UbsYZ4wz7Siv+fhjgrVrXA7sCZCZm4FxVG56XAdcD0zKd28IbUrS/FWQ\nrdtnURkH/AKwAvjrzHyGypWmf6XSRT8CeKilhmbm8p3dp1HbGtpaXJX5NJXx1i9SOd+bgPcWdf8G\n+Iuiu/3rO3gPiq7xccD/oHL1ZyUwsXG9iDiWSkA+KSKqhyr8EHhwR+/Tkoj4U2Ay8EcRcVlEvDcz\nNwFXR8RfRMS3gasz85Xmytvz/lILjDPGGeOPdjXjjHHG+NMGLT5oOCJ+ApwMvJKZHy7KBlL58IcB\n9cDEIrsmKjcTTgbeBs7PzOZuBFQ7RMSvaXr2GqnbioihwD9SGeOewE2Z+cOImAZ8nsp/cADfzMxf\nFfsYc3YR44x6uuKelurel0OozF73T/g9p0MYZ9RT7agH66fAJxuVXQLML8alLijWiYjhVKayHF7s\nc2NUppmUpNbYDPxZZh5BZbaqr0TE4VSSresy8yPFa2tyZcyR1GaZ+R9b4wqV4U2/pTItt99zJLVL\ni4EhMx8EXm1UPB6YUSzPALY+xO0UYHZmbs7Meio35I0ur6mSerLMXJOZjxfLrwPP8O7sVE2N4zfm\nSCrLCcDzmbkSv+dIaqe2XHkZnJUHtkHlJr/BxfJ+wKqqeqvYfupOlSAzD7Y7XT1ZRBxE5QGQi4ui\nr0bEExFxc0T0L8qMObuQcUY15kxgdrHs95wOYpxRT9Wuru1ims4d3by4jX/4h3/Yuo8vX76672uX\niYh+VGZMuqDoyfo74GAqNxu/zLtTBjdlu7YZc3ztzGvatGmd3gZf2712qYjoQ2UCglsbb/N7ji9f\nNflqt7YkWGsjYl+AiBgCbJ0lZDWVB6ltdUBRto36+vo2vKWkWlBMT/tz4J8y8w6AYoaiLL7o/Jh3\nh+QYc1S673znO53dBHW8Pwb+PTO3TqTj9xxJ7dKWBw3fBZwLfLf4eUdV+ayIuI5Kl/mhVOb9l6Qd\niogAbgaWZ+b3q8qHZObLxeqpwFPFsjFHbTb1we2feTp97Nmd0BJ1AWfx7vBA8HuOpHZqMcGKiNnA\nHwKDImIl8G3gKmBOREyhmL4UIDOXR8QcKk+a3gJ8OVuaA16StnUccA7wZEQ8VpR9EzgrIkZS6bb/\nNTAVjDmS2i8i+lKZ4OILVcV+z5HULi0mWJl5VjObTmim/pXAle1tlKTak5kP0fSw5V+1sI8xR1Kb\nFQ9IHdSobAN+z5HUDj6/QZKkdlq3bt2OK0mSaoIJliRJ7XDrrbdy0UUXcfvtt5dyvKVLl3LBBRcw\nc+ZMvvjFL/LCCy80W3fhwoXMnj2bm2++mXPOOYcFCxZss33jxo2cfvrprFy5cpvyyZMn07t3b97z\nnvcwevRoHn300YZtBxxwAAMGDNjm9ZWvfKWUc5OkWtCWSS4kSRLw4osvstdeezFjxgx+8YtfUF9f\nz0EHHdTm47355pucfvrpLFmyhMGDB3P44Ydz1lln8cgjTc+lcMYZZ3DttdcyZcoU+vfvz/jx43nl\nlVfo27cvN998M6tWreL222/nuuuu22a/YcOGsXr1ajKTwYMHN5SvXbuWb3zjG4wfP57KvDNwzTXX\n8Fd/9VdtPidJqjUmWJIktdEhhxzCIYccAsDJJ5/c7uM98MAD9OvXryHpOfroo3nmmWeaTdwWLVrU\nUP7OO++wZcuWhm1TpkwBmp96/v3vf3+T5ZMmTWLAgAEAzJkzh7PPPpv+/fs3WVeStD0TLEmS2mDT\npk3ceOONLF68mMmTJ7NhwwYeffRRxo0bx/HHHw/Aq6++yve+9z1ammyurq6Oyy67jLq6Ourr69ln\nn30atkUEAwYM4Omnn24ywRo+fHjD8ty5c5k2bRp9+/bdYdt/97vfcdNNN/He976XBQsW8PWvf53D\nDz98m96s1atXs3z5ciZOnNiaX4ckqWCCJUlSG8ydO5cvfelL3HnnnbzxxhtMmjSJ0aNHM3HiRJ54\n4gkABgwYwJVXtn7SufXr17PnnntuU7bHHnuwcePGZvd5/PHHWbBgAf369eNrX/taq97nwx/+MBMn\nTqR37968//3vZ8KECTz77LMNwwIBvvWtbzFt2rRWt12SVOEkF5IktcG4cePo3bs3K1as4NRTTwXg\npZdeYv369W0+Zv/+/bfr7Xr99dcZNGhQM3vAyJEjufDCCzn22GMZO3YsmzZt2uH7nHnmmfTu3Ruo\nDHNcsWIFTz31VMP2V155hfvvv79d95NJUq2yB0uSpDbYe++9eeCBBxg9ejS9elWuV86bN4+TTjqp\noc6GDRu45pprWhwiuNtuuzFt2jTq6uo47LDDmD59esO2LVu2sGHDBoYNG7bdfosXL2bChAksWbKE\nYcOGMXbsWKZOncq8efM47bTTmn2/xYsXc+KJJ7J+/Xp23333ht6xPn36NNT51a9+tc1QRUlS65lg\nSZLURgsXLuTII48EKs/Cuvvuu5k/f37D9oEDB+7UEMGxY8eybt06Vq5cydChQ1m0aBFHHHEEhx56\nKAALFixg0KBBHHXUUdTV1TFixAiGDBkCVGY07NOnDyNHjtzuuNUJ3tChQ7n44ovZfffdAXj44Yc5\n7rjjOOywwxrqLFu2bLuhipKk1jHBkiSpjRYtWsSYMWOYPXs2S5Ys4bbbbuPAAw9s8/Hq6uqYOXMm\nV1xxBWPGjOH+++/nlltuadh+ww03MGrUKI466iiOOeYYzjvvPK6//np69erFQw89xD333NMwq+Gs\nWbN46KGHiAguueQSPvaxj/GVr3yF/fffn1GjRnHttdfy9ttvs2LFCubOnbtNO/bee28+9KEPtfk8\nJKmWRUvDFnaFadOmpTfNSt1e7LhK12DMUVOmPjhru7LpY88mIloczlftrbfeYujQoaxZs2abySFU\num71yzXmSN1eu2OOk1xIktQGS5YsYcSIESZXkqRtmGBJkrSTli1bxuWXX8769eu3uedKkiTvwZIk\nqUrj4YPTx569XZ0RI0aYWEmSmmQPliRJkiSVxARLkiRJkkpigiVJkiRJJTHBkjqAU/ZKkiTVBhMs\nqQN85zvf6ewmSJIkqQOYYEmSJElSSZymXeokjaeChqang5YkSVL3YQ+WJEmSJJXEBEvqwtatW9fZ\nTZCkHisi+kfEbRHxTEQsj4iPRsTAiJgfEc9FxL0R0b+q/qURsSIino2Ikzqz7ZK6LhMsqYu69dZb\nueiii7j99ttLOd7SpUu54IILmDlzJl/84hd54YUXmq27cOFCZs+ezc0338w555zDggULGrZNnjyZ\niNgcEb+LiEciYtTWbRGxKiJebfS6oWr7XsWXmaGlnJQktc8PgF9m5uHAkcCzwCXA/Mz8ILCgWCci\nhgOfAYYDnwRujAi/R0najoFB6oJefPFF9tprL2bMmMHuu+9OfX19u4735ptvcvrpp/PNb36TSZMm\n8fnPf56zzjqr2fpnnHEGmzdvZsqUKZx66qmMHz+eTZs2ATBs2DCA/YGDMnN0Zj4KEBGDgauBjwAj\ni9c/Ad8qtk8Bvg6cBkS7TkiS2iki9gbGZuZPADJzS2b+BhgPzCiqzQAmFMunALMzc3Nm1gPPA6M7\nttWSugMnuZC6oEMOOYRDDjkEgJNPPrndx3vggQfo168fgwcPBuDoo4/mmWeeob6+noMOOmi7+osW\nLWoof+edd9iyZcs22zPzlWbeamZmvgoQEROBWZn5WrHPzUX5Ze0+IUlqv4OBdRHxU+Ao4N+BrwGD\nM3NtUWctMLhY3g9YXLX/KioXmyRpGyZYUhezadMmbrzxRhYvXszkyZPZsGEDjz76KOPGjeP4448H\n4NVXX+V73/semdnscerq6rjsssuoq6ujvr6effbZp2FbRDBgwACefvrpJhOs4cOHNyzPnTuXadOm\n0bdvXwB+97vfERH/C/hv4BPAdZn5TNUXEiJif2B4Zs5p1y9DknadOmAU8KeZ+W8R8X2K4YBbZWZG\nRPOBFlraphoybdo0pk2b1tnNUBdhgiV1MXPnzuVLX/oSd955J2+88QaTJk1i9OjRTJw4kSeeeAKA\nAQMGcOWVV7b6mOvXr2fPPffcpmyPPfZg48aNze7z+OOPs2DBAvr168fXvva1hvIPf/jDAD/NzM0R\n8QpwR0Qclttme1cA01rdQEnqeKuAVZn5b8X6bcClwJqI2Dcz10TEEGBrj/1qoPr+0QOKMonvfOc7\nJlhq4D1YUhczbtw4evfuzYoVKzj11FMBeOmll1i/fn2bj9m/f//tertef/11Bg0a1Ow+I0eO5MIL\nL+TYY49l7NixDfdgnXnmmWTm5qLai8ChwIe37hcR7wc+XtyjIEldUmauAVZGxAeLohOAp4G7gXOL\nsnOBO4rlu4AzI6JPRBxMJfY90oFNltRN2IMldTF77703DzzwAKNHj6ZXr8o1kHnz5nHSSe/OCLxh\nwwauueaaFocI7rbbbkybNo26ujoOO+wwpk+f3rBty5YtbNiwYeuEFdtYvHgxEyZMYMmSJQwbNoyx\nY8cydepU5s2bx3777ceJJ57I66+/vntmvgnsVez2VtUh/hj4r3b8CiSpo3wV+FlE9AFeAM4DdgPm\nFBPz1AMTATJzeUTMAZYDW4AvZ0tBWFLNMsGSuqCFCxdy5JFHApVnYd19993Mnz+/YfvAgQN3aojg\n2LFjWbduHStXrmTo0KEsWrSII444gkMPPRSABQsWMGjQII466ijq6uoYMWIEQ4YMASozGvbp04eR\nI0ey++67c/HFF/Ptb3/7zeLQxwEPZ+azVW83AvjtDprkLIKSOl1mPgEc28SmE5qpfyXQ+uArqSaZ\nYEld0KJFixgzZgyzZ89myZIl3HbbbRx44IFtPl5dXR0zZ87kiiuuYMyYMdx///3ccsstDdtvuOEG\nRo0axVFHHcUxxxzDeeedx/XXX0+vXr146KGHuOeeexpmNRw1ahQRcSGVq7yHAqc2ervfAP/RuA0R\ncTbwMSo3hV8VEQ9l5g2N60mSVCvWrVvH+973vs5uhkpmgiV1MW+99RbLli3jvvvuIyJafF7Vzjj+\n+OMbZiH83Oc+t822xg8z/uxnP9uwXD3BBcCnP/1pMvPa5t4nM/+6mfJZwCzgyzvVcEmSeqBbb72V\ne+65h1NOOYXTTjut3cdbunQpM2fO5JhjjuHhhx/m4osv5gMf+ECTdRcuXMjLL7/Mb3/7W+6//37O\nO+88PvGJT+zUcS688EJOO+00jjvuuB0es9Y4yYXUxSxZsoQRI0YQ4Sg6SZJ6ohdffJG99tqLGTNm\nsPvuu1NfX9+u47355pucfvrpfPOb32TSpEl8/vOfb/EC7RlnnMHmzZuZMmUKp556KuPHj2fTpk2t\nPs6iRYv42c9+xttvv73DY9Yie7CkLmTZsmVcfvnlrF+/nvnz53PiiSd2dpMkSVLJDjnkkIah9yef\nfHK7j/fAAw/Qr18/Bg+uPBf76KOP5plnnqG+vr7J510uWrSoofydd95hy5YtrT7Ob37zG5588kkO\nP/zwVh2zFrU5wYqIS4FzgHeAp6jMvNMXuAUYRjHzTma+1v5mSrVhxIgR20xmUUsiYijwj8D7qdyn\ndVNm/jAiBtJMXCni0GTgbeD8zLy3M9ouSVJrbdq0iRtvvJHFixczefJkNmzYwKOPPsq4ceMahvK/\n+uqrfO9732txtuC6ujouu+wy6urqqK+vZ5999mnYFhEMGDCAp59+uskEa/jw4Q3Lc+fOZdq0afTt\n27dVx/nRj37E+eefv93tBc0dsxa1KcGKiIOALwCHZ+abEXELcCZwBDA/M6+OiD+n8kT0S5o9kCS9\nazPwZ5n5eET0A/49IuZTuXizXVyJiOHAZ4DhwP7AfRHxwcx8p7NOQJKkHZk7dy5f+tKXuPPOO3nj\njTeYNGkSo0ePZuLEiTzxxBMADBgwYKdmC16/fj177rnnNmV77LEHGzdubHafxx9/nAULFtCvX7+G\n+613dJy7776bT33qU/Tp06fVx6xFbe3B+m8qX4b2jIi3gT2B/6TyBPQ/LOrMABZigiWpFYqHfq4p\nll+PiGeoJE7jaTqunALMLh56XB8RzwOjgcUd3HRJUg2Z+uCsdu0/btw4evfuzYoVKzj11MpEvC+9\n9BLr169v8zH79++/XW/X66+/zqBBg5rdZ+TIkYwcOZIf/ehHjB07lkWLFrH33ns3e5z//M//5LXX\nXtump6px3aaOWYu9WG1KsDJzQ0RcC7wE/A6Yl5nzI2JwZq4tqq0FBpfUTqlLay7YTh97dge3pGco\nesk/AiwBmosr+7FtMrWKSkImSVKXtffee/PAAw8wevRoevWqzDc3b948TjrppIY6GzZs4Jprrmlx\niOBuu+3GtGnTqKur47DDDmP69OkN27Zs2cKGDRsYNmzYdvstXryYCRMmsGTJEoYNG8bYsWOZOnUq\n8+bN4/DDD+emm25q8jj33nsva9as4bvf/S4Azz33HLNnz2bTpk0MHDiw2WOWMUNid9PWIYIfAL4G\nHETlmTe3RsQ51XUyMyPCJ5xL2inF8MCfAxdk5sbq2RRbEVeMOZKkLm/hwoUceeSRQOVZWHffffc2\n92APHDhwp4YIjh07lnXr1rFy5UqGDh3KokWLOOKIIzj00EMBWLBgAYMGDeKoo46irq6OESNGMGTI\nEKAyo2GfPn0YOXIkBx54YLPH2Xqsrf7+7/+es88+mz/4gz9g6dKlzR6zFrV1iOAxwP/JzP8CiIjb\ngd8H1kTEvpm5JiKGAK+U1E6pQzXVI2Vv1K4XEb2pJFczM/OOonhtM3FlNTC0avcDijJJkrq0RYsW\nMWbMGGbPns2SJUu47bbbOPDAA9t8vLq6OmbOnMkVV1zBmDFjuP/++7nlllsatt9www2MGjWKo446\nimOOOYbzzjuP66+/nl69evHQQw9xzz33NMxq2NJxAFavXs0PfvAD1qxZw7XXXsvGjRs5+eSTWzxm\nrWlrgvUs8JcR8R7gDeAE4BFgE3Au8N3i5x3NHkGSqkSlq+pmYHlmfr9q0100HVfuAmZFxHVUhgYe\nSiUOSZLUZb311lssW7aM++67j4ho8XlVO+P4449vmIXwc5/73DbbGs/499nPfrZhufFkFC0dB2D/\n/ffn6quv5uqrr271MWtNW+/BeiIi/hFYSmWa9keBm4C9gDkRMYViOuWS2imp5zuOyqMfnoyIx4qy\nS4GraCKuZObyiJgDLAe2AF/OlgarS5LUBSxZsoQRI0ZQPQRePUubn4OVmVcDVzcq3kClN0uSdkpm\nPgT0amZzk3ElM68EWj9IXZKkTrRs2TIuv/xy1q9fz/z58znxxBM7u0naBdqcYEmSJElqvREjRmwz\nmYV6puauFkuSJEmSdpI9WNJOcHZBSZIktcQeLEmSJEkqiQmWJEmSJJXEBEuSJEmSSmKCJUmSJEkl\nMcGSJEmSpJKYYEmSJElSSZymXdrFtk7tXj3Fu1O7S1Lni4h64L+Bt4HNmTk6IgYCtwDDgHpgYma+\nVtS/FJhc1D8/M+/tjHZL6trswZIkSbUqgT/KzI9k5uii7BJgfmZ+EFhQrBMRw4HPAMOBTwI3RoTf\noyRtx8AgSZJqWTRaHw/MKJZnABOK5VOA2Zm5OTPrgeeB0UhSIyZYkiSpViVwX0QsjYgvFGWDM3Nt\nsbwWGFws7wesqtp3FbB/xzRTUnfiPViSJKlWHZeZL0fE+4D5EfFs9cbMzIjIFvZvaZukGmUPliRJ\nqkmZ+XLxcx0wl8qQv7URsS9ARAwBXimqrwaGVu1+QFEmSdswwZIkSTUnIvaMiL2K5b7AScBTwF3A\nuUW1c4E7iuW7gDMjok9EHAwcCjzSsa2W1B04RFCSJNWiwcDciIDK96GfZea9EbEUmBMRUyimaQfI\nzOURMQdYDmwBvpyZDhGUtB0TLEmSVHMy89fAyCbKNwAnNLPPlcCVu7hpkro5hwhKkiRJUklMsCRJ\nkiSpJCZYkiRJklQSEyxJkiRJKokJliRJkiSVxARLkiRJkkpigiVJkiRJJTHBkiRJkqSS+KBhqZuY\n+uCs7cqmjz27E1oiSZKk5tiDJUmSJEklMcGSJEmSpJKYYEmSJElSSUywJEmSJKkkJliSJEmSVBIT\nLEmSJEkqiQmWJEmSJJWkzQlWRPSPiNsi4pmIWB4RH42IgRExPyKei4h7I6J/mY2V1LNFxE8iYm1E\nPFVVNi0iVkXEY8Xrj6u2XRoRKyLi2Yg4qXNaLUmS9K72PGj4B8AvM/P0iKgD+gLfAuZn5tUR8efA\nJcVL6lTd6SG93amtu8BPgb8F/rGqLIHrMvO66ooRMRz4DDAc2B+4LyI+mJnvdFRjJUmSGmtTD1ZE\n7A2MzcyfAGTmlsz8DTAemFFUmwFMKKWVkmpCZj4IvNrEpmii7BRgdmZuzsx64Hlg9C5sniRJ0g61\ndYjgwcC6iPhpRDwaET+KiL7A4MxcW9RZCwwupZWSat1XI+KJiLi5aujxfsCqqjqrqPRkSZIkdZq2\nJlh1wCjgxswcBWyi0VDAzEwqQ3skqT3+jspFnZHAy8C1LdQ15kiSpE7V1gRrFbAqM/+tWL+NSsK1\nJiL2BYiIIcAr7W+ipFqWma9kAfgx7w4DXA0Mrap6QFEmSZLUado0yUVmromIlcUN5c8BJwBPF69z\nge8WP+8oraWSalJEDMnMl4vVU4GtMwzeBcyKiOuoDA08FHikE5qoGlDjk8/0aBGxG7CUyoXjcREx\nELgFGAbUAxMz87Wi7qXAZOBt4PzMvLdzWi2pK2vPLIJfBX4WEX2AF4DzgN2AORExhSIotbuFkmpG\nRMwG/hAYFBErgcuAP4qIkVSG//0amAqQmcsjYg6wHNgCfLno5ZKknXEBlTiyV7F+CU3MiOzMpZJa\nq80JVmY+ARzbxKYT2t4cSbUsM89qovgnLdS/Erhy17VIUk8WEQcAnwKuAL5eFI+ncqEHKjMiL6SS\nZDXMXAqzqnlFAAAQpUlEQVTUR8TWmUsXd2SbJXV9bX7QsCRJUjf3v4GLgepeqOZmRHbmUkmtYoIl\nSZJqTkR8GnglMx+j6WfttWZGZIclS9pOe+7BkiSpS3NyCrVgDDA+Ij4F7AG8NyJmAmsjYt9iQq/q\nGZGduVRSq9iDJUmSak5mfjMzh2bmwcCZwL9k5iQqM5SeW1SrnhH5LuDMiOgTEQfjzKWSmmEPliRJ\n0rvD/a6iiRmRnblUUmuZYEmSpJqWmYuARcXyBpqZEdmZSyW1hkMEJUmSJKkkJliSJEmSVBITLEmS\nJEkqifdgSZIkSe3kYyG0lT1YkiRJklQSe7CkHsoraZIkSR3PHixJkiRJKokJliRJkiSVxARLkiRJ\nkkpigiVJkiRJJXGSC6mbczILSZKkrsMeLEmSJEkqiQmWJEmSJJXEBEuSJEmSSuI9WKpp3r8kSZKk\nMtmDJUmSJEklMcGSJEmSpJKYYEmSJElSSUywJEmSJKkkJliSJEmSVBITLEmSJEkqiQmWJEmSJJXE\nBEuSJNWciNgjIpZExOMRsTwi/qYoHxgR8yPiuYi4NyL6V+1zaUSsiIhnI+Kkzmu9pK7MBEuSJNWc\nzHwD+HhmjgSOBD4eER8DLgHmZ+YHgQXFOhExHPgMMBz4JHBjRPg9StJ2DAySJKkmZeZvi8U+wG7A\nq8B4YEZRPgOYUCyfAszOzM2ZWQ88D4zuuNZK6i5MsCRJUk2KiF4R8TiwFrg/M58GBmfm2qLKWmBw\nsbwfsKpq91XA/h3WWEndRl1nN0CSJKkzZOY7wMiI2BuYFxEfb7Q9IyJbOsQubaCkbskeLEmSVNMy\n8zfAL4CjgbURsS9ARAwBXimqrQaGVu12QFEmSdtoV4IVEbtFxGMRcXex3uzMO5K0IxHxk4hYGxFP\nVZU5o5ek0kXEoK3xJCLeA5wIPAbcBZxbVDsXuKNYvgs4MyL6RMTBwKHAIx3bakndQXt7sC4AlvNu\nF3mTM+9IUiv9lMrsXNWc0UvSrjAE+JfiHqwlwN2ZuQC4CjgxIp4Dji/WyczlwBwq33t+BXw5Mx0i\nKGk7bb4HKyIOAD4FXAF8vSgeD/xhsTwDWIhJlqRWyswHI+KgRsXNxZWGGb2A+ojYOqPX4g5prKRu\nLTOfAkY1Ub4BOKGZfa4ErtzFTZPUzbXnau//Bi4G3qkqa27mHUlqK2f0kiRJ3UabEqyI+DTwSmY+\nBkRTdYpuc7vOJZWmFXHFmCNJkjpVW4cIjgHGR8SngD2A90bETIqZdzJzTaOZdySprZqLK87oJUmS\nupw29WBl5jczc2hmHgycCfxLZk6i+Zl3JKmtnNFLkiR1G2U9aHjrsJyrgDkRMQWoByaWdHxJNSAi\nZlOZ0GJQRKwEvk0zcSUzl0fE1hm9tuCMXpIkqQtod4KVmYuARcVyszPvSNKOZOZZzWxyRi9JktQt\n+MwYSZIkSSpJWUMEpQ419cFZ25VNH3t2J7REkiRJepc9WJIkSZJUEhMsSZIkSSqJCZYkSZIklcR7\nsCRJaiPvB5UkNWYPliRJkiSVxARLkiRJkkpigiVJkiRJJTHBkiRJkqSSmGBJkiRJUklMsCRJkiSp\nJE7Trh7FKZMlSZLUmezBkiRJkqSSmGBJkiRJUklMsCRJkiSpJCZYkiSp5kTE0Ii4PyKejohlEXF+\nUT4wIuZHxHMRcW9E9K/a59KIWBERz0bESZ3XekldmQmWJEmqRZuBP8vMI4DfA74SEYcDlwDzM/OD\nwIJinYgYDnwGGA58ErgxIvweJWk7BgZJklRzMnNNZj5eLL8OPAPsD4wHZhTVZgATiuVTgNmZuTkz\n64HngdEd2mhJ3YIJliRJqmkRcRDwEWAJMDgz1xab1gKDi+X9gFVVu62ikpBJ0jZMsCRJUs2KiH7A\nz4ELMnNj9bbMTCBb2L2lbZJqlAmWJEmqSRHRm0pyNTMz7yiK10bEvsX2IcArRflqYGjV7gcUZZK0\njbrOboCkjjX1wVnblU0fe3az5VJ34N+vdlZEBHAzsDwzv1+16S7gXOC7xc87qspnRcR1VIYGHgo8\n0nEtltRdmGBJkqRadBxwDvBkRDxWlF0KXAXMiYgpQD0wESAzl0fEHGA5sAX4cjGEUJK2YYIlSZJq\nTmY+RPO3SpzQzD5XAlfuskZJ6hG8B0uSJEmSSmKCJUmSJEklMcGSJEmSpJKYYEmSJElSSUywJEmS\nJKkkJliSJEmSVBITLEmSJEkqiQmWJEmSJJXEBw2rS5v64KztyqaPPbsTWiJJkiTtmD1YkiRJklSS\nNiVYETE0Iu6PiKcjYllEnF+UD4yI+RHxXETcGxH9y22upFoVEfUR8WREPBYRjxRlxhxJktSltLUH\nazPwZ5l5BPB7wFci4nDgEmB+Zn4QWFCsS1IZEvijzPxIZo4uyow5kiSpS2lTgpWZazLz8WL5deAZ\nYH9gPDCjqDYDmFBGIyWpEI3WjTmSJKlLafc9WBFxEPARYAkwODPXFpvWAoPbe3xJKiRwX0QsjYgv\nFGXGHEmS1KW0axbBiOgH/By4IDM3Rrx7cTkzMyKyne2TpK2Oy8yXI+J9wPyIeLZ6ozFHkiR1BW3u\nwYqI3lSSq5mZeUdRvDYi9i22DwFeaX8TJQky8+Xi5zpgLjAaY44kSepi2jqLYAA3A8sz8/tVm+4C\nzi2WzwXuaLyvJO2siNgzIvYqlvsCJwFPYcyRJEldTFuHCB4HnAM8GRGPFWWXAlcBcyJiClAPTGx3\nCyWpcm/V3GIYch3ws8y8NyKWYsyRJEldSJsSrMx8iOZ7v05oe3MkaXuZ+WtgZBPlGzDmSJKkLqTd\nswhKkiRJkipMsCRJkiSpJCZYkiRJklQSEyxJklRzIuInEbE2Ip6qKhsYEfMj4rmIuDci+ldtuzQi\nVkTEsxFxUue0WlJ3YIIlSZJq0U+BTzYquwSYn5kfBBYU60TEcOAzwPBinxsjwu9QkppkcJAkSTUn\nMx8EXm1UPB6YUSzPACYUy6cAszNzc2bWA89Tedi5JG3HBEuSJKlicGauLZbXUnkGH8B+wKqqequA\n/TuyYZK6DxMsSZKkRjIzgWypSke1RVL3YoIlSZJUsTYi9gWIiCHAK0X5amBoVb0DijJJ2k5dZzdA\nApj64KztyqaPPbsTWiJJqmF3AecC3y1+3lFVPisirqMyNPBQ4JFOaaGkLs8ES5Ik1ZyImA38ITAo\nIlYC3wauAuZExBSgHpgIkJnLI2IOsBzYAny5GEIoSdsxwZIkqWT2ynd9mXlWM5tOaKb+lcCVu65F\nknoK78GSJEmSpJKYYEmSJElSSUywJEmSJKkkJliSJEmSVBITLEmSJEkqiQmWJEmSJJXEBEuSJEmS\nSmKCJUmSJEklMcGSJEmSpJKYYEmSJElSSUywJEmSJKkkdZ3dAEmSWmvqg7O2K5s+9uxOaIkkSU2z\nB0uSJEmSSmKCJUmSJEklMcGSJEmSpJKYYEmSJElSSZzkQh3KG9QlSZLUk5lgSZIkSY14UVht5RBB\nSZIkSSqJPViSJHUQr4hLUs9nD5YkSZIklcQeLEmSJGkXaarnGuy97slK78GKiE9GxLMRsSIi/rzs\n40tSNWOOpI5kzJG0I6X2YEXEbsD1wAnAauDfIuKuzHymzPdR1+H9BOpMxhxJHcmYI6k1yh4iOBp4\nPjPrASLin4FTAAOPpF3BmNNDefFGXZQxR9IOlZ1g7Q+srFpfBXy05PfQLtTcl5qd/bLjlyN1EGOO\negRjZrdhzJG0Q5GZ5R0s4k+AT2bmF4r1c4CPZuZXq+r8mEpAktR91WfmP3R2I4w5Uk3oEvEGjDlS\njWh3zCm7B2s1MLRqfSiNgkxmfr7k95RUu4w5kjqSMUfSDpU9i+BS4NCIOCgi+gCfAe4q+T0kaStj\njqSOZMyRtEOl9mBl5paI+FNgHrAbcLMz60jaVYw5kjqSMUdSa5R6D5YkSZIk1bJShwju6OF7EfHZ\niHgiIp6MiIcj4sjW7tsVtfN864vyxyLikY5tedu14pxPKc75sYj494g4vrX7dkXtPN9u9xm39jOK\niGMjYktxw/dO7duR7TXmdO+YU2vxBow5LdTr9JhTa/EGjDlNbO9RMafW4g10YMzJzFJeVLrKnwcO\nAnoDjwOHN6rz+8DexfIngcWt3bervdpzvsX6r4GBnX0eu+Cc+1Ytf5jK80J68mfc5Pl2x8+4tZ9R\nUe9fgHuAP+msz9eY07NjTq3Fm/aec0/9jKvqdWrMqbV4095z7ql/jz0p5tRavNmZz6mMmFNmD1bD\nw/cyczOw9eF7DTLzXzPzN8XqEuCA1u7bBbXnfLeKXd/MUrXmnDdVrfYD1rd23y6oPee7VXf6jFv7\nGX0VuA1Y14Z9O7S9xpxuHXNqLd6AMacrx5xaizdgzOnpMafW4g10YMwpM8Fq6uF7+7dQfwrwyzbu\n2xW053wBErgvIpZGxBd2Qft2hVadc0RMiIhngF8B5+/Mvl1Me84Xut9nvMPzjYj9qQSUvyuKtt7E\n2RmfrzGnZ8ecWos3YMzpyjGn1uINGHN6esyptXgDHRhzypxFsNWzZUTEx4HJwHE7u28X0p7zBTgu\nM1+OiPcB8yPi2cx8sOxGlqxV55yZdwB3RMRYYGZEHLZrm7XLtOl8gQ8Vm7rbZ9ya8/0+cElmZkQE\n71696ox/w8acZvSQmFNr8QaMOU3pKjGn1uINGHOartRzYk6txRvowJhTZoK1w4fvARQ3QP6IypPQ\nX92ZfbuY9pwvmfly8XNdRMyl0vXY1f8wd+pzyswHI6IOGFjU65Gf8VZbzzci9snM/+qGn3Frzvdo\n4J8rMYdBwB9HxOZW7ls2Y07Pjjm1Fm/AmNOVY06txRsw5vT0mFNr8QY6MuY0d3PWzr6oJGsvULn5\nqw9N3yx3IJUbxH5vZ/ftaq92nu+ewF7Fcl/gYeCkzj6nks75A7w7/f8o4IUe/hk3d77d7jPe2c8I\n+ClwWmd9vsacnh1zai3elHDOPfIzblS/02JOrcWbEs65R/499qSYU2vxpi2fU3tiTmk9WNnMw/ci\nYmqxfTrwbWAA8HdFZrg5M0c3t29ZbdsV2nO+wL7A7UVZHfCzzLy3E05jp7TynP8E+FyR7b8OnNnS\nvp1xHq3VnvOlG37GrTzfndq3C7TXmNNNY06txRsw5tCFY06txRsw5tDDY06txRvo2Jjjg4YlSZIk\nqSSlPmhYkiRJkmqZCZYkSZIklcQES5IkSZJKYoIlSZIkSSUxwZIkSZKkkphgSZIkSVJJTLAkSZIk\nqSQmWJIkSZJUEhMs7TIRMaGz2yCpdhhzJHUkY46aY4KlXSIijgI+EhGf7ey2SOr5jDmSOpIxRy2p\n6+wGqHuLiFHAKcBKYA3wocy8lkry/t/Am53YPEk9jDFHUkcy5qgtIjM7uw3qxiLiOOD3geWZ+cuI\nWJCZn+jsdknqmYw5kjqSMUdt4RBBtUtmPgx8FHggIgLYt5ObJKkHM+ZI6kjGHLWFCZbKsE9mvg4c\nD9zZ2Y2R1OMZcyR1JGOOdor3YKldIuIDQF1EjAOOBS7r5CZJ6sGMOZI6kjFHbeE9WGqXiJgEZGb+\nU2e3RVLPZ8yR1JGMOWoLhwiqzSJiCPB5YP/Oboukns+YI6kjGXPUVvZgSZIkSVJJ7MGSJEmSpJKY\nYEmSJElSSUywJEmSJKkkJliSJEmSVBITLEmSJEkqiQmWJEmSJJXEBEuSJEmSSmKCJUmSJEkl+X9y\nATn1LeTKcQAAAABJRU5ErkJggg==\n",
       "text": [
        "<matplotlib.figure.Figure at 0x109b8ca10>"
       ]
      }
     ],
     "prompt_number": 19
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "# Problem 2\n",
      "\n",
      "Yes, the models are nested. We can express the parameters of $H_0$ in terms of the parameters of $H_A$ ($H_0$ is a subset of $H_A$). Set $\\alpha = (1-\\theta)^2)$, $\\beta = 2\\theta(1-\\theta)$. Thus from the definition we have $P_{H_0}(M) = P_{H_A}(M)$ and $P_{H_0}(MN) = P_{H_A}(MN)$. To prove that $P_{H_0}(N) = P_{H_A}(N),$\n",
      "\n",
      "$$\n",
      "\\begin{align}\n",
      "P_{H_A} &= 1 - \\alpha - \\beta\\\\\n",
      "&= 1 - (1-\\theta)^2) - 2\\theta(1-\\theta)\\\\\n",
      "&= 1 - (1 - 2\\theta + \\theta^2) - 2\\theta + 2\\theta^2\\\\\n",
      "&= \\theta^2\\\\\n",
      "&= P_{H_0}\n",
      "\\end{align}\n",
      "$$\n",
      "\n",
      "Thus the proof is complete.\n",
      "\n",
      "Using the code from HW2, the multinomial distribution for $H_0$ is,\n",
      "\n",
      "$$\n",
      "\\begin{align}\n",
      "P(k_1, k_2, k_3, \\theta) = \\frac{K!}{k_1! k_2! k_3!} \\left((1-\\theta)^2\\right)^{k_1} \\left(2\\theta(1-\\theta)\\right)^{k_2} \\left(\\theta^2\\right)^{k_3}\n",
      "\\end{align}\n",
      "$$\n"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "M = 342.\n",
      "MN = 500.\n",
      "N = 187.\n",
      "\n",
      "# Estimate theta from the data. Self-plagarizing from hw2\n",
      "hat_theta = (MN + 2*N)/(2*(M + MN + N))\n",
      "print 'hat_theta_H0', hat_theta_H0\n",
      "\n",
      "# Estimate alpha, beta from the data.\n",
      "total = M + MN + N\n",
      "hat_alpha = M/total\n",
      "hat_beta = MN/total\n",
      "hat_one_minus_alpha_minus_beta = 1 - hat_alpha - hat_beta"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "hat_theta_H0 0.424684159378\n"
       ]
      }
     ],
     "prompt_number": 38
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%load_ext rmagic"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "The rmagic extension is already loaded. To reload it, use:\n",
        "  %reload_ext rmagic\n"
       ]
      }
     ],
     "prompt_number": 39
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%Rpush M MN N hat_theta total hat_alpha hat_beta hat_one_minus_alpha_minus_beta"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 40
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%%R\n",
      "\n",
      "p1 = (1-hat_theta)^2\n",
      "p2 = 2*hat_theta*(1-hat_theta)\n",
      "p3 = hat_theta^2\n",
      "prob0 = c(p1, p2, p3)\n",
      "\n",
      "probA = c(hat_alpha, hat_beta, hat_one_minus_alpha_minus_beta)\n",
      "\n",
      "x = c(M, MN, N)\n",
      "\n",
      "lik_H0 = dmultinom(x, size=total, prob=prob0)\n",
      "lik_HA = dmultinom(x, size=total, prob=probA)\n",
      "\n",
      "print(paste('likelihood of H0', lik_H0))\n",
      "print(paste('likelihood of HA', lik_HA))\n",
      "\n",
      "lrt = -2*log(lik_H0) + 2*log(lik_HA)\n",
      "print(paste('likelihood ratio test:', lrt))\n",
      "\n",
      "# The LRT is a chi-squared distribution with 1 degree of freedom (2 in HA - 1 in H0)\n",
      "df = 1\n",
      "p_chisq = pchisq(lrt, df)\n",
      "print(paste('chi-squared (df=1) p value:', p_chisq))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "display_data",
       "text": [
        "[1] \"likelihood of H0 0.000887596137156158\"\n",
        "[1] \"likelihood of HA 0.000902136784551232\"\n",
        "[1] \"likelihood ratio test: 0.0324986307405197\"\n",
        "[1] \"chi-squared (df=1) p value: 0.143062350122923\"\n"
       ]
      }
     ],
     "prompt_number": 45
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Since the likelihood ratio test produces a $p$-value greater than our threshold of $\\alpha = 0.05$, we cannot reject the null, and thus $H_0$ is a damn fine model for this study."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "# Problem 3\n",
      "\n",
      "For this problem, $N=100,000$ and $\\bar{x} = \\sum_{i=1}^N x_i$"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "! head BNFO285_HW3_P3.csv"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\"\",\"x\"\r",
        "\r\n",
        "\"1\",23\r",
        "\r\n",
        "\"2\",22\r",
        "\r\n",
        "\"3\",30\r",
        "\r\n",
        "\"4\",17\r",
        "\r\n",
        "\"5\",26\r",
        "\r\n",
        "\"6\",16\r",
        "\r\n",
        "\"7\",37\r",
        "\r\n",
        "\"8\",7\r",
        "\r\n",
        "\"9\",26\r",
        "\r\n"
       ]
      }
     ],
     "prompt_number": 48
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import pandas as pd\n",
      "unknown_discrete_distribution = pd.read_csv('BNFO285_HW3_P3.csv', index_col=0, squeeze=True)\n",
      "\n",
      "# can I tell what distribution this is from just looking at it? :)\n",
      "ppl.hist(unknown_discrete_distribution.values, bins=20)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 65,
       "text": [
        "<matplotlib.axes._subplots.AxesSubplot at 0x1093de6d0>"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAEpVJREFUeJzt3X+s3XV9x/HnaxQYKoE0bpUfFZrlklCCgmzUzRHNJKwu\nC2BGoDBJNzvTpE6Zf2yzLtn8q9EsiiQLhMxiCxsdRBRqgozKTFazYJGBVktHSehGr7aYoqBbBi2+\n98f5VI93t73l3lvv554+H8lJP+f9/X4On3e5Oa9+f5xzU1VIktSbX5rrBUiSNBkDSpLUJQNKktQl\nA0qS1CUDSpLUJQNKktSlIwZUksVJvprkO0m+neTDrf7xJHuSPNEe7xmaszbJriQ7k1wxVL8kyfa2\n7Zah+slJ7mn1R5OccywalSTNL1MdQR0APlJVFwBvBz6Y5HyggE9X1cXt8WWAJEuB64ClwHLg1iRp\nr3UbsKqqxoCxJMtbfRWwv9VvBj45i/1JkuapIwZUVe2tqifb+MfAU8BZbXMmmXIVsKmqDlTVbuAZ\nYFmSM4BTq2pb2+9O4Oo2vhLY2Mb3Ae+eZi+SpBFy1NegkpwLXAw82kofSvLNJOuTnN5qZwJ7hqbt\nYRBoE+vj/CzozgKeA6iqg8CLSRa+tjYkSaPmqAIqyRuAzwM3tSOp24AlwEXA94BPHbMVNhs2bCgG\npxZ9+PDhw8f8eUzblAGV5EQGp97+oaruB6iq56sBPgtc2nYfBxYPTT+bwZHTeBtPrB+a8+b231oA\nnFZVL0xcx+7du4++K0nSvDfVXXwB1gM7quozQ/UzhnZ7L7C9jTcDK5KclGQJMAZsq6q9wEtJlrXX\nvBF4YGjOyja+Bnhkhj1JkkbAgim2vwN4H/CtJE+02seA65NcxODw7VlgNUBV7UhyL7ADOAisqZ99\nXfoaYANwCvBgVT3U6uuBu5LsAvYDK2ajMUnS/HbEgKqqrzH5UdaXjzBnHbBukvrjwIWT1F8Grp1y\npZKk44rfJCFJ6pIBJUnqkgElSeqSASVJ6pIBJUnqkgElSeqSASVJ6pIBJUnqkgElSeqSASVJ6pIB\nJUnqkgElSeqSASVJ6pIBJUnqkgElSeqSASVJ6pIBJUnqkgElSeqSASVJ6pIBJUnqkgElSeqSASVJ\n6pIBJUnqkgElSeqSASVJ6pIBJUnqkgElSeqSASVJ6pIBJUnqkgElSerSgrlegF671VvvntH82y+7\nYZZWIknHjkdQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQu\nHTGgkixO8tUk30ny7SQfbvWFSbYkeTrJw0lOH5qzNsmuJDuTXDFUvyTJ9rbtlqH6yUnuafVHk5xz\nLBqVJM0vUx1BHQA+UlUXAG8HPpjkfOCjwJaqOg94pD0nyVLgOmApsBy4NUnaa90GrKqqMWAsyfJW\nXwXsb/WbgU/OWneSpHnriAFVVXur6sk2/jHwFHAWcCWwse22Ebi6ja8CNlXVgaraDTwDLEtyBnBq\nVW1r+905NGf4te4D3j3TpiRJ899RX4NKci5wMfB1YFFV7Wub9gGL2vhMYM/QtD0MAm1ifbzVaX8+\nB1BVB4EXkyx8LU1IkkbPUQVUkjcwOLq5qap+NLytqgqoY7A2SdJxbMqASnIig3C6q6rub+V9Sd7U\ntp8BPN/q48DioelnMzhyGm/jifVDc97cXmsBcFpVvTCtbiRJI2Oqu/gCrAd2VNVnhjZtBla28Urg\n/qH6iiQnJVkCjAHbqmov8FKSZe01bwQemOS1rmFw04Uk6Tg31W/UfQfwPuBbSZ5otbXAJ4B7k6wC\ndgPXAlTVjiT3AjuAg8CadgoQYA2wATgFeLCqHmr19cBdSXYB+4EVs9CXJGmeO2JAVdXXOPxR1uWH\nmbMOWDdJ/XHgwknqL9MCTpKkQ/wmCUlSlwwoSVKXDChJUpcMKElSlwwoSVKXprrNXMfA6q13T3vu\n7ZfdMIsrkaR+eQQlSeqSASVJ6pIBJUnqkgElSeqSASVJ6pIBJUnqkgElSeqSASVJ6pIBJUnqkgEl\nSeqSASVJ6pIBJUnqkgElSeqSASVJ6pIBJUnqkgElSeqSASVJ6pIBJUnqkgElSeqSASVJ6pIBJUnq\nkgElSerSgrlewHyzeuvdM5p/+2U3zNJKJGm0eQQlSeqSASVJ6pIBJUnqkgElSeqSASVJ6pIBJUnq\nkgElSeqSASVJ6pIBJUnqkgElSeqSASVJ6tKUAZXkjiT7kmwfqn08yZ4kT7THe4a2rU2yK8nOJFcM\n1S9Jsr1tu2WofnKSe1r90STnzGaDkqT56WiOoD4HLJ9QK+DTVXVxe3wZIMlS4DpgaZtza5K0ObcB\nq6pqDBhLcug1VwH7W/1m4JMz6kiSNBKmDKiq2gr8YJJNmaR2FbCpqg5U1W7gGWBZkjOAU6tqW9vv\nTuDqNr4S2NjG9wHvPvrlS5JG1UyuQX0oyTeTrE9yequdCewZ2mcPcNYk9fFWp/35HEBVHQReTLJw\nBuuSJI2A6QbUbcAS4CLge8CnZm1FkiQxzYCqquerAT4LXNo2jQOLh3Y9m8GR03gbT6wfmvNmgCQL\ngNOq6oXprEuSNDqmFVDtmtIh7wUO3eG3GViR5KQkS4AxYFtV7QVeSrKs3TRxI/DA0JyVbXwN8Mh0\n1iRJGi1T/sr3JJuAdwJvTPIc8DfAu5JcxOBuvmeB1QBVtSPJvcAO4CCwph1lAawBNgCnAA9W1UOt\nvh64K8kuYD+wYpZ6kyTNY1MGVFVdP0n5jiPsvw5YN0n9ceDCSeovA9dOtQ5J0vHFb5KQJHXJgJIk\ndcmAkiR1yYCSJHXJgJIkdcmAkiR1yYCSJHXJgJIkdWnKD+pqNK3eeveM5t9+2Q2ztBJJmpxHUJKk\nLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4Z\nUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuGVCS\npC4ZUJKkLhlQkqQuGVCSpC4ZUJKkLhlQkqQuLZjrBWj+Wr317mnPvf2yG2ZxJZJGkUdQkqQuGVCS\npC5NGVBJ7kiyL8n2odrCJFuSPJ3k4SSnD21bm2RXkp1JrhiqX5Jke9t2y1D95CT3tPqjSc6ZzQYl\nSfPT0RxBfQ5YPqH2UWBLVZ0HPNKek2QpcB2wtM25NUnanNuAVVU1BowlOfSaq4D9rX4z8MkZ9CNJ\nGhFTBlRVbQV+MKF8JbCxjTcCV7fxVcCmqjpQVbuBZ4BlSc4ATq2qbW2/O4fmDL/WfcC7p9GHJGnE\nTPca1KKq2tfG+4BFbXwmsGdovz3AWZPUx1ud9udzAFV1EHgxycJprkuSNCJmfJNEVRVQs7AWSZJ+\naroBtS/JmwDa6bvnW30cWDy039kMjpzG23hi/dCcN7fXWgCcVlUvTHNdkqQRMd2A2gysbOOVwP1D\n9RVJTkqyBBgDtlXVXuClJMvaTRM3Ag9M8lrXMLjpQpJ0nJvymySSbALeCbwxyXPAXwOfAO5NsgrY\nDVwLUFU7ktwL7AAOAmvaKUCANcAG4BTgwap6qNXXA3cl2QXsB1bMTmuSpPlsyoCqqusPs+nyw+y/\nDlg3Sf1x4MJJ6i/TAk6SpEP8JglJUpcMKElSlwwoSVKXDChJUpcMKElSlwwoSVKXDChJUpcMKElS\nlwwoSVKXDChJUpcMKElSlwwoSVKXDChJUpcMKElSlwwoSVKXDChJUpcMKElSlwwoSVKXDChJUpcM\nKElSlwwoSVKXDChJUpcMKElSlwwoSVKXDChJUpcMKElSlwwoSVKXDChJUpcMKElSlwwoSVKXDChJ\nUpcMKElSlwwoSVKXDChJUpcMKElSlwwoSVKXDChJUpcMKElSlwwoSVKXDChJUpcWzGRykt3AS8Cr\nwIGqujTJQuAe4BxgN3BtVf2w7b8WeH/b/8NV9XCrXwJsAH4ZeLCqbprJug5n9da7ZzT/9stumKWV\nSJKmMtMjqALeVVUXV9WlrfZRYEtVnQc80p6TZClwHbAUWA7cmiRtzm3AqqoaA8aSLJ/huiRJ89xs\nnOLLhOdXAhvbeCNwdRtfBWyqqgNVtRt4BliW5Azg1Kra1va7c2iOJOk4NaNTfAyOoL6S5FXg9qr6\ne2BRVe1r2/cBi9r4TODRobl7gLOAA218yHira8R5ylXSkcw0oN5RVd9L8ivAliQ7hzdWVSWpGf43\nJEnHoRmd4quq77U/vw98EbgU2JfkTQDt9N3zbfdxYPHQ9LMZHDmNt/FwfXwm65IkzX/TDqgkr0ty\nahu/HrgC2A5sBla23VYC97fxZmBFkpOSLAHGgG1VtRd4KcmydtPEjUNzJEnHqZmc4lsEfLHdiLcA\n+MeqejjJN4B7k6yi3WYOUFU7ktwL7AAOAmuq6tDpvzUMbjM/hcFt5g/NYF2SpBEw7YCqqmeBiyap\nvwBcfpg564B1k9QfBy6c7lokSaPHb5KQJHXJgJIkdcmAkiR1yYCSJHXJgJIkdcmAkiR1yYCSJHXJ\ngJIkdcmAkiR1yYCSJHXJgJIkdWmmvw9KmlP+0kNpdHkEJUnqkgElSeqSASVJ6pIBJUnqkgElSeqS\nASVJ6pIBJUnqkgElSeqSASVJ6pIBJUnqkgElSeqSASVJ6pIBJUnqkgElSeqSASVJ6pIBJUnqkgEl\nSeqSASVJ6pIBJUnqkgElSerSgrlegDTXVm+9e9pzb7/shllciaRhHkFJkrpkQEmSumRASZK6NG+u\nQb36k5/wyqsHpz3/pBPmTauSJOZRQP3r3mf42GMPTGvuklPfyAcveOcsr0iSdCzNm4B65dWD/OjA\ny9Oa+z8HX5nl1UiSjrV5E1BSr2Zymzp4q7p0ON4kIUnqUjcBlWR5kp1JdiX5y7lejyRpbnVxii/J\nCcDfAZcD48BjSTZX1VNzuzLpF8PThNL/18sR1KXAM1W1u6oOAP8EXDXHa5IkzaEujqCAs4Dnhp7v\nAZYN7/Abv3oOH3rL5dN68VNOOHH6K5MkzYlU1VyvgSR/ACyvqg+05+8DllXVh4b2+SyD4JIkzR+7\nq2rDdCb2cgQ1Diweer6YCWFUVX/yC12RJGlO9XIN6hvAWJJzk5wEXAdsnuM1SZLmUBdHUFV1MMmf\nAv8MnACs9w4+STq+dXENSpKkiXo5xXdYx8MHeJMsTvLVJN9J8u0kH271hUm2JHk6ycNJTp/rtc62\nJCckeSLJl9rz46Hn05N8PslTSXYkWTbqfSdZ236+tye5O8nJo9ZzkjuS7Euyfah22B7b38mu9v52\nxdysemYO0/Pftp/tbyb5QpLThra9pp67DqihD/AuB5YC1yc5f25XdUwcAD5SVRcAbwc+2Pr8KLCl\nqs4DHmnPR81NwA7g0KH88dDzLcCDVXU+8BZgJyPcd5JzgQ8Ab6uqCxmcxl/B6PX8OQbvVcMm7THJ\nUgbX2pe2Obcm6fr9+DAm6/lh4IKqeivwNLAWptdz738hx8UHeKtqb1U92cY/Bp5i8NmwK4GNbbeN\nwNVzs8JjI8nZwO8BnwXSyqPe82nAZVV1Bwyuv1bVi4x23y8x+EfY65IsAF4HfJcR67mqtgI/mFA+\nXI9XAZuq6kBV7QaeYfB+N69M1nNVbamqn7SnXwfObuPX3HPvATXZB3jPmqO1/EK0f21ezOB/7KKq\n2tc27QMWzdGyjpWbgT8HfjJUG/WelwDfT/K5JP+e5O+TvJ4R7ruqXgA+BfwXg2D6YVVtYYR7HnK4\nHs/k5z9KM6rvbe8HHmzj19xz7wF1XN3BkeQNwH3ATVX1o+FtNbibZWT+PpL8PvB8VT3Bz46efs6o\n9dwsAN4G3FpVbwP+mwmntkat7yS/BvwZcC6DN6k3tA/j/9So9TyZo+hxpPpP8lfAK1V1pC+aPGLP\nvQfUlB/gHRVJTmQQTndV1f2tvC/Jm9r2M4Dn52p9x8BvAVcmeRbYBPxOkrsY7Z5h8PO7p6oea88/\nzyCw9o5w378O/FtV7a+qg8AXgN9ktHs+5HA/zxPf285utZGQ5I8YnL7/w6Hya+6594A6Lj7AmyTA\nemBHVX1maNNmYGUbrwTunzh3vqqqj1XV4qpawuCC+b9U1Y2McM8wuN4IPJfkvFa6HPgO8CVGt++d\nwNuTnNJ+1i9ncGPMKPd8yOF+njcDK5KclGQJMAZsm4P1zbokyxmcur+qqv53aNNr77mqun4A7wH+\ng8EFtbVzvZ5j1ONvM7gO8yTwRHssBxYCX2FwJ8zDwOlzvdZj1P87gc1tPPI9A28FHgO+yeBo4rRR\n7xv4CwZBvJ3BzQInjlrPDM4EfBd4hcG18z8+Uo/Ax9r72k7gd+d6/bPU8/uBXcB/Dr2X3Trdnv2g\nriSpS72f4pMkHacMKElSlwwoSVKXDChJUpcMKElSlwwoSVKXDChJUpf+D8/muonmMzxpAAAAAElF\nTkSuQmCC\n",
       "text": [
        "<matplotlib.figure.Figure at 0x10df5bc10>"
       ]
      }
     ],
     "prompt_number": 65
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "eh it really could be anything. My guess is poisson.\n",
      "\n",
      "## Binomial\n",
      "For the Binomial distribution, we have\n",
      "\n",
      "$$\n",
      "P(x) = \\binom{n}{x} p^k (1-p)^{n-x}\n",
      "$$\n",
      "\n",
      "our maximum likelihood estimation for $\\hat{n}$ is just the maximum of our data, and then $\\hat{p} = \\frac{\\sum_{i=1}^N x_i}{n\\times N}$."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from scipy.stats import binom\n",
      "\n",
      "bar_x = float(unknown_discrete_distribution.mean())\n",
      "N = unknown_discrete_distribution.shape[0]\n",
      "\n",
      "hat_n = unknown_discrete_distribution.max()\n",
      "hat_p = float(unknown_discrete_distribution.sum())/(hat_n * N)\n",
      "\n",
      "# print hat_n, hat_p\n",
      "\n",
      "binom_loglik = binom.logpmf(unknown_discrete_distribution, hat_n, hat_p).sum()\n",
      "print binom_loglik"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "-546551.738242\n"
       ]
      }
     ],
     "prompt_number": 143
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Poisson\n",
      "\n",
      "For the Poisson distribution, our estimate for the parameter $\\lambda$ is just the mean"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from scipy.stats import poisson\n",
      "\n",
      "hat_lambda = float(unknown_discrete_distribution.mean())\n",
      "\n",
      "poisson_loglik = poisson.logpmf(unknown_discrete_distribution, hat_lambda).sum()\n",
      "print poisson_loglik"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "-484817.630774\n"
       ]
      }
     ],
     "prompt_number": 67
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Geometric\n",
      "\n",
      "For the Geometric distribution, our estimate for the parameter is the inverse of the mean. At first, I was getting loglikelihoods of negative infinity, then I guessed that this question was written by Boyko and thus the data created in MATLAB. Then I realized that MATLAB uses $P(x) = p(1-p)^x$ while Python uses $P(x) = p(1-p)^{x+1}$ and thus for this problem, I added 1 to the `unknown_discrete_distribution` to shift the data to match this parameterization of the Geometric distribution."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from scipy.stats import geom\n",
      "\n",
      "hat_p = 1./float((unknown_discrete_distribution+1).mean())\n",
      "\n",
      "print hat_p\n",
      "\n",
      "rv = geom(hat_p)\n",
      "\n",
      "geom_loglik = rv.logpmf(unknown_discrete_distribution+1).sum()\n",
      "print geom_loglik"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "0.0385212315546\n",
        "-423703.291432"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n"
       ]
      }
     ],
     "prompt_number": 107
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Negative Binomial\n",
      "\n",
      "For the negative binomial distribution, we have that the parameters $p$ and $r$ can be estimated by,\n",
      "\n",
      "$$\n",
      "\\hat{p} = \\frac{\\sum_{i=1}^N x_i}{Nr + \\sum_{i=1}^N x_i}\n",
      "$$\n",
      "\n",
      "and\n",
      "\n",
      "$$\n",
      "\\frac{\\partial \\ell}{\\partial r} = \\sum_{i=1}^N\\psi(x_i + r) - N\\psi(r) + N \\ln\\left(\\frac{r}{r + \\sum_{i=1}^N x_i/N}\\right) = 0\n",
      "$$\n",
      "\n",
      "Where $\\psi(x) = \\frac{\\Gamma^\\prime(x)}{\\Gamma(x)}$ aka the \"digamma\" function.\n",
      "\n",
      "Note: to solve for this, I tried initial values of 10, 20, 30, 40, 50, 60, 70, 80, 90, 100 and when\n",
      "I used values above 10 then I got values of r that were greater than the largest\n",
      "value in the dataset (like in the hundreds of thousands), so there's no way it could represent a number of failures, $r$, in this data.\n",
      "Thus I choose the root found with an initial value of 10, giving an estimate of \n",
      "r of around 6.5, which is consistent with the data."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from scipy.stats import nbinom\n",
      "\n",
      "# Solve for roots using fsolve\n",
      "from scipy.optimize import fsolve\n",
      "\n",
      "from scipy.optimize import root\n",
      "\n",
      "# import \"psi\", the digamma function\n",
      "from scipy.special import psi\n",
      "\n",
      "N = unknown_discrete_distribution.shape[0]\n",
      "\n",
      "# Function to estimate r for \n",
      "def estimate_hat_r(r):\n",
      "    return psi(unknown_discrete_distribution + r).sum() - N*psi(r) + N*np.log(float(r)/(r + unknown_discrete_distribution.mean()))\n",
      "\n",
      "# I tried initial values of 10, 20, 30, 40, 50, 60, 70, 80, 90, 100 and when\n",
      "# I used values above 10 then I got values of r that were greater than the largest\n",
      "# value in the dataset, so there's no way it could represent a number of failures.\n",
      "# Thus I choose the root found with an initial value of 10, giving an estimate of \n",
      "# r of around 6.5, which is consistent with the data.\n",
      "hat_r = fsolve(estimate_hat_r, 10)\n",
      "hat_p = float(unknown_discrete_distribution.sum())/(N*hat_r + unknown_discrete_distribution.sum())\n",
      "\n",
      "nbinom_loglik = nbinom.logpmf(unknown_discrete_distribution, hat_r, hat_p).sum()\n",
      "print nbinom_loglik"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "-2849478.97539\n"
       ]
      }
     ],
     "prompt_number": 140
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Uniform\n",
      "\n",
      "For the uniform distribution, we will use the minimum and maxiumum as the parameters"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from scipy.stats import uniform\n",
      "\n",
      "hat_a = unknown_discrete_distribution.min()\n",
      "hat_b = unknown_discrete_distribution.max()\n",
      "\n",
      "uniform_loglik = uniform.logpdf(unknown_discrete_distribution, hat_a, hat_b).sum()\n",
      "print uniform_loglik"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "-469134.788222\n"
       ]
      }
     ],
     "prompt_number": 141
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def AIC(k, loglik):\n",
      "    \"\"\"Akaike information criterion\"\"\"\n",
      "    return 2*k - 2*loglik\n",
      "\n",
      "logliks = pd.Series({'binomial': binom_loglik,\n",
      "                        'poisson': poisson_loglik,\n",
      "                        'geometric': geom_loglik,\n",
      "                        'negative binomial': nbinom_loglik,\n",
      "                        'uniform':uniform_loglik}, name='loglik')\n",
      "n_parameters = pd.Series({'binomial':2,\n",
      "                'poisson':1,\n",
      "                'geometric': 1,\n",
      "                'negative binomial': 2,\n",
      "                'uniform':2}, name='n_parameters')\n",
      "\n",
      "df = pd.concat([n_parameters, logliks], axis=1)\n",
      "print df\n",
      "\n",
      "print df.apply(lambda x: AIC(*x), axis=1)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "                   n_parameters          loglik\n",
        "binomial                      2  -546551.738242\n",
        "geometric                     1  -423703.291432\n",
        "negative binomial             2 -2849478.975392\n",
        "poisson                       1  -484817.630774\n",
        "uniform                       2  -469134.788222\n",
        "\n",
        "[5 rows x 2 columns]\n",
        "binomial             1093107.476483\n",
        "geometric             847408.582865\n",
        "negative binomial    5698961.950783\n",
        "poisson               969637.261547\n",
        "uniform               938273.576445\n",
        "dtype: float64\n"
       ]
      }
     ],
     "prompt_number": 152
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Thus the model with the minimum AIC is the **Geometric** distribution"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "# Problem 4\n",
      "\n",
      "## Problem 4a\n",
      "\n",
      "* Sensitivity = True positives/all positives\n",
      "* Specificity = true negatives/all negatives\n"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "sensitivity = 0.9\n",
      "specificity = 0.95\n",
      "\n",
      "cancer_patients = 32\n",
      "normal_patients = 103\n",
      "total = cancer_patients + normal_patients\n",
      "\n",
      "# Using sensitivity to calculate true positives and false negatives\n",
      "true_positives = int(round(sensitivity * cancer_patients))\n",
      "false_negatives = cancer_patients - true_positives\n",
      "\n",
      "# Using specificity to calcualte false positives and true negatives\n",
      "true_negatives = int(round(specificity * normal_patients))\n",
      "false_positives = normal_patients - true_negatives\n",
      "\n",
      "print 'true_positives:', true_positives, ' false_negatives:', false_negatives\n",
      "print 'true_negatives:', true_negatives, ' false_positives:', false_positives"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "true_positives: 29  false_negatives: 3\n",
        "true_negatives: 98  false_positives: 5\n"
       ]
      }
     ],
     "prompt_number": 2
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Problem 4b\n",
      "\n",
      "Positive Predictive Value (PPV) = TP/(TP + FP)"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "positive_predictive_value = true_positives/float(true_positives + false_positives)\n",
      "print 'positive_predictive_value:', positive_predictive_value"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "positive_predictive_value: 0.852941176471\n"
       ]
      }
     ],
     "prompt_number": 3
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "For the general population, we would want to have maximum coverage, such that we miss the smallest number of true cases. That way, even if a patient is \"over-diagnosed\" then a CAT scan or other follow up studies will show that the patient doesn't have lung cancer. Thus we want a test that has a very low false negative rate, but PPV doesn't measure that. Instead, sensitivity is a better measure."
     ]
    }
   ],
   "metadata": {}
  }
 ]
 }