Bayesian_Analysis_of_Poisson_Model_with_R.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "kernelspec": {
      "display_name": "R",
      "language": "R",
      "name": "ir"
    },
    "language_info": {
      "codemirror_mode": "r",
      "file_extension": ".r",
      "mimetype": "text/x-r-source",
      "name": "R",
      "pygments_lexer": "r",
      "version": "3.3.1"
    },
    "colab": {
      "name": "Bayesian Analysis of Poisson Model with R.ipynb",
      "provenance": [],
      "collapsed_sections": []
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "pT7yQBgzpM9g"
      },
      "source": [
        "### **Introduction**\n",
        "    \n",
        "* Suppose $Y \\sim p(y|\\theta)$ and $p(\\theta)$ is the prior distribution over $\\theta$\n",
        "\n",
        "* Our objective is to make statistical inference about $\\theta$\n",
        "\n",
        "* The **posterior distribution** is:\n",
        "$$\n",
        "p(\\theta|y)=\\frac{p(y|\\theta)p(\\theta)}{\\int_{\\Theta}p(y|\\theta)p(\\theta)d\\theta}\n",
        "$$\n",
        "\n",
        "* **Posterior mean**: *integration problem*\n",
        "$$\n",
        "\\mathbb{E}(\\theta|y)=\\int_{\\Theta}\\theta p(\\theta|y)d\\theta =g(y)\n",
        "$$\n",
        "\n",
        "* **Posterior median**: *integration problem*\n",
        "$$\n",
        "\\int_{-\\infty}^{\\mu_0} p(\\theta|y)d\\theta =\\frac{1}{2},\n",
        "$$\n",
        "\n",
        "* **Posterior Mode**: *optimization problem*\n",
        "\\begin{eqnarray*}\n",
        "\\hat{\\theta}&=&\\text{argmax}_{\\Theta}~p(\\theta|y)\\\\\n",
        "&=&\\text{argmax}_{\\Theta}p(y|\\theta)p(\\theta)\n",
        "\\end{eqnarray*}\n",
        "\n",
        "* **Posterior predictive distribution**:\n",
        "$$\n",
        "p(\\tilde{y}|y)=\\int_{\\Theta}p(\\tilde{y}|\\theta,y)p(\\theta|y)d\\theta\n",
        "$$\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "kqPzsKcpRy2q"
      },
      "source": [
        "Download the data from the following website: https://www.football-data.co.uk/"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "6F5Meds-pM9m",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 325
        },
        "outputId": "8cec0031-f76e-4af5-a955-5c99e9ef94ab"
      },
      "source": [
        "data<-read.csv(file='https://www.football-data.co.uk/mmz4281/1920/E0.csv',header = TRUE)\n",
        "head(data)"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "  Div Date       Time  HomeTeam       AwayTeam         FTHG FTAG FTR HTHG HTAG\n",
              "1 E0  09/08/2019 20:00 Liverpool      Norwich          4    1    H   4    0   \n",
              "2 E0  10/08/2019 12:30 West Ham       Man City         0    5    A   0    1   \n",
              "3 E0  10/08/2019 15:00 Bournemouth    Sheffield United 1    1    D   0    0   \n",
              "4 E0  10/08/2019 15:00 Burnley        Southampton      3    0    H   0    0   \n",
              "5 E0  10/08/2019 15:00 Crystal Palace Everton          0    0    D   0    0   \n",
              "6 E0  10/08/2019 15:00 Watford        Brighton         0    3    A   0    1   \n",
              "  ⋯ AvgC.2.5.1 AHCh  B365CAHH B365CAHA PCAHH PCAHA MaxCAHH MaxCAHA AvgCAHH\n",
              "1 ⋯ 3.43       -2.25 1.91     1.99     1.94  1.98  1.99    2.07    1.90   \n",
              "2 ⋯ 2.91        1.75 1.95     1.95     1.96  1.97  2.07    1.98    1.97   \n",
              "3 ⋯ 1.92       -0.50 1.95     1.95     1.98  1.95  2.00    1.96    1.96   \n",
              "4 ⋯ 1.71        0.00 1.87     2.03     1.89  2.03  1.90    2.07    1.86   \n",
              "5 ⋯ 1.71        0.25 1.82     2.08     1.97  1.96  2.03    2.08    1.96   \n",
              "6 ⋯ 1.73       -0.50 2.04     1.86     2.05  1.88  2.12    1.91    2.05   \n",
              "  AvgCAHA\n",
              "1 1.99   \n",
              "2 1.92   \n",
              "3 1.92   \n",
              "4 2.02   \n",
              "5 1.93   \n",
              "6 1.84   "
            ],
            "text/latex": "A data.frame: 6 × 106\n\\begin{tabular}{r|lllllllllllllllllllll}\n  & Div & Date & Time & HomeTeam & AwayTeam & FTHG & FTAG & FTR & HTHG & HTAG & ⋯ & AvgC.2.5.1 & AHCh & B365CAHH & B365CAHA & PCAHH & PCAHA & MaxCAHH & MaxCAHA & AvgCAHH & AvgCAHA\\\\\n  & <chr> & <chr> & <chr> & <chr> & <chr> & <int> & <int> & <chr> & <int> & <int> & ⋯ & <dbl> & <dbl> & <dbl> & <dbl> & <dbl> & <dbl> & <dbl> & <dbl> & <dbl> & <dbl>\\\\\n\\hline\n\t1 & E0 & 09/08/2019 & 20:00 & Liverpool      & Norwich          & 4 & 1 & H & 4 & 0 & ⋯ & 3.43 & -2.25 & 1.91 & 1.99 & 1.94 & 1.98 & 1.99 & 2.07 & 1.90 & 1.99\\\\\n\t2 & E0 & 10/08/2019 & 12:30 & West Ham       & Man City         & 0 & 5 & A & 0 & 1 & ⋯ & 2.91 &  1.75 & 1.95 & 1.95 & 1.96 & 1.97 & 2.07 & 1.98 & 1.97 & 1.92\\\\\n\t3 & E0 & 10/08/2019 & 15:00 & Bournemouth    & Sheffield United & 1 & 1 & D & 0 & 0 & ⋯ & 1.92 & -0.50 & 1.95 & 1.95 & 1.98 & 1.95 & 2.00 & 1.96 & 1.96 & 1.92\\\\\n\t4 & E0 & 10/08/2019 & 15:00 & Burnley        & Southampton      & 3 & 0 & H & 0 & 0 & ⋯ & 1.71 &  0.00 & 1.87 & 2.03 & 1.89 & 2.03 & 1.90 & 2.07 & 1.86 & 2.02\\\\\n\t5 & E0 & 10/08/2019 & 15:00 & Crystal Palace & Everton          & 0 & 0 & D & 0 & 0 & ⋯ & 1.71 &  0.25 & 1.82 & 2.08 & 1.97 & 1.96 & 2.03 & 2.08 & 1.96 & 1.93\\\\\n\t6 & E0 & 10/08/2019 & 15:00 & Watford        & Brighton         & 0 & 3 & A & 0 & 1 & ⋯ & 1.73 & -0.50 & 2.04 & 1.86 & 2.05 & 1.88 & 2.12 & 1.91 & 2.05 & 1.84\\\\\n\\end{tabular}\n",
            "text/markdown": "\nA data.frame: 6 × 106\n\n| <!--/--> | Div &lt;chr&gt; | Date &lt;chr&gt; | Time &lt;chr&gt; | HomeTeam &lt;chr&gt; | AwayTeam &lt;chr&gt; | FTHG &lt;int&gt; | FTAG &lt;int&gt; | FTR &lt;chr&gt; | HTHG &lt;int&gt; | HTAG &lt;int&gt; | ⋯ ⋯ | AvgC.2.5.1 &lt;dbl&gt; | AHCh &lt;dbl&gt; | B365CAHH &lt;dbl&gt; | B365CAHA &lt;dbl&gt; | PCAHH &lt;dbl&gt; | PCAHA &lt;dbl&gt; | MaxCAHH &lt;dbl&gt; | MaxCAHA &lt;dbl&gt; | AvgCAHH &lt;dbl&gt; | AvgCAHA &lt;dbl&gt; |\n|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|\n| 1 | E0 | 09/08/2019 | 20:00 | Liverpool      | Norwich          | 4 | 1 | H | 4 | 0 | ⋯ | 3.43 | -2.25 | 1.91 | 1.99 | 1.94 | 1.98 | 1.99 | 2.07 | 1.90 | 1.99 |\n| 2 | E0 | 10/08/2019 | 12:30 | West Ham       | Man City         | 0 | 5 | A | 0 | 1 | ⋯ | 2.91 |  1.75 | 1.95 | 1.95 | 1.96 | 1.97 | 2.07 | 1.98 | 1.97 | 1.92 |\n| 3 | E0 | 10/08/2019 | 15:00 | Bournemouth    | Sheffield United | 1 | 1 | D | 0 | 0 | ⋯ | 1.92 | -0.50 | 1.95 | 1.95 | 1.98 | 1.95 | 2.00 | 1.96 | 1.96 | 1.92 |\n| 4 | E0 | 10/08/2019 | 15:00 | Burnley        | Southampton      | 3 | 0 | H | 0 | 0 | ⋯ | 1.71 |  0.00 | 1.87 | 2.03 | 1.89 | 2.03 | 1.90 | 2.07 | 1.86 | 2.02 |\n| 5 | E0 | 10/08/2019 | 15:00 | Crystal Palace | Everton          | 0 | 0 | D | 0 | 0 | ⋯ | 1.71 |  0.25 | 1.82 | 2.08 | 1.97 | 1.96 | 2.03 | 2.08 | 1.96 | 1.93 |\n| 6 | E0 | 10/08/2019 | 15:00 | Watford        | Brighton         | 0 | 3 | A | 0 | 1 | ⋯ | 1.73 | -0.50 | 2.04 | 1.86 | 2.05 | 1.88 | 2.12 | 1.91 | 2.05 | 1.84 |\n\n",
            "text/html": [
              "<table class=\"dataframe\">\n",
              "<caption>A data.frame: 6 × 106</caption>\n",
              "<thead>\n",
              "\t<tr><th></th><th scope=col>Div</th><th scope=col>Date</th><th scope=col>Time</th><th scope=col>HomeTeam</th><th scope=col>AwayTeam</th><th scope=col>FTHG</th><th scope=col>FTAG</th><th scope=col>FTR</th><th scope=col>HTHG</th><th scope=col>HTAG</th><th scope=col>⋯</th><th scope=col>AvgC.2.5.1</th><th scope=col>AHCh</th><th scope=col>B365CAHH</th><th scope=col>B365CAHA</th><th scope=col>PCAHH</th><th scope=col>PCAHA</th><th scope=col>MaxCAHH</th><th scope=col>MaxCAHA</th><th scope=col>AvgCAHH</th><th scope=col>AvgCAHA</th></tr>\n",
              "\t<tr><th></th><th scope=col>&lt;chr&gt;</th><th scope=col>&lt;chr&gt;</th><th scope=col>&lt;chr&gt;</th><th scope=col>&lt;chr&gt;</th><th scope=col>&lt;chr&gt;</th><th scope=col>&lt;int&gt;</th><th scope=col>&lt;int&gt;</th><th scope=col>&lt;chr&gt;</th><th scope=col>&lt;int&gt;</th><th scope=col>&lt;int&gt;</th><th scope=col>⋯</th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;dbl&gt;</th></tr>\n",
              "</thead>\n",
              "<tbody>\n",
              "\t<tr><th scope=row>1</th><td>E0</td><td>09/08/2019</td><td>20:00</td><td>Liverpool     </td><td>Norwich         </td><td>4</td><td>1</td><td>H</td><td>4</td><td>0</td><td>⋯</td><td>3.43</td><td>-2.25</td><td>1.91</td><td>1.99</td><td>1.94</td><td>1.98</td><td>1.99</td><td>2.07</td><td>1.90</td><td>1.99</td></tr>\n",
              "\t<tr><th scope=row>2</th><td>E0</td><td>10/08/2019</td><td>12:30</td><td>West Ham      </td><td>Man City        </td><td>0</td><td>5</td><td>A</td><td>0</td><td>1</td><td>⋯</td><td>2.91</td><td> 1.75</td><td>1.95</td><td>1.95</td><td>1.96</td><td>1.97</td><td>2.07</td><td>1.98</td><td>1.97</td><td>1.92</td></tr>\n",
              "\t<tr><th scope=row>3</th><td>E0</td><td>10/08/2019</td><td>15:00</td><td>Bournemouth   </td><td>Sheffield United</td><td>1</td><td>1</td><td>D</td><td>0</td><td>0</td><td>⋯</td><td>1.92</td><td>-0.50</td><td>1.95</td><td>1.95</td><td>1.98</td><td>1.95</td><td>2.00</td><td>1.96</td><td>1.96</td><td>1.92</td></tr>\n",
              "\t<tr><th scope=row>4</th><td>E0</td><td>10/08/2019</td><td>15:00</td><td>Burnley       </td><td>Southampton     </td><td>3</td><td>0</td><td>H</td><td>0</td><td>0</td><td>⋯</td><td>1.71</td><td> 0.00</td><td>1.87</td><td>2.03</td><td>1.89</td><td>2.03</td><td>1.90</td><td>2.07</td><td>1.86</td><td>2.02</td></tr>\n",
              "\t<tr><th scope=row>5</th><td>E0</td><td>10/08/2019</td><td>15:00</td><td>Crystal Palace</td><td>Everton         </td><td>0</td><td>0</td><td>D</td><td>0</td><td>0</td><td>⋯</td><td>1.71</td><td> 0.25</td><td>1.82</td><td>2.08</td><td>1.97</td><td>1.96</td><td>2.03</td><td>2.08</td><td>1.96</td><td>1.93</td></tr>\n",
              "\t<tr><th scope=row>6</th><td>E0</td><td>10/08/2019</td><td>15:00</td><td>Watford       </td><td>Brighton        </td><td>0</td><td>3</td><td>A</td><td>0</td><td>1</td><td>⋯</td><td>1.73</td><td>-0.50</td><td>2.04</td><td>1.86</td><td>2.05</td><td>1.88</td><td>2.12</td><td>1.91</td><td>2.05</td><td>1.84</td></tr>\n",
              "</tbody>\n",
              "</table>\n"
            ]
          },
          "metadata": {
            "tags": []
          }
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "GRTdVpljpM93",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 34
        },
        "outputId": "b687e147-7e51-45eb-a041-71d0d287bc59"
      },
      "source": [
        "nrow(data)"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "[1] 380"
            ],
            "text/latex": "380",
            "text/markdown": "380",
            "text/html": [
              "380"
            ]
          },
          "metadata": {
            "tags": []
          }
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "N_l0raSxSjZ6",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 67
        },
        "outputId": "31033269-3a17-421e-e1a3-24800c859ebf"
      },
      "source": [
        "FTHG_freq_table = table(data$FTHG)\n",
        "FTHG_freq_table"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "\n",
              "  0   1   2   3   4   5   8 \n",
              " 84 125  99  47  16   8   1 "
            ]
          },
          "metadata": {
            "tags": []
          }
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "O97vFYe1SnMl",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 437
        },
        "outputId": "cd02c781-0bbe-4367-cf88-e81f212695bc"
      },
      "source": [
        "barplot(FTHG_freq_table,col = 'orange'\n",
        "         ,main='Goal by Home Team in Full Time'\n",
        "         ,ylab='Frequency: No of Matches'\n",
        "         ,xlab='No of Goal in a Match')"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAA0gAAANICAIAAAByhViMAAAACXBIWXMAABJ0AAASdAHeZh94\nAAAgAElEQVR4nOzdeUBVZeL/8ed6ZZNVk1VZFBiUHHNLTSQrDRTUzG/qOKUWYSFqRkpluYBL\nfttmLARz0iQxx1RcErWycEdDDRm3cUVzQ1xQAdm9vz/Od+6PYTsX5G7H9+uvc5/zcO/nLo2f\nOatKo9EIAAAAmL9mxg4AAACApkGxAwAAUAiKHQAAgEJQ7AAAABSCYgcAAKAQFDsAAACFoNgB\nAAAoBMUOAABAISh2AAAACkGxAwAAUAiKHQAAgEJQ7AAAABSCYgcAAKAQFDsAAACFoNgBAAAo\nBMUOAABAISh2AAAACkGxAwAAUAiKHQAAgEJQ7AAAABSCYgcAAKAQFDsAAACFoNgBAAAoBMUO\nAABAISh2AAAACkGxAwAAUAiKHQAAgEJQ7AAAABSCYgcAAKAQFDsAAACFoNgBAAAoBMUOAABA\nISh2AAAACkGxAwAAUAiKHQAAgEJQ7AAAABSCYgcAAKAQFDsAAACFoNgBAAAoBMUOMDnr1q1T\nqVQqlcra2rqpZkLxzO7HUGtgPb0Ls/twgEaj2AHVbdmyJTIyMjAw8LHHHrOwsHB0dOzSpcuE\nCRP27dtn7GhGs3TpUtV/3Llzp+aEAQMGSGsHDhxo+HiNpo0t66uvvjJ2WOOo+tXXpXfv3gZO\nxRcH1KW5sQMAJuTMmTOjR48+fPhw1cF79+5lZ2dnZ2d/9dVXw4YNS05OdnR0NFZCoC5dunT5\n+9//LoRo3pz/Ya+ODwePDn7iwP85ffp079698/PzpYdqtfrxxx93cHC4dOnSH3/8odFohBAb\nN24MDw/fuXMn/zwoQ1BQkJ2dnfbhxYsXjxw5Ii2HhYVZWFhoV7Vr187Q4RrIz8/v7bff1utL\nVPtMtAICAvT6ujU19IszwIcDmAj+cQKEEEKj0YwaNUrb6t56663Zs2e3atVKenj8+PHJkyfv\n2LFDCLFv374vvvhi6tSpRsuKphMfH1/14dKlS8ePHy8tf/fdd05OTsYIZbpM5zPhiwPqwjF2\ngBBCpKWlaf8ff2xs7BdffKFtdUKIxx9/fNu2bT169GjRosWQIUO8vb2r/fm9e/c++uij3r17\nt2rVytLS0tXVNTQ0dPny5ZWVldVmajSa1atXh4aGuri4WFhYODg49OrVKyEhoebMBmnWrJkQ\nYsOGDU8//bSjo6ODg8Ozzz67c+dO7YSIiAjpkKNOnTpV+9sVK1Zojyuv9fi5h6HjJ/PNN99I\nGZ5++mkhxLp167p27dqiRQsfH58ZM2aUl5cLIU6ePDl06NCWLVva2dmFhIQcP3682mudOnUq\nKirK39/f2trawcHhySef/PLLLysqKprw7ejyErp8xU3yfqup9fyAai+0Z8+e559/XnrO4ODg\nX375pQk/HCHEjBkzaj3qbuPGjdK4sTZ16/LhmOyvDmgYDQCNZsyYMdJ/EY6OjkVFRbXOuX79\nenFxcc3xI0eOtGnTptb/vp566qmbN29Wnfzyyy/XOnPw4MEPHjyQ5qxdu1YatLKyqj+2dqad\nnd3ixYurPWezZs3S0tKkmVXP/Dh16lTVJxkxYoQ0PmLEiLpe6Ouvv9b+eX5+fs0J/fv3l9aG\nhoY24pNZvXq1NN6pU6c1a9aoVKqqkydOnHj+/PnWrVtXHXR2dr5z5472GVJTU2s94fG5556r\n9VtrxNvU8SV0+Yof/v3WVOvPpuoL/fTTT5aWllWfU61W//LLLw/zmVTz4YcfSjN79epVdXzD\nhg3aV6w/sO4//gaFlP1wjPirA5oWW+wAIYQ4cOCAtDBgwIAWLVrUOsfFxaXm/4jn5+eHh4df\nuXJFCNGuXbvFixdv3Ljxvffek7ZM7N+/f9y4cdrJaWlp3333nRCiWbNmSUlJR48eXbZsmTQz\nLS1t3bp1jc5fXl7+7rvvRkZG/vOf//zkk0+k0zsePHgwYcKEsrIyIUSfPn0CAwOlyevXr6/6\nhz///LO0XDXqw2vQJ6NWq6WFW7duvfPOO0OGDJk0aZL2JJWlS5dOnDjR0tJyypQpvXr1kgZv\n3LixfPlyaTknJ+eVV14pKSkRQsTGxp46dergwYP9+vUTQqSnp8+ZM+fh346OL6HjV/yQ71d3\nVV/ozTffDAwMnD59ekhIiDRYWVn5wQcfNObjUATT/9UBjWHsZgmYBFtbW+m/iOnTpzfoD7XH\n+jg4OFy5ckU7npKSov2v7NChQ9LgokWLwsPDw8PDp0yZop05dOhQadrYsWOlkUZssRNCREVF\nacerdsStW7dKg59//rk08uSTT2pn/vrrr9Kgi4tLeXl5XS9UdYtI/bRb7Br0yVR9I9oNh6mp\nqdpBa2vrc+fOaTSa0tJS7aH6Q4YMkWZOmjRJGnnmmWe0L3Tjxg3p+Hp7e3sdN5/Us+FHx5do\n6FfcuPdbq/o3gAkhgoODtZ+DdrNis2bNysrKGveZ1GReW+xM5FcHNC222AFCCHH//n1pQdvw\nJBUVFbVeHEt7QJ5269fw4cM9PDy0fzh69OiWLVtKy1u3bpUWJk6cmJaWlpaWtnDhQu1M7RF7\nubm5D/MWJk6cqF0eNmyYg4ODtKzdGDl27FhpT9zBgwcvXbokDaalpUkLL7/8ctMe/9SgT6aq\n2NhYaWHIkCFWVlbS8osvvti+fXshhKWl5ZAhQ6RB7bvYtm2btBAUFFTyH3Z2dt26dRNCFBQU\nPPw1CHV8iUZ8xY14v40THx+v3eT82muvSQsPHjy4evWqjs/w8ssvD6vNqlWrHiaYKTDNXx3Q\nCJwVCwghhJ2dXUFBgRDi3r17uv+VRqPRHkz95z//ueoqtVrdoUOH/fv3CyFOnjypHf/ll18W\nLVp09OjRK1eulJaWVv2Thzl/wsLCQrunVXp1Pz+/33//XQjxxx9/SIOtW7ceNmzYmjVrhBAb\nNmx46623RJViN3bsWB1fy9PTUzpXo6rr169L+6QkjfhktLSnd1hYWDg7O1++fFkI0blzZ+0E\nd3d3aUH6yjQaTU5OjjQyf/78+fPn13zO48ePa48CbIQGvURDv+KGvt9G69Gjh3a56glAuj9t\nrUVcVHkL5ssEf3VA41DsACGE8PLykorIsWPHqo43a9as6pawpKQkjUajfXj//n3t6W9Vr6ol\n0W780/7DuXjx4ujoaO3a9u3bW1paXrly5ebNmw+Z39bWtlrZ0h4pWFxcrB2MjIyUit369evf\neuutM2fOnDlzRgjRuXPnLl266Pha//rXv2peTmLAgAHavbqi4Z+MlpWVlY2NTdWH0kLVV6x2\nBsD9+/cfPHhQf+bbt2/XP6F+ur9EQ7/iRrzfxrGysrK3t9c+5M5aVZnmrw5oHIodIIQQvXv3\nlopdenp6fn6+dl9hs2bNFi1apJ321VdfVd3o0qJFi+bNm0sNpuamPm1rkQ7HLiwsnDZtmjTy\n17/+denSpdK/JRMnTkxKSnrI/NLRPFVP6ysqKpIWqu5cHjBggI+Pz4ULF/bu3Xvjxo3NmzdL\n47pvrtNRgz6Zh38ttVotfS8LFy6cMmXKwz9n415Cr1+xKcjPz9fxEnHVNlXm5eXpJ5HRGOBX\nBzQOx9gBQggxevRoaaGkpOTdd9+tdU7N/4Ne9bJw2qPuJOXl5SdOnJCWpR06WVlZ2iP5YmNj\ntVsITp069fD5S0tLqz5PRUXF2bNnpWUfH5+qgSMiIoQQlZWVmzZtkvbDqtXquq7Q0WgN+mQe\n/rV8fX2lZe3esaal40vo9Ss2C9qts1evXq26bVs6KkBJDPCrAxqHYgcIIUT//v2feuopaXnp\n0qWvvfba9evXtWvLy8t//PHH4ODgmsdIDR8+XFrYuHGjdFyOJDk5WdoupVKphg0bJv57G4Z0\n7VMhxIkTJ6QbWlQdbJx//OMf2uX169drt4r17du36rSIiAjpKg/Jycl79+4VQoSEhLi5uT3M\nS9dK90/m4YWGhkoLa9eu1VarysrKl19+OSIiYvr06dJVV/T9Evr+ik1f27ZtpYW8vLwtW7ZI\ny//+97+rngqtGAb41QGNQLED/k9KSoqLi4u0nJyc3KZNmy5duvTv379Hjx4tW7YcNGhQRkaG\ntHbEiBHaax9MnjxZ+sesqKjo2WefTUpKWr9+/bvvvjt58mRpwuuvv/6nP/1JCBEYGKg9DG7u\n3LnHjh3bvHlzSEiIv7+/NHjkyJEDBw40aKeVdiOiWq1OTEx85513fvjhh7/97W/a2ysFBARI\nV9XXatOmzcCBA4UQ+/btk3pG016+Tkv3T+bhvf3229LmsatXr4aEhGzduvXnn38eOXLkqlWr\nli9fvnnz5odvrrq8hD6+YvPSp08f7fEAf/nLX6KiomJjY4OCgnQ/gtOMGOBXBzSG8a60Apic\n8+fPa69EWitvb+9ly5ZV+6t67q8wfPjwqtey0nYaLQ8Pj5ycnKpXA5k9e7bul/JauXKlNNPF\nxaXmeXnW1tb79u2r+Vfai4oJIRwdHXW52laT33mi2idT11vW7u1avHixdjAhIUEa9PX1rfoM\n2mPeq2rTps3Jkydl36Aub1OXl3jIr1j391tTgy4LV/XKKUePHm30Z1JTVFRUtU/A399f2jYs\nhFCpVPXfYcUo17Ez7q8OaFpssQP+v3bt2h04cCAtLS0iIiIgIMDJycnS0tLDw6N79+7S9clO\nnz4tHaNW1RNPPHHixIl58+Y9+eSTjo6OFhYW7u7uL7744g8//FDtjkOff/75xx9/HBgYaGNj\n06ZNm8jIyIMHD/r4+Hz77bcBAQHNmzdv27Ztx44ddQ+svcJIq1atPvjgg+XLl3fp0sXa2rpl\ny5ZDhw49cOBAnz59av7V4MGDnZ2dpeURI0bo7wRJ3T+Zh/fSSy9lZWVFRES0a9fOysqqRYsW\nnTp1+vDDD48ePdqhQweDvUSTf8VmZ9GiRQsWLPD397e0tGzTps0bb7yxd+9e7YGeGo2m6mna\n5s4AvzqgoVSaKse3AngUXL9+3cvLS7rV2N69e4OCgoydCADQNNhiBzxyZs6cKbW67t270+oA\nQEm4jh3wqFi5cuX169f37dunPcZu3rx5xo0EAGha7IoFHhXPPPPMrl27tA/feOONJUuWGDEP\nAKDJscUOeFQ4OztbWVlVVlb6+vq+8cYbb7/9trETAQCaGFvsAAAAFIKTJwAAABSCYgcAAKAQ\nFDsAAACFoNgBAAAoBMUOAABAISh2AAAACkGxAwAAUAiKHQAAgEJQ7AAAABSCYgcAAKAQFDsA\nAACFoNgBAAAoBMUOAABAISh2AAAACkGxAwAAUAiKHQAAgEJQ7AAAABSCYgcAAKAQFDsAAACF\noNgBAAAoBMUOAABAISh2AAAACkGxAwAAUAiKHQAAgEJQ7AAAABSCYgcAAKAQFDsAAACFoNgB\nAAAoBMUOAABAISh2AAAACkGxAwAAUAiKHQAAgEJQ7AAAABSiubEDAP/n/Pnz+fn5xk5hCCqV\n6oknnlCr1cYOAgBQGpVGozF2BkAIIZycnO7evWvsFAaycePGF154wdgpAABKwxY7mIqysrKt\n74pBTxg7h/65TBBlZWXGTgEAUCCOsQMAAFAIih0AAIBCUOwAAAAUgmIHAACgEBQ7AAAAhaDY\nAQAAKIT5Xe5Eo9Hk5OScP3++oKBACOHo6Ojv7+/p6WnsXAAAAEZmTsUuPz9//vz5KSkpeXl5\n1VZ5eXlFRkZOmzbNxsbGKNkAAACMzmyK3bVr14KCgnJycvz9/cPCwry9vW1tbYUQ9+7dO3fu\n3K5du2bNmpWamrpjx46WLVsaOywAAIARmE2xmzlz5uXLl9esWTNixIiaaysrK5csWTJp0qT4\n+PiFCxcaPh4AAIDRmc3JE1u2bBkzZkytrU4IoVaro6OjR44cuX79egMHAwAAMBFmU+xu3brl\n6+tb/5yOHTtev37dMHkAAABMjdkUOw8Pj+zs7PrnZGVleXh4GCYPAACAqTGbYjds2LC1a9d+\n9tlnpaWlNdcWFRXNnj1706ZNo0aNMnw2AAAAU2A2J0/ExcXt2bMnNjZ2zpw5PXv29PT0tLOz\n02g0hYWFFy9ezMzMvH//fnBw8IwZM4ydFAAAwDjMptg5OTnt378/MTFxxYoVO3furKys1K6y\nsLDo3r17RERERESEWq02YkgAAAAjMptiJ4SwtLSMiYmJiYkpKSm5dOmSdOcJBwcHLy8vS0tL\nY6cDAAAwMnMqdlrW1tb+/v41x/Pz8+/evevj42PwRAAAAMZnNidPCCH+9a9/hYeH+/j4BAcH\nJyUlVd0bK/n444/btWtnlGwAAABGZzZb7Pbt29e/f//S0tIWLVpcvXp17969a9as2bBhAzcQ\nAwAAkJjNFrsFCxY8ePBgw4YNhYWFBQUFf/vb3zIyMkJDQ4uKiowdDQAAwCSYTbH717/+NWrU\nqGHDhqlUKisrq5iYmB9//DE7O3vkyJE198kCAAA8gsxmV2xubm779u2rjjz33HNLly4dO3bs\nO++888UXXzTuafPy8iIjI4uLi+uZU1ZWlpube/LkyWbNzKYHAwCAR5DZFDtXV9cjR45UGxwz\nZszJkycXLFjQtm3b2NjYRjytjY3NE088UV5eXs+cK1eu7N69u6KigouqAAAAU2Y2xW748OEJ\nCQmLFi168803LSwstOPz58+/evXqu+++e/Xq1Ubsk7W3t587d279czIyMlauXNngxAAAAIZl\nNsVu1qxZGzdunDx58qZNm7Zv364dV6lUy5cvd3R0XLhwoRHjAQAAGJ3ZHDT22GOPHT58ODo6\nulOnTtVWqVSqL774IjU11dfX1yjZAAAATIHZbLETQrRu3ToxMbGutcOHDx8+fLgh8wAAAJgU\ns9lip4tbt26dPXvW2CkAAACMQ1HF7tNPP631HrIAAACPAkUVOwAAgEcZxQ4AAEAhzObkiR49\nesjOuXLligGSAAAAmCazKXZZWVlCiKqXJq6poqLCUHEAAABMjtnsio2NjbW1tT127FhJ3aZN\nm2bsmAAAAEZjNsVu7ty5fn5+o0ePrv++rgAAAI8ssyl2FhYW33333fHjxz/44ANjZwEAADBF\nZnOMnRCiY8eOubm59RxIN2jQICcnJ0NGAgAAMB3mVOyEEA4ODvWs7devX79+/QwWBgAAwKSY\nza5YAAAA1I9iBwAAoBAUOwAAAIWg2AEAACgExQ4AAEAhKHYAAAAKQbEDAABQCIodAACAQlDs\nAAAAFIJiBwAAoBAUOwAAAIWg2AEAACgExQ4AAEAhKHYAAAAKQbEDAABQCIodAACAQlDsAAAA\nFIJiBwAAoBAUOwAAAIWg2AEAACgExQ4AAEAhKHYAAAAKQbEDAABQCIodAACAQlDsAAAAFIJi\nBwAAoBAUOwAAAIWg2AEAACgExQ4AAEAhKHYAAAAKQbEDAABQCIodAACAQlDsAAAAFIJiBwAA\noBAUOwAAAIWg2AEAACgExQ4AAEAhKHYAAAAKQbEDAABQCIodAACAQlDsAAAAFIJiBwAAoBAU\nOwAAAIWg2AEAACgExQ4AAEAhKHYAAAAKQbEDAABQCIodAACAQlDsAAAAFIJiBwAAoBAUOwAA\nAIWg2AEAACgExQ4AAEAhKHYAAAAKQbEDAABQCIodAACAQlDsAAAAFIJiBwAAoBAUOwAAAIWg\n2AEAACgExQ4AAEAhKHYAAAAKQbEDAABQCIodAACAQlDsAAAAFIJiBwAAoBDNjR2gwTQaTU5O\nzvnz5wsKCoQQjo6O/v7+np6exs4FAABgZOZU7PLz8+fPn5+SkpKXl1dtlZeXV2Rk5LRp02xs\nbIySDQAAwOjMpthdu3YtKCgoJyfH398/LCzM29vb1tZWCHHv3r1z587t2rVr1qxZqampO3bs\naNmypbHDAgAAGIHZFLuZM2devnx5zZo1I0aMqLm2srJyyZIlkyZNio+PX7hwoeHjAQAAGJ3Z\nnDyxZcuWMWPG1NrqhBBqtTo6OnrkyJHr1683cDAAAAATYTbF7tatW76+vvXP6dix4/Xr1w2T\nBwAAwNSYTbHz8PDIzs6uf05WVpaHh4dh8gAAAJgasyl2w4YNW7t27WeffVZaWlpzbVFR0ezZ\nszdt2jRq1CjDZwMAADAFZnPyRFxc3J49e2JjY+fMmdOzZ09PT087OzuNRlNYWHjx4sXMzMz7\n9+8HBwfPmDHD2EkBAACMw2yKnZOT0/79+xMTE1esWLFz587KykrtKgsLi+7du0dERERERKjV\naiOGBAAAMCKzKXZCCEtLy5iYmJiYmJKSkkuXLkl3nnBwcPDy8rK0tDR2OgAAACMzp2In0Wg0\nV69evXjxovaWYlZWVtxSDAAAwJyKHbcUAwAAqIfZFDtuKQYAAFA/syl23FIMAACgfmZT7HS5\npdju3bvXr1/foGKXn58/c+bM8vLyeuZwNwsAAGAWzKbY6XhLsQ0bNjToaTUaTUFBQXFxcT1z\npLM0AAAATJzZFDs93VKsVatW3377bf1zMjIy0tPTG/S0AAAAhsctxQAAABTCbLbYcUsxAACA\n+plNseOWYng0nT9/fvLkyWVlZcYOYgje3t5Lly41dgoAMGNmU+wEtxTDI+no0aO/bt/69kBj\n59C/izdFcrKaYgcAD8Ocip2WtbW1v79/zfFbt27l5+f7+fkZPhKgP9YW4n//YuwQ+vfrcbE2\n09ghAMDMmc3JE7r49NNPay18AAAAjwJFFTsAAIBHGcUOAABAIczmGLsePXrIzrly5YoBkgAA\nAJgmsyl2WVlZQggLC4t65lRUVBgqDgAAgMkxm12xsbGxtra2x44dK6nbtGnTjB0TAADAaMym\n2M2dO9fPz2/06NHl5eXGzgIAAGCKzKbYWVhYfPfdd8ePH//ggw+MnQUAAMAUmc0xdkKIjh07\n5ubm1nMg3aBBg5ycnAwZCQAAwHSYU7ETQjg4ONSztl+/fv369TNYGAAAAJNiNrtiAQAAUD+K\nHQAAgEJQ7AAAABSCYgcAAKAQFDsAAACFoNgBAAAoBMUOAABAISh2AAAACkGxAwAAUAiKHQAA\ngEJQ7AAAABSCYgcAAKAQFDsAAACFoNgBAAAoBMUOAABAISh2AAAACkGxAwAAUAiKHQAAgEJQ\n7AAAABSCYgcAAKAQFDsAAACFoNgBAAAoBMUOAABAISh2AAAACkGxAwAAUAiKHQAAgEJQ7AAA\nABSCYgcAAKAQFDsAAACFoNgBAAAoBMUOAABAISh2AAAACkGxAwAAUAiKHQAAgEJQ7AAAABSC\nYgcAAKAQFDsAAACFoNgBAAAoBMUOAABAISh2AAAACkGxAwAAUAiKHQAAgEJQ7AAAABSCYgcA\nAKAQFDsAAACFoNgBAAAoBMUOAABAISh2AAAACkGxAwAAUAiKHQAAgEJQ7AAAABSCYgcAAKAQ\nuha7yspK7XJpaelvv/2WlZWl0Wj0kwoAAAANJl/sKisrJ06c+Je//EV6eOHChcDAwN69e3fr\n1u3pp58uLCzUc0IAAADoRL7Yffrpp0lJSV5eXtLDiRMn5uTkTJgwITo6OiMjY9GiRXpOCAAA\nAJ00l53x3XffDR8+/PPPPxdCXLlyZdu2bREREUlJSUKIkpKS77///v3339d7TAAAAMiR32J3\n4cKFkJAQafmnn37SaDSjR4+WHnbv3v3ChQv6CwcAAADdyRc7lUqlXf7ll19sbW2Dg4OlhxqN\npry8XF/RAAAA0BDyxc7b23v37t1CiOvXr2/evDkkJMTS0lJalZ2d3bZtW/0GBAAAgG7ki91f\n//rXVatW9enTp1u3boWFhVOmTJHGV6xY8e233w4dOlTPCQEAAKAT+ZMnYmJiTp8+/f3331ta\nWn755Zf9+vWTxt9///2AgIDp06frOaHSlJSUXL161dgpDMTNza1FixbGTgEAwKNCvthZW1sv\nX758+fLl1cbXr1/fo0eP5s3lnwFVTZ06VTqn+FEwduzYb7/91tgpAAB4VDSglhUUFPzxxx9t\n2rRxcnISQvTu3VtvqZTs/v37/9NTfDLa2Dn0Ly5V3L9/39gpAAB4hOhU7Hbt2jV16tTDhw8L\nIbZt2zZw4EAhxNChQ6dMmdK/f3/9BlQie2vR3sXYIfTPsYUoNnYGAAAeKfInT2RmZoaEhJw+\nfTo0NFQ7eOPGjYMHD4aFhUltDwAAAEYnX+zmzJnj5uZ24sSJ5ORk7aCzs3N2drabm9vcuXP1\nmA4AAAA6ky92Bw4cmDBhQs3r1bm4uERFRUmXuAMAAIDRyRe7u3fvenp61rrK3d29sLCwqSMB\nAACgMeSLnZub28mTJ2tdtXv3bg8Pj6aOBAAAgMaQL3ZhYWFJSUm///571cH8/PwPP/xw+fLl\n4eHhessGAACABpC/3El8fPy2bdt69erVuXNnIcT06dOnT59+8uTJ0tJSLy+vWbNm6T/kf9Fo\nNDk5OefPny8oKBBCODo6+vv717WzGAAA4NEhX+zc3NwOHToUFxe3Zs0aIcSRI0eEEK1bt46I\niIiLi3NxMdwF2fLz8+fPn5+SkpKXl1dtlZeXV2Rk5LRp02xsbAyWBwAAwKTodIFiFxeXpKSk\nxMTEvLy8goICe3t7V1dXfSer5tq1a0FBQTk5Of7+/mFhYd7e3ra2tkKIe/funTt3bteuXbNm\nzUpNTd2xY0fLli0NnA0AAMAUNOCWYiqVytXV1fCVTjJz5szLly+vWbNmxIgRNddWVlYuWbJk\n0qRJ8fHxCxcuNHw8AAAAo5MvdhqNZtWqVd9///3Vq1dLSkpqTjh27JgeglW3ZcuWMWPG1Nrq\nhBBqtTo6Onr37t3r16+n2AEAgEeTfLGbO3fu7NmzhRBqtdrOzk7/kWp369YtX1/f+ud07Nhx\nw4YNhskDAABgauSL3dKlS728vDZu3NilSxeVSmWATLXy8PDIzs6uf05WVhbX1QMAAI8s+evY\n5ebmTpo0qWvXrkZsdUKIYcOGrV279rPPPistLa25tqioaPbs2Zs2bRo1apThswEAAJgC+S12\n7u7uGo3GAFHqFxcXt2fPntjY2Dlz5vTs2dPT09POzk6j0RQWFl68eDEzM/P+/S+3Ew4AACAA\nSURBVPvBwcEzZswwdlIAAADjkC92UVFRa9eujYmJsbCwMECgujg5Oe3fvz8xMXHFihU7d+6s\nrKzUrrKwsOjevXtERERERIRarTZiSAAAACOqvdidPXtWuzxy5MjTp08/99xzU6dO9ff3t7Ky\nqjbZz89PjwGrsLS0jImJiYmJKSkpuXTpknTnCQcHBy8vL0tLy8Y9Z2Fh4aefflrr7l2tK1eu\nNO7JAQAADKn2Yufv719zcO/evbVONvyOWmtra23C8vLy06dPl5SUdOrUqWbplFVUVHTo0KGy\nsrJ65ty9e1cY420CAAA0SO3F7vXXXzdwDl2kp6fPnTv3woULHTt2nD17dq9evX766aeIiIir\nV68KIRwcHBYsWBAdHd2g53R1dd2yZUv9czIyMoKCgox77ggAAICs2ovd0qVLDZxD1v79+0ND\nQysqKhwcHH766afdu3dv37595MiRjo6OY8eOLSkp2b59+8SJE9u3bz9w4EBjhwUAADAC+cud\nSI4fP37z5s2qD7OysvQTqXYLFixo3bp1dnb23bt3c3Nze/bs+Ze//KVdu3anTp369ttvv//+\n+3Pnzvn4+HzxxReGTAUAAGA65ItdeXn566+/3qlTp6q3DtuxY0e3bt1ee+21qmen6lVGRsbE\niRM7d+4shHB2dv7kk0/++OOPmJgYGxsbaULLli0jIyMzMzMNkwcAAMDUyBe7hISEb775Jjw8\n3NvbWzv4/PPPjxo1Kjk5edGiRfqM9//dvXu3aoA2bdoIIZydnavOcXd3v3fvnmHyAAAAmBr5\nYpecnDx48OC0tLR27dppBwMCAlavXh0WFmawYvfYY4+dO3dO+/D06dPivy/LIoQ4d+7cY489\nZpg8AAAApka+2J09e/bZZ5+tddUzzzxz8eLFpo5Uu2effTYhIWHHjh1lZWVHjx6dPHlyx44d\nP//8c+1F5k6ePLl48eLg4GDD5AEAADA18neecHBwuHDhQq2rLly40KpVqyZOVIfZs2dv2bLl\nueeekx62atVq7969gwYN+tOf/tSrV6+SkpKDBw9qNJrY2FjD5AEAADA18lvswsPDly1btnXr\n1qqD5eXlX3/99T/+8Y+QkBC9ZfsvHTp0yMjIGD16dK9evV599dWMjIyOHTumpaV16tRp586d\n+/fv9/LySk1N7dmzp2HyAAAAmBr5LXbz5s3btm1beHi4l5dXQECAlZXVnTt3Tpw4cfv2bXd3\n93nz5hkgpaRTp06rVq2qNvLbb78VFhYWFxdXO5ECAADgUSO/xc7d3T0rKysqKqqoqGj79u1p\naWl79+5Vq9Xjx48/ePCgl5eXAVLWz87OjlYHAAAgv8VOCOHq6rp48eKkpKRr164VFxe7ubnZ\n2trqOxkAAAAaRH6L3d69e2/fvi2EUKlUHh4evr6+2laXmZmZmpqq34AAAADQjXyxCw4O3r17\nd62r9uzZM378+KaOBAAAgMaoc1fs2bNntZf/zcrKsra2rjahuLh4zZo1paWlekwHAAAAndVZ\n7NatWzd9+nRpec6cOXVNe+mll5o+FAAAABquzmL3/vvvjxs37uDBgy+88MKYMWMCAwOrTVCr\n1e3btx86dKieEwIAAEAn9Z0V6+7uPnTo0PDw8Ojo6N69e9ecUFRUdOvWLTc3N73FAwAAgK7k\nT55IS0urtdUJITZt2tStW7emjgQAAIDG0Ok6djdv3ly9evWFCxcqKiq0gyUlJWlpaYWFhXrL\nBgAAgAaQL3YXLlzo2bPnjRs3avnj5s1nzpyph1QAAABoMPliN2PGjJKSkkWLFnXs2LF///5L\nly5t27btzp07U1JSli1bFhoaaoCUAAAAkCVf7Pbs2TNx4sSJEyeWlJQIIR5//PHevXuHhoaO\nGjWqf//+P/zwQ1BQkP5zAgAAQIb8yRPXrl1r3769EKJZs2ZCiLKyMmm8S5cuEydOnD17tl7z\nAQAAQEfyxc7e3v769etCCEtLSzs7u/Pnz2tXBQYGHjp0SI/pAAAAoDOd7hX71Vdf7dy5Uwjx\n5z//OTExUXsmbHp6upWVlV7zAQAAQEfyxe6DDz64devWtGnThBDjx48/dOhQYGDg8OHDu3bt\n+vXXXz///PP6DwkAAAB58idP9OzZc+/evZmZmUKIV1999cyZMwsXLtywYYNKpRo6dOjChQv1\nHxIAAADydLpAcffu3bt37y6EUKlUH3300axZs3Jzc11dXW1sbPQcDwAAALrSqdhVY21t7ePj\n09RJAAAA8FDqLHZV7x4m8xTNG9MOAQAA0LTq7GQWFhY6PoVGo2miMAAAAGi8+ja2qdXqzp07\nBwQEUN0AAABMX53F7u233161alVWVtbt27dHjhw5bty4xx9/3JDJAAAA0CB1Xsfu73//+5Ur\nVzZv3vzkk09++eWXnTp16tGjR0JCws2bNw2ZDwAAADqq7wLFzZs3Hzx48Nq1a3Nzc7/66itL\nS8u33nrLw8PjxRdf3LhxY3l5ucFSAgAAQJb8nSeEEE5OTm+++WZGRsaZM2fef//9I0eOvPji\nix4eHm+99dbhw4f1HREAAAC60KnYafn5+c2ZM+f8+fO7du0KDAxMSEjo0aOHnpIBAACgQRp2\nCboHDx5s3749OTl58+bNRUVFfn5+Y8eO1VMyAAAANIiuxe7UqVPJyckpKSlXrlyxt7cfOXLk\na6+9FhwcrNdwAAAA0J1Msbtz587333+fnJx84MABlUrVr1+/jz766KWXXmrRooVh8gEAAEBH\ndRa7n376KTk5eePGjSUlJe3atYuLixs3bhy3iAUAADBZdRa7gQMHqtXqnj17Dhw4sHv37hqN\n5tixY8eOHas5c/DgwfpMCAAAAJ3Utyu2srJy//79+/fvr/8puOEYAACAKaiz2KWkpBgyBwAA\nAB5SncXulVdeMWQOAAAAPKSGXaAYAAAAJotiBwAAoBAUOwAAAIWg2AEAACgExQ4AAEAhKHYA\nAAAKUXuxmzp1anp6urQcFRWVnZ1twEgAAABojNqL3cKFCzMzM6XlJUuW5OTkGDASAAAAGqP2\nCxS7urp+/PHHly5dsre3F0KkpKQcOHCgrqf43//9X32lAwAAgM5qL3affPLJ+PHjk5KSpIfr\n16+v5ykodgAAAKag9mL3yiuvDB48+OzZsyUlJcHBwR999FFwcLCBkwEAAKBB6rxXrJOTU48e\nPYQQoaGhzzzzzFNPPWXAVAAAAGiwOoud1o8//igt3Lx588yZM0VFRfb29gEBAU5OTnrOBgAA\ngAbQ6Tp2e/fu7d27t7Ozc58+fZ5//vnevXu3atVqwIABx44d03c+AAAA6Eh+i11mZuaAAQMq\nKir69u0bEBBgY2NTVFR04sSJ9PT0oKCgzMzMgIAAAwQFAABA/eSL3bx585ydnbdv396hQ4eq\n41lZWQMHDoyPj1+1apXe4gEAAEBX8rtiMzIyoqOjq7U6IUTXrl2jo6O1N6gAAACAcckXu7t3\n77Zt27bWVT4+Prdv327qSAAAAGgM+WLn4uJy8uTJWledOHHCxcWlqSMBAACgMeSLXUhISEJC\nwqZNmzQajXZQo9Fs2LAhMTFx0KBB+owHAAAAXcmfPDF79uytW7cOGzbMzc0tMDDQ1tZWOis2\nNzfX3d199uzZBkgJAAAAWfJb7Hx8fA4dOjRu3Lji4uL09PTNmzenp6eXlZVFRkYePny4rsPv\nAAAAYGDyW+yEEJ6ensnJyRqNJjc3t6ioyM7Ozs3NTd/JAAAA0CA6FTuJSqVyd3fXXxQAAAA8\nDJ1uKQYAAADTR7EDAABQCIodAACAQlDsAAAAFIJiBwAAoBANOCv25s2bZ86cKSoqsre3DwgI\ncHJy0l8sAAAANJROW+z27t3bu3dvZ2fnPn36PP/88717927VqtWAAQOOHTum73wAAADQkfwW\nu8zMzAEDBlRUVPTt2zcgIMDGxka6pVh6enpQUFBmZmZAQIABggIAAKB+8sVu3rx5zs7O27dv\n79ChQ9XxrKysgQMHxsfHr1q1Sm/xAAAAoCv5XbEZGRnR0dHVWp0QomvXrtHR0enp6foJBgAA\ngIaRL3Z3795t27Ztrat8fHxu377d1JEAAADQGPLFzsXF5eTJk7WuOnHihIuLS1NHAgAAQGPI\nF7uQkJCEhIRNmzZpNBrtoEaj2bBhQ2Ji4qBBg/QZDwAAALqSP3li9uzZW7duHTZsmJubW2Bg\noK2trXRWbG5urru7++zZsw2QEgAAALLkt9j5+PgcOnRo3LhxxcXF6enpmzdvTk9PLysri4yM\nPHz4cF2H3wEAAMDAdLrzhKenZ3Jyskajyc3NLSoqsrOzc3Nz03cyAAAANEgDbimmUqnc3d31\nF0VHGo0mJyfn/PnzBQUFQghHR0d/f39PT09j5wIAADCyOotdly5ddHyKI0eONFEYGfn5+fPn\nz09JScnLy6u2ysvLKzIyctq0aTY2NoYJAwAAYGrqLHaFhYV1rWrWrJmlpaVGozlx4oR+UtXi\n2rVrQUFBOTk5/v7+YWFh3t7etra2Qoh79+6dO3du165ds2bNSk1N3bFjR8uWLQ2WCgAAwHTU\nWezOnj1bz59t2LDhrbfeEkI8++yzTR+qNjNnzrx8+fKaNWtGjBhRc21lZeWSJUsmTZoUHx+/\ncOFCw0QCAAAwKfJnxVZz8eLFIUOGDB8+vKys7NtvvzXYLcW2bNkyZsyYWludEEKtVkdHR48c\nOXL9+vWGyQMAAGBqGlDsKioqPv3008DAwC1btowfP/7f//732LFj9Zesmlu3bvn6+tY/p2PH\njtevXzdMHgAAAFOj61mxGRkZUVFRR48e7dy581dfffXUU0/pNVZNHh4e2dnZ9c/Jysry8PAw\nTB4AAABTI7/FLj8//8033+zbt29OTs6nn356+PBhw7c6IcSwYcPWrl372WeflZaW1lxbVFQ0\ne/bsTZs2jRo1yvDZAAAATIHMFruVK1dOnTo1Ly/vhRdeSEhIMOLl4uLi4vbs2RMbGztnzpye\nPXt6enra2dlpNJrCwsKLFy9mZmbev38/ODh4xowZxkoIAABgXHUWu9OnT0+YMCE9Pd3b2/uH\nH34YMmSIIWPV5OTktH///sTExBUrVuzcubOyslK7ysLConv37hEREREREWq12oghAQAAjKjO\nYte5c+fS0tLg4OCpU6eq1eq0tLS6Zg4ePFg/2aqztLSMiYmJiYkpKSm5dOmSdOcJBwcHLy8v\nS0vLxj1naWnpypUrq9bEms6dO9e4JwcAADCkOouddCjbnj179uzZU/9TaDSaJg4lx9ra2t/f\nXwhRWVl54sSJQ4cOeXp6Nm438c2bN7/++uuKiop65tRzrWYAAADTUWexS0lJMWQOXWRkZKxa\ntWrRokXSw5UrV06bNk17fZMnnnjiyy+/fPrppxv0nG3atDlw4IDs6wYFBTUiMAAAgCHVWexe\neeUVQ+aQtXPnztDQUEtLy4SEBJVKtW7dujFjxtjZ2Y0YMcLZ2fnMmTO//vprSEjIvn37unfv\nbuywAAAARqDrdeyMLj4+3snJad++fSqVSgjx7rvvent779+/393dXZrw22+/Pfvss/Hx8T/8\n8INRkwIAABhHg28pZiy///772LFj/fz8hBB3797Nycl55513tK1OCNGrV69XXnlF9ohAAAAA\npTKbYldZWWljYyMtW1tbq1Sqtm3bVpvTtm3bkpISg0cDAAAwCWZT7Lp06bJ69er79+8LIays\nrJ566qn9+/dXnVBaWrp+/fqAgAAjBQQAADAysyl277///pkzZ4KDg3/++eeKioqEhITvvvtu\nxYoV9+/fLy8v/+2338LCwrKzs6Ojo42dFAAAwDjM5uSJwYMHf/3112+//XZoaKiNjU27du0s\nLS3HjRsXEREhhKisrFSpVO+888748eONnRQAAMA4zKbYCSEiIyOHDBmSkpLyyy+//Pvf/759\n+7aVlZWdnZ2Pj09QUNC4ceO6detm7IwAAABGY07FTgjh6uo6bdq0adOmGTsIAACAyZE/xm7A\ngAEDBgwoKiqqNl5UVCSt0k8wAAAANIz8Frtff/1VCFFeXl5tvLy8XFoFAAAAUyBf7I4ePSqE\ncHBwqDbu4OAgrQIAAIApkC92nTp1qnW8WbNmda0CAACA4ckfY/fUU08tWbLk7t27BkgDAACA\nRpMvdocOHYqKinJ3d//rX/+6ffv2Bw8eGCAWAAAAGkq+2OXm5i5ZsqRPnz5r1qwJCQnx8fGZ\nMWPG2bNnDRAOAAAAupMvdo899tgbb7zxyy+/XLt2bfHixX5+fgsWLPD39w8ODl62bFlBQYEB\nUgIAAEBWA+4V6+zsHBUVlZ6efvny5b///e8FBQWRkZFubm4TJkw4ffq0/iICAABAFw0odpLi\n4uJ9+/bt3btXKnOtW7detmxZp06d4uPjNRqNHhICAABAJw0odvv27Rs/frybm9uIESO2bt06\nfPjwHTt2XLx48dy5c0OHDo2Li4uPj9dfUAAAANRP/jp2ly5dWrFixbfffnvmzBkhRNeuXV9/\n/fWXX37ZyclJmuDp6bl27dqQkJDFixfHxcXpNS4AAADqIl/sfHx8Hjx44OjoGBUVFRkZ2b17\n95pzVCrVsGHDuMMYAACAEckXu6CgoNdff33kyJE2Njb1TAsNDU1NTW26YAAAAGgY+WK3e/du\nIcTx48ddXV1bt24tDR4/frysrKxr167aaX5+fn5+fnpKCQAAAFnyJ0+Ul5e//vrrnTp1Onbs\nmHZwx44d3bp1e+211yorK/UZDwAAALqSL3YJCQnffPNNeHi4t7e3dvD5558fNWpUcnLyokWL\n9BkPAAAAupIvdsnJyYMHD05LS2vXrp12MCAgYPXq1WFhYRQ7AAAAEyFf7M6ePfvss8/WuuqZ\nZ565ePFiU0cCAABAY8gXOwcHhwsXLtS66sKFC61atWriRAAAAGgU+WIXHh6+bNmyrVu3Vh0s\nLy//+uuv//GPf4SEhOgtGwAAABpA/nIn8+bN27ZtW3h4uJeXV0BAgJWV1Z07d06cOHH79m13\nd/d58+YZICUAAABkyW+xc3d3z8rKioqKKioq2r59e1pa2t69e9Vq9fjx4w8ePOjl5WWAlAAA\nAJAlv8VOCOHq6rp48eKkpKRr164VFxe7ubnZ2trqOxkAAAAaRKdiJ1GpVB4eHvqLAgAAgIch\nX+w0Gs26detWrFhx+fLl8vLymhOq3pECAAAAxiJf7D7//PPY2FghRIsWLSwsLPQfCQAAAI0h\nX+y++OKL0NDQpKSk9u3bGyAQAAAAGke+2F2/fn3dunW0OgAAABMnf7kTV1dXjUZjgCgAAAB4\nGPLFbvTo0SkpKQaIAgAAgIchvyt21qxZL7300ssvvzx27FgvL6+a50/4+fnpJxsAAAAaQL7Y\n2dvbSwurVq2qdQI7agEAAEyBfLEbPXq0paVl8+YNuJQxAAAADE++rtW1oQ4AAAAmRf7kCa2C\ngoLjx4/fuXNHf2kAAADQaDoVu127dvXo0cPBwaFTp04HDhyQBocOHfrrr7/qMxsAAAAaQL7Y\nZWZmhoSEnD59OjQ0VDt448aNgwcPhoWFHT58WJ/xAAAAoCv5Yjdnzhw3N7cTJ04kJydrB52d\nnbOzs93c3ObOnavHdAAAANCZfLE7cODAhAkT2rZtW23cxcUlKipq9+7d+gkGAACAhpEvdnfv\n3vX09Kx1lbu7e2FhYVNHAgAAQGPIFzs3N7eTJ0/Wumr37t0eHh5NHQkAAACNIV/swsLCkpKS\nfv/996qD+fn5H3744fLly8PDw/WWDQAAAA0gX+zi4+Pt7Ox69eoldbjp06d37drV3d39o48+\n8vLymjVrlv5DAgAAQJ5Ou2IPHTo0fvz4ixcvCiGOHDly5MgRe3v7CRMmHDx40NXVVf8hAQAA\nIE+nO8C6uLgkJSUlJibm5eUVFBTY29vT5wAAAEyNTsVOolKpXF1dqXQAAACmSb7YDRgwoJ61\nZWVlXMoOAADAFMgXu3puCGtvb29vb9+keQAAANBI8sWuvLy82khZWVlOTk5ycnJmZubmzZv1\nEwwAAAANI39WbPMaWrRo8fjjj3/66ad9+vR57733DJASAAAAsuSLXT1eeOGFH374oamiAAAA\n4GE04KzYmgoKCu7cudNUUQBAd2vWrFm3bp2xUxjISy+9NHLkSGOnAGAG5ItdrdWtvLz8+PHj\n7777brt27fSQCgBkbNq06VjG2uAAY+fQvz2nhIWFBcUOgC7ki13Lli3rWZuSktJ0YQCgAYID\nxJLXjR1C/95cJgqNnQGAuZAvdtItYquxsLBwd3f/n//5n/79++shFQAAABpMvtilpaUZIAcA\nAAAe0kOdFQsAAADTIb/FrkuXLlZWViqVSpenO3DgwENHAgAAQGPIF7vc3Nx79+4VFxdLD1Uq\nlUajkZZtbGzKysr0mA4AAAA6k98Ve/Lkye7du0+cOPH3338vLi5+8ODB3bt3d+3aNXz48ODg\n4Nu3b1dUYYDEAAAAqJV8sZs6daqfn9+iRYu6du1qbW0thHBwcHj66adTU1ObNWs2depU/YcE\nAACAPPlil5aWFhwcXOuqAQMGcEsxAAAAEyFf7O7du5ebm1vrqry8vLt37zZ1JAAAADSGfLEL\nDAxMTEz87bffqo3v27fvm2++6dChg36CAQAAoGHkz4qNi4sbPnx4796927Vr5+vra2NjU1xc\nfP78+fPnz6tUqq+++soAKQEAACBLvtgNHTr0119/XbBgwa5du3JycqRBS0vL5557bvr06QMG\nDNBzQgAAAOhEvtgJIfr169evX78HDx5cu3bt/v37NjY27u7uarVa3+EAAACgO52KnaSoqOjO\nnTtt2rRxcnLSXyAAAAA0jk73it21a1ePHj0cHBw6deqkvWmYtItWn9kAAADQAPLFLjMzMyQk\n5PTp06GhodrBGzduHDx4MCws7PDhw/qMBwAAAF3JF7s5c+a4ubmdOHEiOTlZO+js7Jydne3m\n5jZ37lw9pgMAAIDO5IvdgQMHJkyY0LZt22rjLi4uUVFRu3fv1k8wAAAANIx8sbt7966np2et\nq9zd3QsLC5s6EgAAABpDvti5ubmdPHmy1lW7d+/28PBo6kgAAABoDPliFxYWlpSU9Pvvv1cd\nzM/P//DDD5cvXx4eHq63bAAAAGgA+evYxcfHb9u2rVevXp07dxZCTJ8+ffr06SdPniwtLfXy\n8po1a5b+Q/4XjUaTk5Nz/vz5goICIYSjo6O/v39dO4sBAAAeHfLFzs3N7dChQ3FxcWvWrBFC\nHDlyRAjRunXriIiIuLg4FxcXvWf8j/z8/Pnz56ekpOTl5VVb5eXlFRkZOW3aNBsbG4PlAQAA\nMCk63XnCxcUlKSkpMTExLy+voKDA3t7e1dVV38mquXbtWlBQUE5Ojr+/f1hYmLe3t62trRDi\n3r17586d27Vr16xZs1JTU3fs2NGyZUsDZwMAADAF8sXuhx9+8PX1ffzxx1Uqlaurq+ErnWTm\nzJmXL19es2bNiBEjaq6trKxcsmTJpEmT4uPjFy5caPh4AAAARid/8sSoUaPS0tIMEKV+W7Zs\nGTNmTK2tTgihVqujo6NHjhy5fv16AwcDAAAwEfLFrm/fvrt27Xrw4IEB0tTj1q1bvr6+9c/p\n2LHj9evXDZMHAADA1Mjvil25cmVMTEx4ePjYsWP/9Kc/OTo6Vpvg5+enn2z/xcPDIzs7u/45\nWVlZXFcPAAA8snQ6K1Za+PHHH2udoNFomjJRHYYNG/bll18++eSTkydPtrKyqra2qKjok08+\n2bRp03vvvWeAMAAAACZIvtiNGjXK0tLSwsJCpVIZIFBd4uLi9uzZExsbO2fOnJ49e3p6etrZ\n2Wk0msLCwosXL2ZmZt6/fz84OHjGjBlGDAkAAGBE8sVu9erVBsghy8nJaf/+/YmJiStWrNi5\nc2dlZaV2lYWFRffu3SMiIiIiItRqtRFDAgAAGFGdxW7RokVdunTp27dv1cEjR444Ozu3adNG\n/8FqYWlpGRMTExMTU1JScunSJenOEw4ODl5eXpaWlo17zsrKyrS0tLKysnrmnDp1qnFPDgAA\nYEh1FrvJkydPmTKlWrHr2rXrxIkTFy1apP9g9bG2tvb399c+vHfv3qxZs1599dUOHTo09Kku\nX748YcKEkpKSeuZUVFQIQx1KCAAA0Gg63XnCxN27d+/jjz/u27dvI4qdt7f31atX65+TkZER\nFBRk3EMMAQAAZJlNsYuMjKxr1f3794UQCQkJGzduFEIsXbrUcLEAAABMhtkUu2XLltU/4eef\nf5YWKHYAAODRJH/nCRMRExOjVqu7dOny448/5v+348ePCyFWr14tPTR2UgAAAOMwm2L3t7/9\n7cCBA0KIQYMGffDBByqVyuk/HBwchBC2trbSQ2MnBQAAMA6zKXZCiB49ehw8eHDBggXJycmB\ngYGpqanGTgQAAGBC6jvG7sCBA3FxcdUGMzMzqw3WnKM/zZs3f++991566aU333zzpZdeGjJk\nSGJiIuerAgAAiPqL3W+//fbbb79VGzx48ODBgwerjhiy2El8fX1/+eWX5OTkqVOnBgYGTp06\n1cABAAAATFCdxS4lJcWQORrh1VdfDQsLmzJlSnx8vLGzAAAAGF+dxe6VV14xZI7GcXFx+ec/\n/zl27Nhff/3V19fX2HEAAACMyWyuY1ePQYMGDRo0yNgpAAAAjMyczooFAABAPSh2AAAACkGx\nAwAAUAiKHQAAgEJQ7AAAABSCYgcAAKAQFDsAAACFoNgBAAAoBMUOAABAISh2AAAACkGxAwAA\nUAiKHQAAgEJQ7AAAABSCYgcAAKAQFDsAAACFoNgBAAAoBMUOAABAISh2AAAACkGxAwAAUAiK\nHQAAgEJQ7AAAABSCYgcAAKAQFDsAAACFoNgBAAAoBMUOAABAISh2AAAACkGxAwAAUAiKHQAA\ngEJQ7AAAABSCYgcAAKAQFDsAAACFoNgBAAAoBMUOAABAISh2AAAACkGxTwxDRAAAHzJJREFU\nAwAAUAiKHQAAgEJQ7AAAABSCYgcAAKAQFDsAAACFoNgBAAAoBMUOAABAISh2AAAACkGxAwAA\nUAiKHQAAgEJQ7AAAABSCYgcAAKAQFDsAAACFoNgBAAAoBMUOAABAISh2AAAACkGxAwAAUAiK\nHQAAgEJQ7AAAABSCYgcAAKAQFDsAAACFoNgBAAAoBMUOAABAISh2AAAACkGxAwAAUAiKHQAA\ngEJQ7AAAABSCYgcAAKAQFDsAAACFoNgBAAAoBMUOAABAISh2AAAACkGxAwAAUAiKHQAAgEJQ\n7AAAABSCYgcAAKAQFDsAAACFoNgBAAAoBMUOAABAISh2AAAACkGxAwAAUIjmxg7QYBqNJicn\n5/z58wUFBUIIR0dHf39/T09PY+cCAAAwMnMqdvn5+fPnz09JScnLy6u2ysvLKzIyctq0aTY2\nNkbJBgAAYHRmU+yuXbsWFBSUk5Pj7+8fFhbm7e1ta2srhLh37965c+d27do1a9as1NTUHTt2\ntGzZ0thhAQAAjMBsit3MmTMvX768Zs2aESNG1FxbWVm5ZMmSSZMmxcfHL1y40PDxAAAAjM5s\nTp7YsmXLmDFjam11Qgi1Wh0dHT1y5Mj169cbOBgAAICJMJtid+vWLV9f3/rndOzY8fr164bJ\n8//au/eoqMvEj+PPVxhQHK6KGjDIRSrTzQTTPGKaWil6FMVVLNBFKdAkNTPTVUOtzN1Mk7Da\n4ymVjrdUVEqzMmU1NcqSTMkUkbxAbArIgMhtfn/MrzkTILeE78zD+/UXPt/L85nRwQ/fGwAA\nAJbGaoqdh4dHenp63ev88MMPHh4eLZMHAADA0lhNsQsNDf3444/ffPPN27dv11xaXFz8yiuv\n7NmzZ+LEiS2fDQAAwBJYzc0T8fHxR44cmTdv3rJly/r27avT6bRarcFg0Ov12dnZaWlpJSUl\nAwcOXLRokdpJAQAA1GE1xc7FxeX48eOJiYmbNm06fPhwZWWlaZFGowkKCpo6derUqVNtbGxU\nDAkAAKAiqyl2Qgg7O7s5c+bMmTOntLT08uXLxt884eTk5O3tbWdn17R9GgyGr7/+urS0tI51\nzpw507SdAwAAtCRrKnYmbdu2DQgIqDmen59fWFjo4+PT8F1lZWUNGTKkvLy83jUNBkPDdwsA\nANDyrObmCSHEjz/+OHLkSB8fn4EDB65bt878bKzRypUrfX19G7VPPz+/srIyQ52+/vprIYSi\nKHftlQAAADQDqzli9/XXXw8dOvT27dsODg7Xrl07evTo9u3bk5OT+QViAAAARlZzxG7FihVV\nVVXJycl6vb6oqOitt946duzYk08+WVxcrHY0AAAAi2A1xe7HH3+cOHFiaGiooij29vZz5sz5\n7LPP0tPTJ0yYUPOcLAAAQCtkNcUuNzfXz8/PfGTIkCHr16/ft2/fCy+8oFYqAAAAy2E119h1\n7tz51KlT1QYjIyMzMjJWrFjh5eU1b948VYIBAABYCKspduPGjUtISHjnnXdiYmI0Go1p/LXX\nXrt27dpLL7107do1zskCAIDWzGqK3ZIlS3bv3h0XF7dnz54vvvjCNK4oyocffujs7LxmzRoV\n4wEAAKjOaq6x69Chw8mTJ2fMmNGzZ89qixRFefvtt3fu3Onv769KNgAAAEtgNUfshBAdO3ZM\nTEy809Jx48aNGzeuJfMAAABYFKs5YgcAAIC6UewAAAAkQbEDAACQBMUOAABAEhQ7AAAASVDs\nAAAAJEGxAwAAkATFDgAAQBLW9IBiAEC9Fi9enJeXp3aKlqAoyqJFi7y8vNQOAlgQih0ASOX1\n118P6VXl4ap2jub34X9FSEgIxQ4wR7EDANnMGSGG9FA7RPPbelztBIDl4Ro7AAAASVDsAAAA\nJEGxAwAAkATFDgAAQBIUOwAAAElQ7AAAACRBsQMAAJAExQ4AAEASFDsAAABJUOwAAAAkQbED\nAACQBMUOAABAEhQ7AAAASVDsAAAAJEGxAwAAkATFDgAAQBIUOwAAAElQ7AAAACRBsQMAAJAE\nxQ4AAEASFDsAAABJUOwAAAAkQbEDAACQBMUOAABAEhQ7AAAASVDsAAAAJEGxAwAAkATFDgAA\nQBIUOwAAAElQ7AAAACRBsQMAAJAExQ4AAEASFDsAAABJUOwAAAAkQbEDAACQBMUOAABAErZq\nBwAAoNkVFRVt3769srJS7SAtwdPTc+TIkWqngDoodgAA+X3++ecxz0Z37ah2juZXUiaKytvr\n9Xq1g0AdFDsAgPyqqqrctCJztdo5mt/+dBGWUKV2CqiGa+wAAAAkQbEDAACQBMUOAABAEhQ7\nAAAASVDsAAAAJEGxAwAAkATFDgAAQBIUOwAAAElQ7AAAACRBsQMAAJAExQ4AAEASFDsAAABJ\nUOwAAAAkQbEDAACQBMUOAABAEhQ7AAAASVDsAAAAJEGxAwAAkATFDgAAQBIUOwAAAElQ7AAA\nACRBsQMAAJAExQ4AAEASFDsAAABJUOwAAAAkQbEDAACQBMUOAABAEhQ7AAAASVDsAAAAJEGx\nAwAAkATFDgAAQBK2agcAAACqqaysvHnzptopWoizs3ObNpIf0rK+YmcwGLKysi5evFhUVCSE\ncHZ2DggI0Ol0aucCAMD6xMbGrl+/Xu0ULWT27NmrV69WO0XzsqZil5+f/9prryUlJeXl5VVb\n5O3tHR0d/eKLL7Zr106VbAAAWKOCgoIJj4iXRqmdo/nF7xQFBQVqp2h2VlPscnJyBgwYkJWV\nFRAQEBIS0rVr1/bt2wshbt68mZmZmZqaumTJkp07dx46dMjV1VXtsAAAWI1OTiLIV+0Qza+j\no9oJWoTVFLvFixdfuXJl+/btf//732suraysfP/992fOnLl06dI1a9a0fDwAAADVWU2x+/TT\nTyMjI2ttdUIIGxubGTNm/Pe//921a1dji92vv/5aUVFRxwrXrl1r1A7rVVQqLlY/mSyhwhIh\nnBq3SW5Bq3hnKqsat35VVat4W3LyG73JzVut4p25eavRzy/IaR0fpSpD49avbB0fpdzGn2ws\nLGkV70xRqWgNx+yspthdv37d39+/7nW6d++enJzcqN1mZmYGBAQYDPV8e1AU5W7dR+Pk5LQh\nTexMuys7s3RTpzai2Tk5OU39z63mC2NRnJwa+s44OTkVlQr/Oc0ax1K4ujbuH8x7x8XW480X\nx4LE9mrcOxOxTv4LiYwa9VG6oW8tH6XOnRv3D+aDHSLpaPPFsSDPP9LI4w1WSKm301gIHx+f\nfv36bdu2rY51QkND09PTs7KyGrXnmzdvVlZW1r1OVVVVhw4dGrXbO6moqDDez9saaLVajUbT\nwJVLSkpu377drHkshKIoLi4uDV+/oKDAWj6nf5G9vb2Dg0MDV25VHyVHR0db24b+HN56Pkpt\n2rRxdnZu+Pp8lGpVXl6u1+ubNY/laNRHyUpZzcsLDQ1du3btww8/HBcXZ29vX21pcXHxv/71\nrz179syfP7+xe274D3x3ha2tLbd31MrBwaHh34lalUa1wNaDj9Kd8FG6Ez5KtdJoNHyUZGI1\nR+wKCgqGDh36/fffOzo69u3bV6fTabVag8Gg1+uzs7PT0tJKSkoGDhy4b98+rVardlgAAAAV\nWE2xE0KUlZUlJiZu2rTp9OnT5idPNRpNUFDQ1KlTp06damNjo2JCAAAAFVlTsTMpLS29fPmy\n8fIaJycnb29vOzs7tUMBAACozCqLHQAAAGqS/FfhAgAAtB4UOwAAAElQ7AAAACRBsQMAAJAE\nxQ4AAEASFDsAAABJUOwAAAAkQbEDAACQBMUOAABAEhQ7AAAASVDsAAAAJEGxAwAAkATFDgAA\nQBIUOwAAAElQ7AAAACRBsQMAAJAExQ4AAEASFDsAAABJUOwAAAAkQbEDAACQBMUOAABAEhQ7\nAAAASVDsAAAAJEGxAwAAkATFDgAAQBIUOwAAAElQ7AAAACRBsQMAAJAExQ4AAEASFDvJFRQU\nzJ4928fHx87OzsPDIzo6OicnR+1QFqG8vHzBggU2NjZ9+vRRO4sFyc/Pf/HFF7t27Wpvb+/r\n6xsaGnrixAm1Q1mEixcvPvvss/7+/vb29u7u7qGhoWlpaWqHsiwvvPCCoijR0dFqB1Hfhg0b\nlNq8+uqrakdT2c8//xwZGXnPPfdoNBp3d/exY8fyObrrbNUOgGZUVlY2dOjQ77//PiwsLDAw\nMDMzc9OmTV999dXJkyddXV3VTqemjIyMiIiI8+fPqx3Esty4cSMoKOjSpUsjR46cMmXKxYsX\nt23bduDAgbS0tL/97W9qp1PTuXPnBgwYUFRUNGHCBH9//wsXLmzfvn3fvn2pqan9+/dXO51F\n+O6779auXat2CktRUFAghJg0aZK3t7f5+IABA1RKZBHOnDnTv39/jUYzc+bMbt26ZWdnJyYm\nDhgw4MCBA0OGDFE7nUQMkNdbb70lhFi5cqVpZNu2bUKIuXPnqphKdYWFhe3atevTp8/58+ft\n7e2DgoLUTmQpnnvuOSFEQkKCaWTnzp1CiJCQEBVTWYLHH39cUZTU1FTTyK5du4QQEyZMUDGV\n5SgvL3/ooYd69eolhJg2bZracdT3yiuvCCG+/fZbtYNYlqeeekoI8dVXX5lG0tPThRCDBw9W\nMZV8OBUrs02bNjk6Os6aNcs0MmHChG7duiUlJRkMBhWDqauiomLGjBnHjh3r1q2b2lksi0aj\nGTp0aExMjGlk7Nix7dq1O3PmjIqpLEG/fv1efvnlRx991DQyevRojUZz7tw5FVNZjlWrVqWn\np7/xxhtqB7EUxiN2Li4uagexLJmZmUKI4OBg08iDDz7o5OR06dIl1TLJiGInrdLS0tOnT/ft\n29fe3t58PDg4OC8vLysrS61gqnNzc3vzzTc1Go3aQSzO6tWrv/zyS/N3pqysrKKiwsvLS8VU\nlmD58uWvv/66+Uhubm55ebmvr69akSxHZmbm0qVLY2NjH3nkEbWzWApTsausrLxy5crvv/+u\ndiKLcP/99wshzH8c+v333/V6fffu3dULJSGKnbQuX75cWVmp0+mqjXft2lUIcfHiRTVCwcq8\n//775eXl4eHhagexICUlJYcPHw4JCXF0dPznP/+pdhz1xcTEuLi4rFixQu0gFqSwsFAIsWbN\nGnd3d51O5+7uft99923evFntXCqbP3++q6trRETE0aNHc3Nzf/jhh/Dw8LZt2xrPXONu4eYJ\naRUVFQkh2rdvX21cq9WalgJ1SE1NnTdvXnBwcGxsrNpZLIWLi4vx/+yIiIjk5GQ/Pz+1E6ls\nw4YNBw8e3LFjh7Ozs/EwFcQfR+y2bNny0ksveXp6ZmRkJCYmPv3000VFReaXOrQ23bt3P378\n+Lhx4wYOHGgc8fb2/vLLL/v166duMMlQ7CSnKEq1EePVdTXHAXNbtmyJiorq2bPnnj17bG35\nRvH/pk+ffuPGjZ9++mnz5s2XLl3auHFja+52eXl5c+fOHTVqVFhYmNpZLMvixYtnzpw5fPhw\n04/WERERgYGBCxcujIqKsrOzUzeeWjIyMkaOHFlRUbFq1ap77703Ly/vrbfeGjFixI4dO4YN\nG6Z2Onnw/VpaTk5OorYjczdv3hRCODo6qpAJ1sBgMMTHxy9btmz48OHbt2/nn4o50wnHw4cP\njxo1auzYsT/88EObNq30mpZZs2aVlZUlJiaqHcTi1Hx4xwMPPBASEpKcnJyenv7www+rkkp1\nU6dO/e2333755RdPT0/jSHh4+L333vuPf/wjKyuL657vllb6/ag18Pb2trW1zc7OrjZuvC8p\nICBAjVCwdAaDITo6etmyZXFxcZ988gmt7k4GDx48ZsyYH3/8sdXeGLt///6tW7fOmTOnTZs2\nV65cuXLlyrVr14QQJSUlV65cMf4ACXOdOnUSQuj1erWDqEOv13/zzTf9+vUztTohhIODw9Ch\nQ69evfrLL7+omE0yFDtp2dnZBQUFpaWllZSUmAarqqpSU1N1Ol21x2YCRnPmzPnggw9ef/31\ntWvX2tjYqB3HIly9erVXr16TJ0+uNl5aWiqEKC4uViOU+g4ePCiEWL58ue4PPXr0EEJs2bJF\np9NVu4m4VdHr9e++++6WLVuqjRsfG2S8fa0VunXrlsFgMH5qzBlHao6j6VR9ih6a13/+8x8h\nRHx8vGnk3XffFUIsXbpUxVQWhQcUmzM+jnjWrFlqB7E4Xl5ednZ2J06cMI2cO3dOq9VqtVrj\nf1et0NmzZ1P+bOvWrUKIJ554IiUlJSMjQ+2AqqmsrPT09NRqteZvwu7du4UQvXv3VjGY6nx9\nfY1PfzSN5Ofnu7m5OTk5lZaWqhhMMoqhFT+oVnqVlZWPPfbYkSNHxowZExgYmJGRsW3btp49\ne544ccLBwUHtdKpJTU3dv3+/8es333zT3d19ypQpxj/OmzevQ4cO6kVTWbdu3TIzM+Pi4mr+\n8zA+p0CVVJZg9+7d48ePb9OmTVhYmL+//9WrVz/++OPi4uJ33nnH+Os6IIQoKChwdXWdNm3a\n+vXr1c6isr1794aGhjo4OISHh3t4ePz000+7d+92dHQ8dOhQYGCg2ulUk5ycPH78eFdX19jY\nWH9//5ycnPXr12dlZSUmJs6YMUPtdBJRu1mieRUVFRl/p7tGo/H09HzuueeuX7+udiiV1fHA\nrfPnz6udTk11fKPIyspSO53KTpw4ERoa6u7ubmNj4+LiMmzYsL1796odyrLk5+cLfqXYH44d\nOzZixAgXFxdbW1sPD4/Jkye38m8vRseOHTN+jmxtbV1dXYcNG/bpp5+qHUo2HLEDAACQBDdP\nAAAASIJiBwAAIAmKHQAAgCQodgAAAJKg2AEAAEiCYgcAACAJih0AAIAkKHYAAACSoNgBAABI\ngmIHAAAgCYodAACAJCh2AAAAkqDYAQAASIJiBwAAIAmKHQAAgCQodgAAAJKg2AEAAEiCYgcA\nACAJih0AAIAkKHYAAACSoNgBAABIgmIHAAAgCYodAACAJCh2AAAAkqDYAQAASIJiBwAAIAmK\nHQAAgCQodgAAAJKg2AEAAEiCYgcAACAJih0AAIAkKHYAAACSoNgBsCybN2/28vKytbWdN2+e\n2ln+X3h4uKIoV65cqWNpbm5uC6dSV3R0tKIoFy5cUDsIgD+h2AH4k48++khRlLZt2/7yyy81\nl3br1u2hhx5qvtkLCwujo6P1ev3y5cuffPLJOtY0GAy7du0KCwvT6XRt27Zt27atr6/vtGnT\nvvnmm+aLV6uHHnroySeftLe3b7EZjX9HiqIcOHCg1hVmz55tXKGioqLhu33jjTcoaoC1o9gB\nqMXt27dnzJjR8vOeP3/+1q1bTz/99IIFC4YNG3an1fLz8x9//PGwsLDPP/+8d+/ezzzzzKRJ\nkzp16vTBBx/079//jTfeaMnML7/88meffebq6tqSkwohbGxsNmzYUHO8oqJiy5YtNjY2jdpb\nTk7OggULKHaAtaPYAajFwIEDDx48+NFHH7XwvKWlpUIIR0fHOtYxGAwTJ048ePDgxIkTs7Oz\n9+7dm5CQ8OGHH37zzTdpaWleXl4LFixISUlpqciq6du37+7duwsKCqqNf/bZZ3l5eQ8++GCj\n9vbtt9/evWgAVEOxA1CLuXPndu3ade7cufn5+XWslp2dHRUV5enpaWdn17Fjx9GjR6elpdW9\n5zo2GT58+MCBA4UQK1euVBQlNja21j188sknX3zxRf/+/Tdv3uzm5ma+6OGHH961a9czzzzj\n5OTUqJxpaWljx47t2LGjnZ2dj49PZGTkpUuX6n4hJubX2D311FOKouj1+vnz5/v4+Njb2+t0\nutWrVxsMhjr20LTZhw8fXlpaum3btmrjGzdu9Pb2vu+++xo+y6hRo8aMGSOEGDFihKIoR48e\nNY7n5uZGR0d7enq2b9++V69eb7/9drVzu23atFm5cqWfn5+9vb23t/fy5cvrfqUAmput2gEA\nWCKNRpOQkDB69OiXX375/fffr3Wdy5cv9+3bt6SkZPr06T169Lh69eq6deseffTRL7/8Mjg4\nuAmbvPLKK4MGDVq4cOG4ceMiIyN9fX1r3UlSUpIQYsGCBW3a1PKjaZ8+ffr06dOonCdPnhw0\naJCbm9usWbO6dOly8eLFxMTEzz///OzZsx06dGjU+2ZnZyeEGD9+vK+v79atW6uqqpYuXfrC\nCy+4uLhERUXVukmTZw8MDPT19d2wYUNMTIxpsKCgICUl5fnnn//1118bPsuiRYvc3NySkpKW\nLFnSu3fvBx54QAjxv//9r0+fPnq9fvLkyV27dj18+PDs2bNPnz69fv16025fffXVU6dOPfvs\nszY2NgkJCUuWLOnWrdukSZMa9aYBuJsMAGDGWJtSUlIMBsOYMWMURTl27Jhpqb+/f69evYxf\nT5kyRQixa9cu09KzZ8/a2Ng88sgjd9p5vZscOXJECDF//vw6Enp7eyuKUlRU1MBXVO+k69at\nCwwMPHTokGmFhIQEIURCQoLxjxMnThRCXL58udb9G5fm5OQYDIZp06YJISZNmmRampmZKYQY\nNWrUneLVO3tNpr+j+Ph4IcTPP/9sWvTee+8JIU6fPm1MVV5e3sBZVqxYIYTYv3+/aYXp06cL\nIQ4cOGAaGTlypBDip59+Mr3S4ODgsrIy49KTJ08KIUaPHn2n2ABaAKdiAdxRQkKCg4NDTExM\nzZsrDQbD7t27O3fuHBoaahrs3r17//79T5w4cf369Zp7a8Imtfrtt9+cnZ21Wm1DVm7IpNOn\nTz958uTgwYOFEOXl5aWlpcZDVg0/G1uNsUoa+fn5OTg43OlRKX9x9ilTpiiKYn4LxcaNG4OC\ngnr27PkXZzEYDNu3b9fpdI8//rhpcO3atV999VXnzp1NI3PnztVoNMave/fubWNjc+3atXpj\nA2g+FDsAd6TT6eLj40+fPr169epqi3JzcwsLC3v06KEoivm48dKuWh+V0oRNamVjY1NVVVVt\nMDg4WPkz410FDZw0KSlp0KBBrq6udnZ27dq1Gzp0qBCiUc8KMeft7W3+R41GU15eXsf6TZ7d\nx8dn0KBBSUlJxjfkwoULx48fN6+VTZ4lJyfn+vXr999/v/n75ufn99hjj3Xs2NE0EhAQYPpa\nURStVnvr1q16YwNoPhQ7AHWZPXv2gw8+GB8fn52dbT5eXFwshGjfvn219du1a2daWk0TNqmV\nh4fHzZs3b9y4YT4YGhoa8wd/f/9GTbpw4cLJkyeXlJSsXr368OHDx48fN7+MrAlMB7Ea4i/O\nHhUVdfXq1S+++EIIsWnTJo1GU+slbo2dxdjP6n04X0s+vQ9AQ3DzBIC62NravvfeewMGDIiL\ni9u7d6/pfgXjmdCabcw4UuvzSpqwSa2Cg4MvXLiwb9++iIgI0+CLL75o+jo8PNx4ZVtDJi0t\nLV2zZo1Opzt06JDp9G5hYWEDw/xFf3328ePHz5w5c8OGDU888URSUtLIkSPNj6g1eZYuXboI\nIWo+SwWAheOIHYB69O/fPzo6OiUlJTk52XjXpxCiS5cubm5uGRkZhj8/3uLs2bOKotR81kbT\nNqlVdHS0EGLZsmUlJSX1rlzvpLm5ubdu3erTp4/5RXupqakNDPMX/fXZHRwcxo8fn5KScvTo\n0UuXLtV6HrYJs7Rv397d3T0jI8P8JPK5c+feeeedM2fONDwegBZGsQNQv5UrV7q7uz///PPm\n9WjcuHE5OTl79uwxjZw6dSotLW3IkCEuLi617qcJm9Q0YMCAyMjI8+fPDx8+vNoJ4tu3b69d\nuzYlJcXR0dHUQeuetHPnzoqimN9DcOrUqU2bNok/npbcrO7K7FFRUcXFxYsWLerYsaPxxtUm\nzGL8TRXmV8iNGTPm+vXrGzduNI3Ex8fHxcXdvn27oS8PQMtT74ZcAJbI/HEn5kz/wZsed3L1\n6tUuXbpotdqFCxdu3Lhx6dKlnTp1cnR0TE9Pv9PO692kIY87MRgMpaWlxvOwdnZ2w4YNmzFj\nRmxsbEhIiPF8bmBg4OnTpxs+6ahRo4QQMTExW7ZsWbx4saur6759+2xtbb28vDZv3qzX6xv7\nuJPz58+br+Ds7NyjR487vZZ6Z6+5Sc2/I+NlhXFxcdVSmR53Uu8sO3bsEEL07dt31apVaWlp\nBoPh8uXLXbp0sbW1nTlz5r///W/jHiZPnmzcYRNeKYAWQLED8Cd3KnYGg+Gxxx4zL3YGg+HX\nX3+Nioq65557bG1tO3XqFB4efvbs2br3X/cmDSx2RocOHYqIiPDx8Wnbtq1Wqw0ICIiMjNyz\nZ09VVVWjJs3Ly3vqqafc3d2dnZ2HDBly5MgRg8GwdOlSrVbbpUuXnJycZi129c5ec5Oaf0fL\nly8XQnz33XfVUpmKXb2zlJWVhYWFtWvXztXV9eOPPzZudenSpYiIiE6dOmk0Gj8/v1WrVlVU\nVBgXUewAy6QY+PUvAAAAUuAaOwAAAElQ7AAAACRBsQMAAJAExQ4AAEASFDsAAABJUOwAAAAk\nQbEDAACQBMUOAABAEhQ7AAAASVDsAAAAJEGxAwAAkATFDgAAQBIUOwAAAElQ7AAAACRBsQMA\nAJAExQ4AAEASFDsAAABJUOwAAAAkQbEDAACQBMUOAABAEhQ7AAAASVDsAAAAJEGxAwAAkATF\nDgAAQBIUOwAAAElQ7AAAACRBsQMAAJAExQ4AAEASFDsAAABJUOwAAAAk8X/u8LLoVQMriQAA\nAABJRU5ErkJggg==",
            "text/plain": [
              "Plot with title “Goal by Home Team in Full Time”"
            ]
          },
          "metadata": {
            "tags": [],
            "image/png": {
              "width": 420,
              "height": 420
            }
          }
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 34
        },
        "id": "aJnRpvzE1TOT",
        "outputId": "da147171-e97b-47d0-c9cf-55308a7b2b8a"
      },
      "source": [
        "c(round(mean(data$FTHG),3), round(var(data$FTHG),3))"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "[1] 1.516 1.559"
            ],
            "text/latex": "\\begin{enumerate*}\n\\item 1.516\n\\item 1.559\n\\end{enumerate*}\n",
            "text/markdown": "1. 1.516\n2. 1.559\n\n\n",
            "text/html": [
              "<style>\n",
              ".list-inline {list-style: none; margin:0; padding: 0}\n",
              ".list-inline>li {display: inline-block}\n",
              ".list-inline>li:not(:last-child)::after {content: \"\\00b7\"; padding: 0 .5ex}\n",
              "</style>\n",
              "<ol class=list-inline><li>1.516</li><li>1.559</li></ol>\n"
            ]
          },
          "metadata": {
            "tags": []
          }
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "umZa-ETH2CX9"
      },
      "source": [
        "* Number of goals scored by home team in English premier league (EPL) may follow Poisson distribution\n",
        "\n",
        "* Suppose $X$: Number of goals scored by home team in EPL\n",
        "\n",
        "* $D=\\{X_1,X_2,\\cdots,X_n\\} \\sim_{iid} Poisson(\\lambda)$\n",
        "\n",
        "* $T=\\sum_{i=1}^n X_i \\sim Poisson(n\\lambda)$ is a **sufficient statistics**\n",
        "\n",
        "* The concept of **sufficient statistics** is accepted and important in both classical frequentist and Bayesian statistics\n",
        "\n",
        "* One can show:\n",
        "Fisher's Information in $D$ is equal to Fisher's information in $T$\n",
        "\n",
        "* Under the exact same prior $p(\\lambda)$, one can show:\n",
        "$$\n",
        "p(\\lambda | x_1,x_2,\\cdots,x_n ) = p(\\lambda|t) =\\frac{f(t|\\lambda)p(\\lambda)}{\\int f(t|\\lambda)p(\\lambda)d\\lambda}\n",
        "$$\n",
        "\n",
        "* Consider conjugate prior for Poisson: $\\lambda \\sim Gamma(a,b)$ where $\\mathbb{E}(\\lambda)=\\frac{a}{b}$\n",
        "\n",
        "\\begin{eqnarray*}\n",
        "p(\\lambda|t)&\\propto& f(t|\\lambda)p(\\lambda)\\\\\n",
        "&=& e^{-(n\\lambda) }\\frac{(n\\lambda)^t}{t!}\\times \\frac{b^a}{\\Gamma(a)}e^{-b\\lambda}\\lambda^{a-1}\\\\\n",
        "&=& e^{-(n+b)\\lambda}\\lambda^{t+a-1}\\times \\frac{n^tb^a}{t!\\Gamma(a)}\\\\\n",
        "&\\propto& e^{-(n+b)\\lambda}\\lambda^{t+a-1}\n",
        "\\end{eqnarray*}\n",
        "\n",
        "* The **posterior distribution** of $\\lambda$ is $Gamma(a+t,b+n)$\n",
        "\n",
        "* The posterior mean is:\n",
        "\\begin{eqnarray*}\n",
        "\\mathbb{E}(\\lambda|t)&=&\\frac{a+t}{b+n}\\\\\n",
        "&=&\\frac{b}{b+n}\\frac{a}{b}+\\frac{n}{b+n}\\frac{t}{n}\\\\\n",
        "&=&\\frac{b}{b+n}\\mathbb{E}(\\lambda)+\\frac{n}{b+n}\\bar{X},\n",
        "\\end{eqnarray*}\n",
        "\n",
        "* Posterior mean is weighted average of prior mean and MLE, where $\\bar{X}=\\frac{t}{n}=\\frac{\\sum_{i=1}^nx_i}{n}$ is the MLE of $\\lambda$\n",
        "\n",
        "* Suppose we choose $Gamma(a=10,b=10)$ as prior distribution"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 34
        },
        "id": "tU_rgRNC2A3f",
        "outputId": "ac58e994-a84b-4aab-eeec-cc540cd76169"
      },
      "source": [
        "t = sum(data$FTHG)\n",
        "t"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "[1] 576"
            ],
            "text/latex": "576",
            "text/markdown": "576",
            "text/html": [
              "576"
            ]
          },
          "metadata": {
            "tags": []
          }
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "V7W2P6AUMRzk"
      },
      "source": [
        "* The posterior is then $Gamma(a+t=10+576,b+n=10+380)$"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "d2M287FKMd16",
        "outputId": "06d5f160-6cb8-4e3f-eba5-00cd5417a0e4"
      },
      "source": [
        "t = sum(data$FTHG)\n",
        "n = nrow(data)\n",
        "a = 10\n",
        "b = 10\n",
        "## MLE\n",
        "lambda_hat = t/n\n",
        "cat('MLE of lambda = ',lambda_hat,'\\n')\n",
        "se_mle = sqrt(lambda_hat/n)\n",
        "\n",
        "\n",
        "## Prior mean of lambda\n",
        "cat('Prior mean of lambda = ',a/b,'\\n')\n",
        "\n",
        "\n",
        "## Posterior mean\n",
        "cat('Posterior mean of lambda = ',(t+a)/(n+b),'\\n')\n",
        "\n"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "MLE of lambda =  1.515789 \n",
            "Prior mean of lambda =  1 \n",
            "Posterior mean of lambda =  1.502564 \n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 67
        },
        "id": "Mu3YKQirQfUU",
        "outputId": "6841915a-b6bc-4266-8d00-3c5e997c7d9c"
      },
      "source": [
        "## 95% Bayesian credible interval\n",
        "Bayesian_CI = qgamma(c(0.025,0.5,0.975),shape = t+a,rate=n+b)\n",
        "cat('Bayesian 95% CI = ', Bayesian_CI,'\\n')\n",
        "\n",
        "## 95% Frequentist Interval\n",
        "Frequentist_CI = c(lambda_hat-1.96*se_mle,lambda_hat+1.96*se_mle)\n",
        "cat('Frequentist 95% CI = ', Frequentist_CI,'\\n')\n",
        "\n",
        "## 95% prior Bayesian confidence interval\n",
        "qgamma(c(0.025,0.5,0.975),shape = a,rate=b)"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "Bayesian 95% CI =  1.383355 1.501709 1.62663 \n",
            "Frequentist 95% CI =  1.392 1.639579 \n"
          ],
          "name": "stdout"
        },
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "[1] 0.4795389 0.9668715 1.7084803"
            ],
            "text/latex": "\\begin{enumerate*}\n\\item 0.479538869613243\n\\item 0.966871461471413\n\\item 1.70848034514192\n\\end{enumerate*}\n",
            "text/markdown": "1. 0.479538869613243\n2. 0.966871461471413\n3. 1.70848034514192\n\n\n",
            "text/html": [
              "<style>\n",
              ".list-inline {list-style: none; margin:0; padding: 0}\n",
              ".list-inline>li {display: inline-block}\n",
              ".list-inline>li:not(:last-child)::after {content: \"\\00b7\"; padding: 0 .5ex}\n",
              "</style>\n",
              "<ol class=list-inline><li>0.479538869613243</li><li>0.966871461471413</li><li>1.70848034514192</li></ol>\n"
            ]
          },
          "metadata": {
            "tags": []
          }
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "oEis9UxSVwI8"
      },
      "source": [
        "Sampling distribution\n",
        "$$\n",
        "\\hat{\\lambda}\\approx N(\\lambda,\\frac{\\lambda}{n})\n",
        "$$\n",
        "Frequentist 95% CI\n",
        "$$\n",
        "\\hat{\\lambda} \\pm 1.96 \\sqrt{\\frac{\\hat{\\lambda}}{n}}\n",
        "$$"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "yOIlFuYCSI3f"
      },
      "source": [
        "**Criticism of Bayesian Statistics** How to choose the hyperparameters? In this case a and b of the Gamma distribution"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "_o3a_uMpSXi2",
        "outputId": "2b364712-4393-4761-bbc3-e77f716aea5a"
      },
      "source": [
        "a=1\n",
        "b=1\n",
        "## 95% Bayesian credible interval\n",
        "Bayesian_CI = qgamma(c(0.025,0.5,0.975),shape = t+a,rate=n+b)\n",
        "cat('Bayesian 95% CI = ', Bayesian_CI,'\\n')"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "Bayesian 95% CI =  1.393371 1.513561 1.640472 \n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "RRPzpPoBSySO",
        "outputId": "ea82db1c-6aba-424f-a9b1-18e68a505018"
      },
      "source": [
        "a=10\n",
        "b=10\n",
        "## 95% Bayesian credible interval\n",
        "Bayesian_CI = qgamma(c(0.025,0.5,0.975),shape = t+a,rate=n+b)\n",
        "cat('Bayesian 95% CI = ', Bayesian_CI,'\\n')"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "Bayesian 95% CI =  1.383355 1.501709 1.62663 \n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "eLx8b5Q-S0x-",
        "outputId": "19ebd157-93f2-4e69-81a3-c9b3f3d954f9"
      },
      "source": [
        "a=100\n",
        "b=100\n",
        "## 95% Bayesian credible interval\n",
        "Bayesian_CI = qgamma(c(0.025,0.5,0.975),shape = t+a,rate=n+b)\n",
        "cat('Bayesian 95% CI = ', Bayesian_CI,'\\n')\n"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "Bayesian 95% CI =  1.304155 1.407639 1.516457 \n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "4G1wwoHV2n06",
        "outputId": "8869b835-c223-4afe-bb25-c60507a6b62c"
      },
      "source": [
        "a=1000\n",
        "b=1000\n",
        "\n",
        "## 95% Posterior Bayesian credible interval\n",
        "Bayesian_CI = qgamma(c(0.025,0.5,0.975),shape = a,rate=b)\n",
        "cat('Prior Bayesian 95% CI = ', Bayesian_CI,'\\n')\n",
        "\n",
        "\n",
        "## 95% Posterior Bayesian credible interval\n",
        "Bayesian_CI = qgamma(c(0.025,0.5,0.975),shape = t+a,rate=n+b)\n",
        "cat('Posterior Bayesian 95% CI = ', Bayesian_CI,'\\n')\n"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "Prior Bayesian 95% CI =  0.938973 0.9996667 1.062921 \n",
            "Posterior Bayesian 95% CI =  1.086336 1.141787 1.199095 \n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "8T2VDaQqbYJZ"
      },
      "source": [
        "$x \\sim f(x)$ $y=g(x)$ $x=g'(y)$ $J(x,y)=\\frac{\\partial x}{\\partial y}$\n",
        "\n",
        "$\n",
        "h(y)=f(g'(y))\\frac{\\partial x}{\\partial y}\n",
        "$\n",
        "\n",
        "$\n",
        "p(\\theta)\\propto 1 ~~if ~~\\theta \\in \\mathbb{R}\n",
        "$\n",
        "\n",
        "$\\lambda \\in R^+$ \n",
        "\n",
        "$\\theta = \\log(\\lambda)$ $J =\\frac{1}{\\lambda}$\n",
        "\n",
        "$p(\\lambda) \\propto \\frac{1}{\\lambda}$\n",
        "\n",
        "$Gamma(a=1,b=0)$\n",
        "\n",
        "$Gamma(a+t,b+n)$"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "JoaZRwVoUzo9"
      },
      "source": [
        "One approach to avoid this criticism is **Jeffrey's Prior**\n",
        "\n",
        "* The Jeffreys prior for the rate parameter $\\lambda \\geq 0$ is\n",
        "\\begin{eqnarray*}\n",
        "p(\\lambda) &\\propto& \\frac{1}{\\sqrt{\\lambda}}\\\\\n",
        "&\\equiv& Gamma(a=\\frac{1}{2},b=0)\n",
        "\\end{eqnarray*}\n",
        "\n",
        "* So **Jeffrey's Prior** for Poisson is an improper prior.\n",
        "\n",
        "* What about posterior distrinution under Jeffrey's prior?\n",
        "\\begin{eqnarray*}\n",
        "p(\\lambda|t,n) &\\propto&  Gamma(a+t=\\frac{1}{2}+t~,~b+n=0+n)\n",
        "\\end{eqnarray*}\n",
        "\n",
        "So if $n>0$, which is trivially true, the __posterior distribution is proper__.\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "dMDfLC1eXEwJ",
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "outputId": "399b91f1-2f08-4731-f3d8-a7e4fd9dc886"
      },
      "source": [
        "a=0.5\n",
        "b=0\n",
        "\n",
        "## Posterior mean\n",
        "cat('posterior mean = ',(a+t)/(b+n),'\\n')\n",
        "\n",
        "## 95% Bayesian credible interval\n",
        "Bayesian_CI = qgamma(c(0.025,0.5,0.975),shape = t+a,rate=n+b)\n",
        "cat('Bayesian 95% CI = ', Bayesian_CI,'\\n')"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "posterior mean =  1.517105 \n",
            "Bayesian 95% CI =  1.395775 1.516228 1.64342 \n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Od7c4UfqVZhQ"
      },
      "source": [
        "$$\n",
        "\\mathbb{P}(a\\leq \\lambda \\leq b | D) =0.95\n",
        "$$"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "mJ3DTv1aauYj"
      },
      "source": [
        "**Alternate approach**: \n",
        "\n",
        "-  Choose Cauchy distribution for location parameter  (Berger, 1985)\n",
        "-  Choose Half-Cauchy distribution for scale parameter (Gelman, 2006)\n",
        "\n",
        "__Why__?\n",
        "\n",
        "* Mean, variance and other higher order moments does not exists for neither Cauchy nor Half-Cauchy distribution\n",
        "\n",
        "* The pdf of Half-cauchy distribution:\n",
        "\\begin{eqnarray*}\n",
        "p(\\theta)=\\frac{2}{\\pi(1+\\theta^2)}~~~\\theta>0\n",
        "\\end{eqnarray*}\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "okvpf0eKhaoD"
      },
      "source": [
        "dHalfCauchy= function(theta,log=FALSE){\n",
        "  if(log==FALSE){\n",
        "    return(2/(pi*(1+theta^2)))\n",
        "  }else{\n",
        "    return(log(2)-log(pi)-log(1+theta^2))\n",
        "  }\n",
        "}\n",
        "\n",
        "## Test\n",
        "# dHalfCauchy(theta=1,log=FALSE)\n",
        "# dHalfCauchy(theta=1,log=TRUE)\n",
        "# log(dHalfCauchy(theta=1,log=FALSE))"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "z6GhtaIJjFre"
      },
      "source": [
        "As $T\\sim Poisson(n\\lambda)$ and $\\lambda\\sim \\text{Half-Cauchy}(0,1)$ the posterior will be proper, i.e.,\n",
        "\n",
        "\\begin{eqnarray*}\n",
        "p(\\lambda | t,n) &\\propto& l(\\lambda|t,n)p(\\lambda)\\\\\n",
        "&\\propto& e^{-n\\lambda}{n\\lambda}^{t}\\times \\frac{2}{\\pi(1+\\lambda^2)}\n",
        "\\end{eqnarray*}\n",
        "\n",
        "We will try to estimate the **posterior mode**\n",
        "$$\n",
        "\\max_{\\lambda} p(\\lambda | t,n)\n",
        "$$"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "bgGZR_5oiY3Q",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 34
        },
        "outputId": "960a14b8-fe83-443d-d725-895d353e567b"
      },
      "source": [
        "logLike = function(lambda,t,n){\n",
        "  l = dpois(t,lambda = n*lambda,log=TRUE)\n",
        "  return(l)\n",
        "}\n",
        "\n",
        "logPosterior = function(lambda,t,n){\n",
        "  lpost = logLike(lambda=lambda,t=t,n=n)+dHalfCauchy(theta=lambda,log=TRUE)\n",
        "  return(lpost)\n",
        "}\n",
        "\n",
        "lambda.init = 1\n",
        "logPosterior(lambda=lambda.init,t=t,n=n)"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "[1] -48.82124"
            ],
            "text/latex": "-48.8212379191782",
            "text/markdown": "-48.8212379191782",
            "text/html": [
              "-48.8212379191782"
            ]
          },
          "metadata": {
            "tags": []
          }
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "eWbdB6w-mMpS"
      },
      "source": [
        "one-dimensional optimisation"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "NloJxi4Mlsxt",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 112
        },
        "outputId": "2b4bffb5-2b04-4e56-b287-826a0484bf2d"
      },
      "source": [
        "optimize(f=logPosterior,lower = 0,upper=10,t=t,n=n,maximum=TRUE)"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "$maximum\n",
              "[1] 1.512117\n",
              "\n",
              "$objective\n",
              "[1] -5.740236\n"
            ],
            "text/latex": "\\begin{description}\n\\item[\\$maximum] 1.51211696386756\n\\item[\\$objective] -5.74023555223261\n\\end{description}\n",
            "text/markdown": "$maximum\n:   1.51211696386756\n$objective\n:   -5.74023555223261\n\n\n",
            "text/html": [
              "<dl>\n",
              "\t<dt>$maximum</dt>\n",
              "\t\t<dd>1.51211696386756</dd>\n",
              "\t<dt>$objective</dt>\n",
              "\t\t<dd>-5.74023555223261</dd>\n",
              "</dl>\n"
            ]
          },
          "metadata": {
            "tags": []
          }
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "dWBOSqAHmaaE"
      },
      "source": [
        "**Question:** How to estimate posterior mean? Bayesian interval?"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "zvcwAAPxmtZU"
      },
      "source": [
        "**Question:**  What is posterior mean?\n",
        "\n",
        "$$\n",
        "\\mathbb{E}(\\lambda|D)=\\int \\lambda p(\\lambda|D)d\\lambda\n",
        "$$\n",
        "\n",
        "We don't know $p(\\lambda|D)$ fully. We know upto its kernel."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "MhxZQIv9nRg-"
      },
      "source": [
        "**Question:** What is Bayesian interval?\n",
        "\n",
        "$$\n",
        "\\mathbb{P}(a \\leq \\lambda \\leq b | D) = 0.95\n",
        "$$\n",
        "\n",
        "$$\n",
        "\\int_a^b p(\\lambda |D) d\\lambda  = 0.95\n",
        "$$"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "_EkCbflSoBYR"
      },
      "source": [
        "* So finding posterior mode is an optimisation problem\n",
        "\n",
        "* Estimation of posterior mean, posterior variance or 95% CI is an integration problem"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "bwTYo9AggAu8"
      },
      "source": [
        "**Reference**\n",
        "\n",
        "1.   A. Gelman (2006), \"Prior distributions for variance parameters in hierarchical models\" *Bayesian Analysis*, Vol. 1, N. 3, pp. 515–533\n",
        "http://www.stat.columbia.edu/~gelman/research/published/taumain.pdf\n",
        "\n",
        "2.   N. G. Polson and J. G. Scott (2012), \"On the Half-Cauchy Prior for a Global Scale Parameter\", *Bayesian Analysis*, Vol. 7, No. 4, pp. 887-902 https://projecteuclid.org/euclid.ba/1354024466\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 470
        },
        "id": "czd485Ut5JxE",
        "outputId": "2dd7246f-a6ab-4d63-a71d-05ceae79799f"
      },
      "source": [
        "n1 = 300\n",
        "length(data$FTHG[1:n1])\n",
        "t = sum(data$FTHG[1:n1])\n",
        "t\n",
        "## lambda ~ Gamma(a+t,b+n)\n",
        "a = 1\n",
        "b = 1\n",
        "lambda = seq(1,2,length.out = 100)\n",
        "plot(lambda,dgamma(lambda,a + t, b+n1),lwd=2,type='l')"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "[1] 300"
            ],
            "text/latex": "300",
            "text/markdown": "300",
            "text/html": [
              "300"
            ]
          },
          "metadata": {
            "tags": []
          }
        },
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "[1] 450"
            ],
            "text/latex": "450",
            "text/markdown": "450",
            "text/html": [
              "450"
            ]
          },
          "metadata": {
            "tags": []
          }
        },
        {
          "output_type": "display_data",
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAA0gAAANICAMAAADKOT/pAAADAFBMVEUAAAABAQECAgIDAwME\nBAQFBQUGBgYHBwcICAgJCQkKCgoLCwsMDAwNDQ0ODg4PDw8QEBARERESEhITExMUFBQVFRUW\nFhYXFxcYGBgZGRkaGhobGxscHBwdHR0eHh4fHx8gICAhISEiIiIjIyMkJCQlJSUmJiYnJyco\nKCgpKSkqKiorKyssLCwtLS0uLi4vLy8wMDAxMTEyMjIzMzM0NDQ1NTU2NjY3Nzc4ODg5OTk6\nOjo7Ozs8PDw9PT0+Pj4/Pz9AQEBBQUFCQkJDQ0NERERFRUVGRkZHR0dISEhJSUlKSkpLS0tM\nTExNTU1OTk5PT09QUFBRUVFSUlJTU1NUVFRVVVVWVlZXV1dYWFhZWVlaWlpbW1tcXFxdXV1e\nXl5fX19gYGBhYWFiYmJjY2NkZGRlZWVmZmZnZ2doaGhpaWlqampra2tsbGxtbW1ubm5vb29w\ncHBxcXFycnJzc3N0dHR1dXV2dnZ3d3d4eHh5eXl6enp7e3t8fHx9fX1+fn5/f3+AgICBgYGC\ngoKDg4OEhISFhYWGhoaHh4eIiIiJiYmKioqLi4uMjIyNjY2Ojo6Pj4+QkJCRkZGSkpKTk5OU\nlJSVlZWWlpaXl5eYmJiZmZmampqbm5ucnJydnZ2enp6fn5+goKChoaGioqKjo6OkpKSlpaWm\npqanp6eoqKipqamqqqqrq6usrKytra2urq6vr6+wsLCxsbGysrKzs7O0tLS1tbW2tra3t7e4\nuLi5ubm6urq7u7u8vLy9vb2+vr6/v7/AwMDBwcHCwsLDw8PExMTFxcXGxsbHx8fIyMjJycnK\nysrLy8vMzMzNzc3Ozs7Pz8/Q0NDR0dHS0tLT09PU1NTV1dXW1tbX19fY2NjZ2dna2trb29vc\n3Nzd3d3e3t7f39/g4ODh4eHi4uLj4+Pk5OTl5eXm5ubn5+fo6Ojp6enq6urr6+vs7Ozt7e3u\n7u7v7+/w8PDx8fHy8vLz8/P09PT19fX29vb39/f4+Pj5+fn6+vr7+/v8/Pz9/f3+/v7////i\nsF19AAAACXBIWXMAABJ0AAASdAHeZh94AAAgAElEQVR4nO3dCbxNVfsH8N+58zVT5ozR25xQ\nKQmhSdFEoTdyS0VS6R8NUqlIc9IkzdGguaRkKM1I06tyk0ITIcrs3v0/507ueIa7n73W2mv/\nvp9P5x7O3ut5uuf83Hv22XstOETkGnQ3QGQDBolIAINEJIBBIhLAIBEJYJCIBDBIRAIYJCIB\nDBKRAAaJSACDRCSAQSISwCARCWCQiAQwSEQCGCQiAQwSkQAGiUgAg0QkgEEiEsAgEQlgkIgE\nMEhEAhgkIgEMEpEABolIAINEJIBBIhLAIBEJYJCIBDBIRAIYJCIBDBKRAAaJSACDRCSAQSIS\nwCARCWCQiAQwSEQCGCQiAQwSkQAGiUgAg0QkgEEiEsAgEQlgkIgEMEhEAhgkIgEMEpEABolI\nAINEJIBBIhLAIBEJYJCIBDBIRAIYJCIBDBKRAAaJSACDRCSAQSISwCARCWCQiAQwSEQCGCQi\nAQwSkQAGiUgAg0QkgEEiEsAgEQlgkIgEMEhEAhgkIgEMEpEABolIAINEJIBBIhLAIBEJYJCI\nBDBIRAIYJCIBDBKRAAaJSACDRCSAQSISwCARCWCQiAQwSEQCGCQiAQwSkQAGiUgAg0QkgEEi\nEsAgEQlgkIgEMEhEAhgkIgEKgvTlIiJf+TLxV7n3QVoIIp9ZmPDL3PsgfYTtntcgErQdHyW8\nD4NEVAqDRCSAQSISwCARCWCQiAQwSEQCGCQiAQwSkQAGiUgAg0QkgEEiEsAgEQlgkIgEMEhE\nAhgkIgEMEpEABolIAINEJIBBIhLAIBEJYJCC48P9z92iuwdrMUiB8WYVoOM63V3YikEKimdS\nI5Ov7fez7j4sxSAFxP1JSJ/eH2j4he5O7MQgBcMEoNYHTu5YoNrbunuxEoMUBLmXA/WXRO7d\nm4S0Z3W3YyMGKQguAlpk5999Ph2hWXq7sRKDFAArk3Dgr4V/mFcdx+lsxlIMUgDcChQ7xHAp\nklbq68VWDFIA7IsDiv3pK2CctlasxSDZ71PgjuJ/boMWubp6sRaDZL+hSPmt+J/vBebr6sVa\nDJL1tu+Jk0r8xbp0DNTTisUYJOvNAJ4r+Td9UHWjnl7sxSBZrxdqljrpeybwqJ5e7MUg2W5N\nKoaU+qucJuiopReLMUi2uwdln+Krge909GIxBsl2bdGqzMHuZSGM1tGLxRgky30L3FT2b49G\ng53qe7EZg2S5KxFaXvZvpwJvqu/FZgyS3XY1Rpdy/vrf6jhDeS9WY5DsNhN4rLy/Pw9pa1T3\nYjUGyW5no8qm8v7+A+Ae1b1YjUGy2t+ZGFDuA7mt0UZxL3ZjkKz2JPBu+Y/cgNCfanuxG4Nk\ntSzU2FX+I+8DL6ntxW4MktX+g+MreGRrBi5T2orlGCSbrQ2V92lsvo5op7IV2zFINnsNmFPR\nY6OR/LfKXizHINnsKqT8U9FjbwGclksOg2SzaL++/Z2MaxW2YjsGyWLbMnBpxY+2wTHqWrEe\ng2Sxj4DnK350ONK3quvFdgySxSYCqyp+9AXgQ3W92I5BslhvNIvy6O/Arao6sR+DZK/ceugf\n7fFWOFFVK/ZjkOz1AzA52uODKzx/iBLGINnrMeDLaI8/DixR1Yv1GCR7VXzGar4fgftU9WI9\nBsle+8ZaCKkR+qjpJAAYJGutC+GG6FuchfpqWgkABslarwOzo29xP7BMTS/2Y5CsNQrJ5U7X\nsNtXwFQ1vdiPQbLW0Tg0xhY5dTBISSsBwCDZansmLom1TU+0VNFKEDBItvoYmB5rmwlRT8aj\nBDBItroDiLl4+cdl1iCjSmKQbHUaGsfcZkcVDFPQShAwSLZqgLNjb9QFB3vfSSAwSJbKBibF\n3moMQuu87yUIGCRLPQksjr3VbOAN73sJAgbJUkNQLY6lxDYlYaznrQQCg2SpA9Etns1a41Sv\nOwkGBslO21LjWyX2DLTwupVgYJDs9CXwbDzb3YTQRq97CQQGyU5PA1/Hs92rnEpIBoNkp1FI\njetbuCLGvA4UJwbJTj1xYFzb5dbChR63EgwMkp2axnNeQ0QndPC2k4BgkKy0MYSb49tyGKrm\neNtLMDBIVvoIeC2+LR/h5eYiGCQrPQQsj2/LT4EXve0lGBgkK10S9y9sm5MxxttegkFHkLZ/\nPven6FswSC51xhHxbroPennZSVAoDdK4uZHbh2oDaBd1slwGyaU9cX68m/ZBcy87CQqlQcKo\n8M2bSD/two6o+WOUDRkkd34D7ol323EIbfCyl4BQH6TWNZeGb18KnRdlQwbJnXeA9+Ld9nXg\nAy97CQjlQVqDa/Lun1p6RoEtd00ocgGD5MqdwJ/xbvtLXJfSUgzKg7QST+fdvy611IO/dmxX\npBkqXNWe4jAIdePfuDYu8K6TwFAepF01x+fdH1wnyoYPMUiutI/vqr58x+Bw7zoJDLVB6rcw\ne+3VrTaH735X9ZQoGzJIruRUxYj4tx6OTC7c55raIOWb4TjPVk36PMqGDJIry4Ap8W89BfjB\nu16CQmmQHr977IiBp3aZ4ziTG0edvIZBcuVl4NP4t/4ceMG7XoJC0ylC/0Q/gYVBcuVGhGIs\n6FLclmRc610vQWHmuXYMkit9EpvRZF9Ee79KcWGQLLRvYqfP9UVTrzoJDgbJPttSEvtd7RZg\nvVe9BAaDZJ8lcayMVNybwHyvegkMBsk+TwHfJLL9SuA+r3oJDAbJPlfFORVXkT2Q5VErwcEg\n2edEHJTYDl1wmDedBAiDZJ8m6JfYDiN4kpBrDJJ1/g7h1sT2mAp8500vwcEgWWcB8Hpieyzi\nmsyuMUjWeRBYkdge21IKrrakSmOQrDMU1XIT3GU/9PSklQBhkKzTOfHZvM/mSUJuMUjW2SPx\nS8dvRIjfcHcYJNusBu5NdJ/ngKjzDFJMDJJtZgFzEt3nywTPzqMyGCTb3Av8lug+W5Iw1oNW\ngoRBss0wVE98p+aJngxBpTBItumBdonvdDwOle8kUBgk2zSPd9HL4kagaqKfPVEJDJJltiXj\n+sT3egD4Rb6XIGGQLPMtCuaETsgc4F35XoKEQbLMKwnNaVfoV14k6xKDZJnbgHWV2K0Ghom3\nEigMkmXOx56V2e0wdJfuJFgYJMt0xpGV2e0c7CXdSbAwSJZpiIGV2W0cT1t1h0Gyyz8h3FyZ\n/V4EFkn3EigMkl2+AJ6vzH5fA89K9xIoDJJdKntBROU+x6UiDJJdxiW0oksxLdFXuJVgYZDs\nci4aVW7Hk3CIbCcBwyDZpQO6VG7Hy5EZffE3iopBskslJmzI91DCk3hRcQySVdYDEyu353zg\nbdlegoVBssonwKuV2/N34B7ZXoKFQbLKU8D/KrlrbVws2krAMEhWGYOkrZXc9Qh0FW0lYBgk\nq5yd2HrmxQ2s7IFzimCQrNIWx1V211uBvyVbCRgGySours97CfhcspWAYZBs4ubQ2/+ApyR7\nCRgGySbvAzMru+/2FFwr2UvAMEg2eRTIrvTOrXGmYCtBwyDZZBRSd1R655NxoGArQcMg2eQ0\n/KfyO1+JdK5tXmkMkk0OxMmV33kK8KNcK0HDIFkkJxNXVH7vD4C35HoJGgbJIj8DD1Z+7zXA\nnXK9BA2DZJHZwHsudt8DQ8RaCRwGySKTgZUudj8KnaU6CR4GySKXubtc/DzUF2slcBgki/TE\nQW52nwCsl2olcBgki7TGGW52fxX4RKqVwGGQ7LEjFaPd7P8d8IRUL4HDINljGTDVzf7bkzFG\nqpfAYZDs8SawwNUAzdFPqJXgYZDscTfwh6sBuuMwoVaCh0Gyx1DUcDfARagl00kAMUj26IH2\n7ga4A/hLppXgYZDs0RJnuxvg1UqtiE4RDJI1dqa6vVb8G642VmkMkjWWA4+6G2FLEm6U6SV4\nGCRrzAbmuhxiL5wr0koAMUjWeBj42eUQXXCUSCsBxCBZYzRS3c65cD7qibQSQAySNfqitdsh\nxnPa4spikKzRvvLzfhd6EVgs0UoAMUjWqIOL3A6xBHheopUAYpBs8Tdwm9sx/gnhFoleAohB\nssUXwIuuB6mPwQKtBBGDZIsZEu9vOuIYgVaCiEGyxURgg+tBuGxfZTFItrgYddwPchNC/7of\nJYgYJFsc7/YiiohpwFfuRwkiBskWrdHX/SCfAy+5HyWI3AdpxezpD02fvUKon3wMUsJy0jHK\n/SgbBI6hB5PLIC0f3gr5Wl36k1xXDFLCfgEeEhiG839XkqsgrRmcgnpnXXffU/ddd1ZdpAxe\nI9UVg5SwecC7AsMcgWMFRgkgN0GaVy95wOLcgj/kLh6QXG++UFcMUsKmyiwTNgBNBUYJIDdB\nSj2m5FOXfUyqQEcRDFLCrkNy5ZeP3W0skrYIDBM8boI0pvTlL7ukJupkkBLWHy0khnkK+J/E\nOIEjdvj7r8qvS18Wg5SwDugmMczHwGsS4wSOWJBGSX66xCAlrB4ukBiG619WEoNkh3+AW0UG\nqoWhIuMEDYNkh6+B6SIDtXN/nW0guQpSu2IaMEg6vQp8JjLQWWgpMk7QuApSUlJ6kWQGSae7\ngLUiA12L5O0iAwWMqyCNqr77UB1/tdNqOKrLDPQY8IPMSMHiKkg7Dm1f9CEgg6RVT7SRGegD\n4C2ZkYLF3cGGpZlXFt5NKEgbR30X9XEGKVH74XSZgX4D7pUZKVhcHrXbuK7w3vzxCYywCm9E\nfZxBSlBuJq6MvVVcquNSoZECRemFfVmF+uG4rKwoGzJICVoNTBYa6hCcJDRSoCgNEkqIsiGD\nlKDwO5u3hYY6A/sIjRQoboOU+8LJbQ7IF3vHy5PbzNoQ8T88t6HMlDc/Ly9yM4OUmCeB74WG\nGuV+Lv4gchuk24EqNfPFsefCNqGLI7O0l/ce6cdQ8R9XDFJCxiJpm9BQjwCC1zoHhtsg7XX8\n8kR23Tkhs9GMCg42bFxf5C4GKTHnoonUUHNlLrUNGrdBSk109d4fu+GUlTxqJ+xodJYaaiXw\ngNRYAeL6J9InCe/+eJ1qYxkkWY1wntRQuZkYKTVWgLgN0v9V4qT7P88GgyRqawg3iQ12AHqL\njRUcboP0z/H9Zy3NzpPACDNHLo36OIOUmKXA02KD9UYcB2CpFLdBiu9zoUQxSIl5A/hYbLDL\nkZkbeysqyW2Q+g0sOl1BrikGKUH3Ab+LDTYJ+E1ssMDg3N82uBxV5H6IzAQWiA0WGAySDXrj\nQLnBvgeelBstKBgkGxyEXnKDbUvGWLnRgoJBskE1XCY4WlOcIzhaQDBIFvhD9mK8LjhKcLSA\nYJAs8AnwuuBwWagvOFpAMEgWmAZ8KzjcLTz5PnEMkgWEX/nTga8FhwsG90H6vbfch+qFGKSE\nDEE9yeE+A16RHC8Q3AcpGy9KNVOEQUpIDxwuOdxfnEg/cQySBVrjLNHxamGY6HhBwCD5n8yC\n5sUcihNFxwsCl0EaOXJkFk4J38peC8YgJWI18KDogGdyIqGEuQxSs2bNGqNu+LaZXEsOg5SY\nBXJzceW7CmmcSChB/NXO/54Gos8AnaiHgF9EBwwABsn/bkJos+iAs4G5ogMGAIPkf1loKDvg\ncuBR2RHtxyD537E4UnbAnam4RnZE+wmcIrRT/gp/BikRLdFfeMS9hT+YCgCea+d7u+R/fvTA\nYcIjWo9B8r2fgUeEh7wIdYRHtB6D5Hvz5Sfrvg0os1gIRcUg+d4TQCKTc8ZjBrBYeEjbMUi+\nJ7ikS6EvgBeEh7Qdg+R7g7CX9JAbgQnSY1pOIEgrF8z+SG6ezzwMUgI642jxMffEEPEx7eY6\nSFNa5E38vf9zYi05DFJCmuG/4mMeju7iY9rNbZDuR1rXrOFZndNEZ+dkkOK3MwVjxAfthxbi\nY9rNbZBat/w57+uq/+wv1FEEgxS/5cBU8UGvRfIO8UGt5nrpy7EFd+5Il2inAIMUvznAHPFB\npwI/ig9qNbdBajS+4M6DjUX6yccgxW+qF6uQzwfeER/Uam6DdEWnnPw7vSQvNmeQ4jfGi9/C\nVkpfvW49N0FaFba8Z+dp3/zy/SunHC95TgmDFL//opn8oDkZ+D/5UW3mJkgoSbArBil+R6Oz\nB6Pui9M9GNViboJ0VgmSS2EzSPHbC4M8GPUktPFgVIvxFCGf257kybJgw1Hdg1EtxiD53DLg\nCQ+GvRtY48Gw9mKQfO5dYL4Hw74GfOrBsPZikHzuEW/moPsGeNaDYe3FIPncNUj1YlbUzSGM\n82BYezFIPtcfLT0ZtyHO82RcWzFIPnckjvVk3I6efDxlLwbJ5xoiy5Nxz5W/7tZqDJK/bQ3h\nJk8GvgGhrZ4MbCmpIGV3PKb37WKn2zFI8foOeNqTgZ+SXuLCclJBWgJMPf0G9/3kY5DiNQtY\n4MnAHwFvejKwpaSCtD0723G2uO8nH4MUrweB1Z4M/DtwnycDW4rvkfxtFNJyvBm5Gi7zZmA7\nuQ3SgnUFdz6bIdJPPgYpXn3R2qORD8IpHo1sJbdBwisFd+6oLdJPPgYpXoejh0cj98YBHo1s\nJVdByn77bVz/dp6XD68i2BWDFK96ns3keAUy5Re+sperII0vfoHsmYJdMUhx+jeEWz0a+n7g\nV4+GtpG7X+1+ew3/HZ9n4gzJGTgYpDh9C0zzaOiZXh1Yt5Pb90g9P5HrZTcGKU5vAh97NPQP\n3lwxaCse/va18O9fwusXFNmejOs9GtpGDJKvXYkMz44INMU5Xg1tIQbJ187Evp6N3RVHeja2\nfRgkX2uPEzwbOwv1PRvbPgySr+2Biz0b+xbwWYgfg+Rnm4DbPBt8OvC1Z4Nbh0Hys6+A5z0b\n/HMUnf9FMbkP0u+95T/JYJDi8xrwmWeD/wXc6dng1nEfpGy8KNVMEQYpPvd4Oh9qLQzzbnDb\nMEh+djmqejh6Ww8PCVqHQfKzUz291OFM7OPh6JZxGaSRI0dm4ZTwreR6fQxSvNrgZA9HH4U0\nLyZxtZPLIDVr1qwx6oZvm8m15DBI8aqFSzwc/WHgZw+Htwt/tfOx9cAdHg4/G5jr4fB2YZB8\n7AtAcqaM0n4Cpng4vF0YJB97GVjs4fC7UnG1h8PbhUHysTuBdbG3qry9cZaXw1tF4BShnfJX\nxDBIcfF6odfj0N7T8W3Cc+187BQc4un4F6GOp+PbhEHysYPQ29PxJwLrPS1gEQbJx6pjhKfj\nzwAWeVrAIgySf60F7va0wBLgBU8LWIRB8q+FwKueFtgEjPe0gEXEgvSv5LRQDFI8XgC+9LbC\nnp5NiGwdsSA929B1L7sxSPGYCPztbYUj0M3bAvZwHaS1k0aOCLuwseRHGgxSPIZCcgWQ8vRD\nC48rWMNtkFbULZhDP+VGuaYYpLichLYeV7gOyds9LmELt0EaUP3+OXh01ujGs+R6YpDisx9O\n97jCY0C2xyVs4TZITUc7W/GJ4yyp86FcUwxSPHKrQPZyyrLmA+94XMIWboOU+kh4iPfDd8ZI\nvi1lkOLwh/fLJa8CHvC4hC3cBqnOOMep9nj4zvSaUi05DFJcPgHe8LhEbgau9LiELdwGqXfj\nec6R7cMv+wvqyTXFIMVjOvCN1zX29fxtmC3cBumzjHbht6RNTmuDAXJNMUjxGA9s8rpGT4/P\nL7eH68+RFj3g5F6diVCvtWI9MUhxuRB1Pa8xHNU8r2EHmTMbtq7YksgAu5Yu3Bp1AwYpDsfh\nMM9r3O3pVK42UXvS6kd9Djl1sZN9IFB9crTtGKQ47IM+ntd4HfBkkWD7KA3Sp6lIRY3lHasO\nOL0aXo+yIYMUW24GrvK8yLfAs54XsYLSIJ2c+vKu1Qedk7zAcX6o2j3KhgxSbL+q+IxnSwg3\neV7ECkqDtEdkdd85OCZyf1C0Ey4ZpNg+BGZ6X6URBnlfxAZKg5Q6NnzzLy6K3L8mpdSDK/dt\nWWRP7w/s+t4zwFLvq3RCJ++L2EBpkFqcG7mtOTpye1bplX63P/FwkQH8iRTTOIQ2e19lEBp5\nX8QGSoOUlb6g8O4nqWdE2ZC/2sV2PhooqBKO678Kyvif0iBl1w4VzIF7TmrK51E2ZJBi64YO\nCqpM44rM8ZEKUnbHY3rfviHWnku7X5d/56Amr0XbjkGKbW/0U1Dlc+BlBWX8TypIS4Cpp98Q\n9xC/Rn+YQYppV5qSGe43ABMVlPE/qSBtz852nITOE4qGQYrpF+BhFXXq4EIVZXyP03H51PuK\nLl49nBMJxYXTcfnUk8AyFXX6o7mKMr7H6bh86gaEop9BL+R6JG1TUcfvOB2XTw1CYyV1wj/5\nvlNSyOc4HZdPdcHRSup86P3MEFbgdFw+1RznKKnzB3CPkkI+x+m4/GlnCq5TU6kGhqsp5G+c\njsuffgIeVVPpUJyoppC/cTouf5oDvKemUh+0VlPI3zgdlz9NAX5SU2k0UnaqqeRrnI7Ln65R\n9vIOR3a5mkq+pmU6rpgYpFjOxt6KKs3jRPrx4Bqy/nQ4ok0eI2kVEHXqNMrjJkhHlHCoYFcM\nUix1la3umpuJKxSV8jM3QUqOSAUQCv9Xs4lgVwxSDCrXGz8AvVWV8jG3v9qt7zTsy63Opg/P\nPlZyXWAGKYYvgedV1eqNA1SV8jG3QRrct+BOzyyRfvIxSDG8DCxUVesKZOSoquVfboNUd2rB\nndsll0ZgkGK4E/hLVa3JwCpVtfzLbZDSJxTcGZ0u0k8+BimGS1BDWa13gHnKivmW2yAd2jh/\nWq3P6kmuSMUgxXASJA+SRrccmKKsmG+5DdIbyWjV45QerRB6Qa4pBimWfRFtfk1Zu9IwWlkx\n33L9geyCEzIApHXhhX0KqV0kuTXOVFfMrwTObMhZvWyV8HlfDFJ0q5WebXAi2qgr5lcCQVq5\nYPZHknNxOQxSLB8Ab6urxoVk4+A6SFNa5M19sv9zYi05DFIsTwA/qKt2D/CHumo+5TZI9yOt\na9bwrM5peFKuKQYphrFK58h6E5CckMNOboPUuuXPeV9X/Wd/oY4iGKTozoXkmY2xfAfRfyXt\n5Hryk7EFd+7gB7LqHI3OCqttS8b1Csv5k9sgNSo8CflByfkKGaToGuE8leWao7/Kcr7kNkhX\ndCo4obHXSJmG8jBIUW1VvNZ4NxyuspwvuQnSqrDlPTtP++aX71855fiYq4wlgEGKainwtMp6\nF6KOynK+5CZIKEmwKwYpqjdRiWfNhYnAepX1/MhNkM4qQfIySgYpqvsA4U/Ao3sZiLbiLzmc\n/MSXLkeVXJX1vgamqaznRwySD6m++HtzCOOUFvQh10Ha9dGM6fnEemKQYjgYp6gt2AiD1Bb0\nH7dBWtScBxuUq44Ragseo2gxJh9zG6QOtUY8MCWfXFMMUlR/Kl+y6Dw0UFvQf9wGqeorcr3s\nxiBF8wnwutqKNwOb1Fb0HbdBqr9YrpfdGKRopgHfqK34PLBEbUXfcRuk4VfL9bIbgxRN+OeD\n4m/PF4DklBw2chukzaf0eWbugjxyTTFIUWUpf8fybwg3Ky7pN64XGmvCo3aqdcVRqkvuhXNV\nl/QZt0E6PKPvtWPzifXEIEXXTHR1xLgciw6qS/qM2yBleHIaMoMUxY5kjFFd8yLUVl3SZ9wG\naQ9PjuYwSFFkA4+prnkXsEZ1TX9xG6TzPTkJi0GK4h1gvuqanP8kFrdB2tB96Oyl2XnkmmKQ\nonkIWKm65jINPwX9xW2QeGGfcqOQtkt1zZ2c/zsGt0HqNzCrkFxTDFI0fbCP+qL74nT1Rf1E\n7HqkfyWv2WSQomiHE9QX7YUD1Rf1E7EgPdvQdS+7MUhR1MbF6oteiXTlv0/6iusgrZ00ckTY\nhY2ri/XEIEXzNzBRfdVHgJ/UV/URt0FaUbfgUEPKjXJNMUhRLAZmqK86HxBdAcs6boM0oPr9\nc/DorNGNudCYIi8CX6iv+jtwn/qqPuI2SE1HO1vxieMsqSP5gR2DVLHbAMm5OONVE5doqOof\nrifRfyQ8xPvhO2O6ifXEIEVzkZ5pT9ujh46yvuE2SHXGOU61x8N3pteUaslhkKLpgfY6yg5A\nMx1lfcNtkHo3nucc2T78sr+gnlxTDFIUzdFPR9kbkLRFR12/cH1hX0Y75zE0Oa2N6DUyDFKF\ntmm4iCJiGvC1jrp+4fpzpEUPOLlXZyLUa61YTwxSFEs1LZ+3SMtRd/+QObNh6wrZH/sMUoVe\nBz7WUXdTCLfqqOsXbtdHKkawKwapQnfqusSO0xZHxfWRfOZiSB4eTUAX9VOu+Inc+khnCXbF\nIFVI09FvxxnCZfui4bIuPtMcZ+spfAfwl57KvuAmSIvKPFL2byqHQarIdk1Hv/Ud5fAJN0FK\nK30a431prvvJxyBV5DtNR78d5wfgCT2VfcFNkC7BkTOL/fXMIzFMpCcGqWKvK16HebedqbhG\nT2VfcPUeadoe2G/Eq9/+sfWPb18dsR/2EFtolEGqiMYJ5vbBmZoq+4G7gw3/jK9fdPS7wQS5\nFz+DVJGhqKGr9Mk4SFdpH3B71C7n8/GDe3XpNXj8why5phikCvVAO12lRyJT8jm2DGcR8pcW\nuo5+501M+bOu2ubjLEK+ou/ot+PMBd7VVdt8nEXIV/Qd/XacX4FJumqbj7MI+cob2o5+h1XH\npdpqG4+zCPmK1uVV2uqY4tUvOIuQr2g8+u04Z6OFvuKm4yxCvqLx6LfjjEXSVn3VDcdZhHxF\n49Fvx3kG+FZfdcNxFiE/2Z6M6/RVXwi8rK+64TiLkJ98p/UM7E3AeH3VDcdZhPzkDb1LuTbg\ntA0V4ixCfnI38KfG8l1whMbqZuOl5n6i9eh3ZOKVGrk665vMdZB2fTRjej6xnhikihyn8+i3\n40zSsJ66X7gN0qLmnI5LnRaQnKspYXO42liF3AapQ60RD0zJJ9cUg1Q+vUe/HecP4C6d9U3m\nNkhVX5HrZTcGqVzf655/ZE+cr7W+wdwGqf5iuV52Y5DKpfnot+N04myrFXEbpOFXV7Lw+hVR\nHmSQyqX56LfjXKhrvmTzuQ3S5lP6PDN3QZ449vzqpGZHT96Vd3dUtIMTDFK5huk9+u049wK/\n6u3AWK5PEWqSwFG7D9NRJZPeQSoAABnISURBVBWd10fuM0iJ03z023Fm82rzirgN0uEZfa8d\nmy/2jj1TX8nddlfqYf86DFJltNR79DvvavN79HZgLLdByng6gR2bnBO5nZN20q5ygrTxmlFF\njmeQyrE9Re/R77DaGKK5A1O5DdIeSxLYMfX6vC9PRa79LxOkNf36FGnHIJVD+9FvxzkKR2vu\nwFRug3T+uAR23KtX/terMZG/2iXuTd1Hvx3nAtTW3IGp3AZpQ/ehs5dm54m946WhSTsiX3MH\n4rLhDFKitB/9zpt8RXIiUIu4DVJCS1/+1RTd8+7kXhp9ewapPNqPfjvOLGCO7h7M5DZI/QZm\nFYpjz7VDLyu499LeDFKijkNb3S2sAkovikV5eD2Sf2g/+h1WExfrbsFMDJJvbE/Btbp7cDqg\ns+4WzMQL+3zje+Bx3T04g7Gn7hbMxAv7fOMV4FPdPUTWNtd95NBMvLDPN24F/tbdgzMTmKe7\nByPxwj7f+C8a6W7BcX4GJuvuwUi8sM832kNyevVKyq0utnS9XbRd2BcVg1RWbg1coruHsMPQ\nVXcLRlJ7YV+8GKSyVprxS9Ug1NfdgpGUXtgXNwaprHfMODvnNkBycmprKL2wL24MUln3Ar/p\n7sHJOwX9A909mEjphX1xY5DKugi1dLcQ8RPwoO4eTKT0wr64MUhldUYH3S1E5FTlkszlUXph\nX9wYpLLq4TzdLeRpZ8JRePMovbAvbgxSGeuAibp7yPNfNNTdgomUXtgXNwapjA+BN3T3kGc8\nsE53DwZSe2FfvBikMqYAP+ruIc9r+meOMJHY9Uj/Sl7LzyCVcQUydunuIU828IjuHgwkFqRn\nJX9zZpDKOAkH624hX04mLou9VeC4DtLaSSNHhF3YuLpYTwxSOVqgr+4WChyK43S3YCC3QVpR\nt+BQQ8qNck0xSGVsScJY3T0UGIC9dLdgILdBGlD9/jl4dNboxqJrIjJIpS0BntPdQ4FbgA26\nezCP2yA1He1sxSfhJ7qO5KEcBqm0acCXunso8BrPtiuH2yClPhIe4v3wnTGSn3czSKWNQdIW\n3T0U+Jlz25XDbZDqjHOcao+H70yXXMuNQSrtTOytu4UidSD5maEl3Aapd+N5zpHtwy/7C+rJ\nNcUglXEATtbdQpGuutc7M5HrC/sy2jmPoclpbTBArikGqbSd6fg/3T0UuRzpO3T3YBzXnyMt\nesDJvToToV6S100ySKX8AEzV3UORJ4CvdfdgHJkzG7aukH0nzCCV8hrwse4einwJPKW7B+Nw\n7m9fmACs191DkR1pGKm7B+O4DdKhRxQ6qtdEsc/pGKRSBqKB7haKacNr+8pwG6S9agJIDv+X\nngY0+1WoKwaplMONmk1uEPbQ3YJx3M9rd+ysTc7mOccN3LnxrmSpzxcYpFJqYqjuFoq5G1il\nuwfTuA3SsK45eV9zjr3ecYZInc3IIJW0Gpiku4di5plyta5B3AapXuHsnw81d5xHUkV6YpBK\nmw28p7uHYjaE4MmcN37mel67wqsnbkt3nLFSF/cxSCXdB0i9/RTRHGfqbsE0boPUtkH+chTf\nNd/XWVhP6jQWBqmkoaiRq7uH4k5FK90tmMZtkF5Pxr4n9+11cAhTnWPSEx+rfAxSSV1xhO4W\nShiL0EbdPRjG9Qey83tkRA6AH/GS4zz2uVRXDFJJDTBQdwslvAJILj5iA4kzG9b/+Mt2mW4K\nMUglbAAm6O6hhJ+A+3X3YBieIuQDHwOv6e6hhNzauEB3D4ZxE6QjSjhUsCsGqYSpwDLdPZTU\nGYfpbsEwboKUHJEafocUCv9Xs4lgVwxSCVcifafuHkoagUzDOtLN7a926zsN+3Krs+nDs4+V\nXLqeQSrheFMmhyzyGPCt7h7M4jZIgwunLezJub89Ux/n6m6hlC+AZ3T3YBa3QapbeOHm7XVF\n+snHIBX3G3Cn7h5K2Z5m0KXvRnAbpPTC47Kj00X6yccgFfe2WWfa5TkEPXS3YBbXF/Y1zv8U\n9rN6hwh1FMEgFTfBwIXEB0LyNxALuA3SG8lo1eOUHq0QekGuKQaphH5orLuFMu4CVuvuwSiu\nP5BdcELkFKG0Lpz72zP74yTdLZQxB3hLdw9GETizIWf1slXCHyowSMVsScE1unsoY30It+ju\nwSg8Rch4C4HndfdQVlNj1msyA4NkvCnA97p7KKsX9tHdglEYJONdgkwzVo8t4Xok8TkqhkEy\nXifDrurL9xIq8cqxGINkutxaGKK7h3IsBybH3io4GCTTrTDzIrpwvs/X3YNJGCTTvWroZd3d\ncJDuFkzCIJnuBkMnGrkGSUb2pQmDZLrTDFr0srjXgDm6ezAIg2S6FjhDdwvl+hO4VXcPBmGQ\nDLcxhJt091C+FuituwWDMEiG+8C0GYSK9IPk+tt+xyAZbhLws+4eyncPsEJ3D+ZgkAx3PmoZ\nNe33bp8C03X3YA4GyXCHobPuFiqwLR2X6e7BHAyS2XZVwaW6e6hIB3TQ3YI5GCSzLQWmxt5K\njxFI36a7B2MwSGabDizW3UNFpgGf6u7BGAyS2UYjZavuHiryE3Cv7h6MwSCZ7QQcqLuFitVH\nf90tGINBMltDDNDdQsV6oaXuFozBIBltDTBRdw8VuwX4Q3cPpmCQjPYu8I7uHio2B3hddw+m\nYJCMdrvR/+ZvSsa1unswBYNktHPQQHcL0RyEbrpbMAWDZLQDcbzuFqI5H9UNnClMCwbJZNtS\ncZXuHqJ5FPhGdw+GYJBM9gkguciHuG+AR3X3YAgGyWR3ASt19xBNTg1coLsHQzBIJuuLRrpb\niI5zchVikEzW1NCJT4pcyzm5CjBIBvsVuF13D9G9zjm5CjBIBpsBfKi7h+j+BNcby8cgGexK\npG7R3UMMLdFLdwtmYJAMdjQO091CLJyTqwCDZK4dmRiuu4dY7gV+0N2DERgkc30OTNPdQyxf\nAQ/p7sEIDJK57gN+0t1DLLn1uShzHgbJXP398P6jL/bI0d2DCRgkc7XAqbpbiO0h4CvdPZhA\nW5D+yo7yIIPk5H1GM0F3D7EtA+7W3YMJtAVpVLRRGCQnb83L+bp7iEMznKK7BRMwSMYajRQ/\nfBcGovpO3T0YgEEyVme01d1CPJ7kfKsRSoPUrpgGpUfJfX92kREMkrOrGobq7iEev3IJzAil\nQUpKSi+SXHqU5ekoZlNla1jjC+Ap3T3EZR90192CAZQGaVT13Yfq+KtdDJOBaAc2zXERMoyd\nnlwdpUHacWj7HYX3GaQY/os9DF2qr5QXgHm6e9BP7cGGpZlXFt5lkGJojZN1txCfv5IwRncP\n+ik+ardxXeG9+eOjbMYgOX+FcLPuHuJ0CDrqbkE/niJkqDeB93T3EKcrkMJjQwySoa7zz7Qi\nbwAzdfegHYNkqG44WHcL8dqUiitjb2U5BslMOTUwRHcPcevgj3MwPMUgmelr4DHdPcTtWiSt\n1d2DbgySmR4BluruIW5zgBm6e9CNQTJTf9Txx8exEVsz/HFaoJcYJCPl1sOZuntIQFfsp7sF\n3RgkI30BPKy7hwSMA1br7kEzBslIt/lgAqFiPgae1t2DZgySkbqjle4WErGzBvrr7kEzBslE\nWzJ99u79LNTcrrsHvRgkE80CXtHdQ0KeA2bp7kEvBslEI5Hyt+4eEvJPBi7S3YNeDJKJDvLd\nhQknoX6wJ1xlkAz0ewg36O4hQY8AH+vuQSsGyUBPohLPil5/JuMq3T1oxSAZ6BzU2BF7K7Mc\n7a8D9uIYJPPkNsRpuntI2B3At7p70IlBMs9XwAO6e0jYCvhmjglPMEjmuR34UXcPiTsY7XS3\noBODZJ7j0Fx3C5Vwvb9OD5TGIBlnaxVcqLuHSlgC3Ku7B40YJOPM9un1pi3QVXcLGjFIxrkK\nyetib2Wey5Ac4JkbGCTjtEEH3S1UyvvA47p70IdBMs1av06lvasueuvuQR8GyTTPAh/o7qFy\nBiPzX909aMMgmWYQqvvu/KB8bwAv6e5BGwbJMLmNfbtK+NZq+K/uHrRhkAzzoR/PDypwJmoH\ndu0+BskwlyDlT909VNYLwLO6e9CFQTLLrvo4XncPlba9Ho7V3YMuDJJZ3vXT5PlljERome4e\nNGGQzJKFNF+e1pDv+xCu1t2DJgySUXbU8feHmh3RwKfH7t1ikIzyOjBddw9uPA68qrsHPRgk\no/RHFV//n2+phZN196AHg2SSzdVwtu4e3LkIyb/o7kELBskkz/v+N6NFwE26e9CCQTLJaajh\n91MD2qJ5IKdcZZAMsjET5+nuwa3JwLu6e9CBQTLIExas6fB3FfTV3YMODJJBTkTdnbp7cO1c\npK3R3YMGDJI51qb6bHmxcn0A3Km7Bw0YJHM8CLyvuwcB++E/ubp7UI9BMkdXNLbhgNcdwALd\nPajHIBnj1yRcrrsHCWvSg3h2A4NkjInAZ7p7EHERQot096Acg2SKHU1wgB3vLX5J8+20E5XH\nIJniaWCq7h6EXAh8rrsH1RgkUxyKen4/PahQ+EeSr6+qqgwGyRDv2bRQ1xBgoe4eFGOQDHEi\nqvyluwcx4R9Jp+ruQTEGyQzfJ2GY7h4EXYBQwH4kMUhmyEJStu4eBIV/JPlvPWlXGCQj/JmB\n03X3IOr8oP1IYpCMcB0q8TyY7Oeg/UhikEyweU8cprsHYVkIfam7B5UYJBM8ALyouwdhP6UG\n67MkBskAOa3RcpfuJqRlAa/p7kEhBskALwOTdPcg7vc6aLhedxPqMEj65XZAbQv/fx8Dhuju\nQR0GSb8ngWt09+CB3O4IzdPdhDIMknYbG6LB37qb8MLP1dAiMKszM0jaDQWe092DN+4CrtLd\ngyoMkm6LktHZjgv6ysg5CilBOb+BQdIspwPSlupuwivfpePggKyXxCBp9qCdRxoKjAVu0d2D\nGgySXmvqoNlm3U14Z9sByPhadxNKMEh6nef7hVyi+yQFzX7X3YQKDJJWH4ZsnwNuInD4Ft1N\nKMAg6bTzYGT+pLsJj50PnGHDBLIxMEg6jQTG6e7Bazu6A6N0N+E9BkmjKUCbbbqb8NyG/YEp\nupvwHIOkz/tp2PNH3U0osKI+Uq1fxY9B0uaHOsj4WHcTSnycgVrWfuhcgEHSZd0+CD2ruwlF\nng+hvh0LBFSIQdJkR1fget1NKDMBqGL152UMki4XAWdaeq5qeR5NRWis7ia8xCDpMQ44Igif\nUxZ5qxowwrqJKXZjkHTYngU0DcSZM7t90Qg41d7TChkkDdZ0CufoW91dqLbyIOCwH3R34RUG\nSb1vWgBH/aG7C/X+7gFkjNuuuw1vMEjKvVEdONf+ExrKseOKJOAAu+ZmLsQgKZYzIQlJt+nu\nQpfPDgaSLrZxqhcGSa05bYDqr+vuQp8dEzKBRg/bssjnbgySSj/0AtAqGJeMVuTH8Dsl1Ltx\nre4+hDFI6qwbkQpUuzlQHx+V56mm4ShVuXiZ7j5EMUiqfH9NnfD7g6zfdPdhgJ3T2oajlNTz\ncXtWzWWQ1Fj/4JHhlw6OXaK7EVPM7RkKfz9Sjp20SncnQlQHKXf57JdfnrMyxlZ2BWnty30z\nIjE6KkjLnMT0v4vqRr4pocOufGm17l4EqA3S+pH1kKfpTVHfKdgTpF+evnD/UN7/8XXWfqhf\nWbvmj2iW/3Jo0ufOuT5Pk9Ig/dYCrQeNnTjxun6NcEi0tXMsCFLOT+/cP/z4JvkvlKrnzgnA\n/B+Vseja9in53yNUO/Ss656au3Sj7pYqR2mQslJfKLi3a3JoRJQN/Rmk3PW//O+z2S9OvmHo\naR1bp6Pw9dH9xrn2nqspYPMHt5/eCLtVaXX0aYNHjrv/mbfeW/Tt8r988kpQGqQGg3ffP6tJ\nlA0rDtIDQ4zQr0++k7t3796pXbtDWrasUzsNpVVtc/Zdn+2s7HcrUFbPmnTJcS2SynwL89Ss\n3aBly/3atWvXsXvEaX0KDVD/zD9Qwf+A0iClFpsH+oa0Ug/+VLd2kSqoYF2dJeV/pw1Td/8u\n/UY+NNeWA1LqbPtu/jN3XTGg88FNqul+DitUwYFXpUFq1nf3/d7NSz2YM292kXtQwTnC/3ap\nrVHzlsW1bhdxVPifyJP69Bky5IpRN9z58AszP/jy14Csv+CtnWt+WPTeWy9MnTzh5lGjhg0Z\nXPDjv1u7fP9pWVpzNa+BLhX8E680SCNCtxec9fzv9VEnDfyooiARmUlpkDa0RfVugy4ZNrBL\nFXSK9iaSQSKfUfs50va72iRHfs9M7fBI1Mv3GSTyGeWnCG1dtnhxdqyYMEjkM2aea8cgkc8w\nSEQCGCQiAQwSkQAGiUgAg0QkgEEiEsAgEQlgkIgEMEhEAhgkIgEMEpEABolIAINEJIBBIhLA\nIBEJYJCIBJgZpIW651wiStTChF/m3gfJ+XJRBU445mmtjmH9YNc/oaJX5peJv8oVBKlCgwZp\nLM76rC9Zn0FifdYXwCCxPusLYJBYn/UFMEisz/oCGCTWZ30BDBLrs74ABon1WV8Ag8T6rC+A\nQWJ91hegM0hDhmgszvqsL1lfZ5DWr9dYnPVZX7K+ziARWYNBIhLAIBEJYJCIBDBIRAIYJCIB\nDBKRAAaJSACDRCSAQSISwCARCWCQiAQwSEQCGCQiAQwSkQAGiUiA+iDtGJ3UrvifN4xoltow\n6zdt9dePbJrWvPcn2upHXI4sffVnHlOtZtd52up/d06DlD1P/UxV+TJPt9DrT3mQlratXuIb\nub0tzrhlcGoLVVdLlq6/rjl6jhmQkvG1pvoRC5PVBalM/cew93VX1k1LfEUgmfrfVq9z/VPj\nGqTMUVO+zNMt9fpTHaSNme2z04t/I+/CbeHb5zFSU/1hmBS+fQknaaoftrPNIcqCVKb+n9UO\n/ddxsqsN1VS/P+aGb79CFzX1yzzdUq8/1UFaN3KHU+Ib2ab6tsiXVvVy9dS/rNuO8G1uZjMl\n5cvWD5sQeltZkMrUvx2zIl/UfPfLqX8EIt9/p0ZzNfXLPN1Srz8dBxuKfyO3JnfL+zoIy7XU\nL7AttaOy8qXr/5h58QaF75FK1T8+c4ezbaPC6qXqD8Q34du1SSeq7GD30y32+tMdpGXIn1xs\nLGZrqV/g3ryf+Hrqd2v4t8YgNdv/i44h7P24rvpLax+y4PcvulX5VGEDxZ5usdef7iAtxrC8\nr7fjZS31881PO3qnsvKl6j+OGY7GIFVv1nDkjHub4llN9Z3v9wfQ9GN15Us83WKvP/1BuiTv\n60S8oqV+nmnpbdcpq16q/p91TnZ0BikdT4Zvf6vWYJee+ktbNLnzjakH1FT3C0mJp1vs9ac7\nSNkYmPf1OrynpX5Y7vU4YZOy4qXrn13tF61B2iN5c+RLH6g6/l+qfocqq8O3mxs33qGoesmn\nW+z1pztI21Pyj3v2wy9a6oe/sYMxXN2/xqXrz8SYVatW/Q/9Vql7w1/i/79dct4reChUfZBU\nsv4/oa55X8/Ft2qKl3q6xV5/uoPkHFEl8i9iTqMmmuo7I3CrutJl6o9EoVFa6juXIO9t/nFY\nqaX+GhyZ97UvFqkpXvrplnr9aQzS1iU/hm8fwQ3h2wdxo6b6L2GEuspl6y99I+I5HPfGd1rq\nO4tCx25znIVJBysrX7J+i9Qfwrcb6tTYpqR2sadb9vWnOkjzR40aldwgfPOX8w0ih/B3dULv\nG88OHbRZU/29MXxUHjXnKJWpn0fde6Sy9S9DmxsvyEybp6n+y0l7XPvYLS0wWU39Yk+37OtP\ndZDGF/4mk134RP5zZbPUxsNUHTUrU7/oV6sVeurnUReksvVzHzoko+ZJn2ur//GpdVNqd39L\nUf1iT7fs64+XURAJYJCIBDBIRAIYJCIBDBKRAAaJSACDRCSAQSISwCARCWCQiAQwSEQCGCQi\nAQwSkQAGiUgAg0QkgEEiEsAgEQlgkIgEMEhEAhgkIgEMEpEABolIAINEJIBBIhLAIBEJYJCI\nBDBIRAIYJCIBDBKRAAaJSACDRCSAQSISwCARCWCQTHYWViX4ePIR3nVDUTBIJmOQfINBMhmD\n5BsMkskYJN9gkEyWF5TPTt0jtdk5K8J/7IcNQ+plHvHZ5hGNqh65OO/x5Vc0SvvP5Mi2b7XN\nqJu1IS9Iu/cgVRgkk0WCtCij0U2PjK5e7y/HGYjuN37xREbTk0ctmlGr/o7I4z073Xp9S0xx\nnAXJjW6dck6n1HCQiu1BqjBIJosE6YG288L3JmGS42Th4vDdvjgzfDsCH0Ue75TjOD+ntXCc\nE/B5+G+HIhykYnuQKgySyQrfA+3YOgcjI0GaHf7DtXg6fPsAZkQefzbycFeszMncO3JvCQre\nIxXsQaowSCbLC9JTx9RC2IhIkJaG/3Is5oZvp2B65PFvIptlYcFq9Ijc25oXpN17kCoMkski\nQboa7R+f/8mj+UHKdiJBWuAUBemXyGbDMXsZTsnbIxQOUrE9SBUGyWThIG3NbPJP+N6sCoL0\nfWSzLHy0Kv8n0j/hn0jF9yBVGCSThYO0AqdF7l1dQZBeiTzYBb/vTGsVufdROEjF9yBVGCST\nhYO0JXRo+M6Sxriw3CBFfqFblbZ/JEyRo3b9w0EqvgepwiCZLPIe6WRcOH1M7Zkpe037t5wg\nHXfqw3fvF7k7M1Rv9O0nH1sz/B6p2B66/weCg0EyWSRIa/rXrXnsAufGag1+LxOk3lh/WcO0\n/R6PbPvcQWl1B29oEv5pVGwPrd0HCoNEJIBBIhLAIBEJYJCIBDBIRAIYJCIBDBKRAAaJSACD\nRCSAQSISwCARCWCQiAQwSEQCGCQiAQwSkQAGiUgAg0QkgEEiEsAgEQlgkIgEMEhEAhgkIgEM\nEpEABolIAINEJIBBIhLAIBEJYJCIBDBIRAIYJCIBDBKRgP8HYpU8NcDzv40AAAAASUVORK5C\nYII=",
            "text/plain": [
              "plot without title"
            ]
          },
          "metadata": {
            "tags": [],
            "image/png": {
              "width": 420,
              "height": 420
            }
          }
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 101
        },
        "id": "LQlMlKiI7Saj",
        "outputId": "93057e1d-3865-4431-af68-048a0f3ed651"
      },
      "source": [
        "lamda_star = rgamma(1000,a + t, b+n1)\n",
        "\n",
        "shapiro.test(lamda_star)"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "\n",
              "\tShapiro-Wilk normality test\n",
              "\n",
              "data:  lamda_star\n",
              "W = 0.99788, p-value = 0.2356\n"
            ]
          },
          "metadata": {
            "tags": []
          }
        }
      ]
    }
  ]
}