Calc confidence interval.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# ベルヌーイ分布の信頼区間計算手法の比較\n",
    "\n",
    "Author: hagino3000  \n",
    "Date: 2017-02-01\n",
    "\n",
    "## 内容\n",
    "\n",
    "次の3つを比較する、実装はstatsmodelsを使う\n",
    "\n",
    "- 正規近似\n",
    "- Clopper Peason (EXACT法)\n",
    "- Wilson Score Interval\n",
    "\n",
    "## Why\n",
    "\n",
    "ネット広告のクリック率とコンバージョン率の値は正例が少なくゼロに近い値(0.005〜0.1)になる。  \n",
    "この時、二項分布の正規近似では信頼区間の下界がマイナスになってしまい、計算に支障が出る。  \n",
    "平均がゼロ付近にのケースに強い方式を検証したい。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Trial = 100, Success = 3 の確率分布を見る"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5,0,'P(success)')"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABIkAAAE+CAYAAAAJRUAyAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzs3Xl0VEX6xvGnOvvCFpIQBEkYFgWRTRRZVHYRgoDghowygoqDODr+GMBlQFQEZdRRB1SWKIjsgiOKLCqMooKoiAooILtsAUQgCSTp+v0R0hJJmu7QyU3C93MOx/Tt6rpPw8s58p6qusZaawUAAAAAAIDzmsvpAAAAAAAAAHAeTSIAAAAAAADQJAIAAAAAAABNIgAAAAAAAIgmEQAAAAAAAESTCAAAAAAAAKJJBAAAAAAAANEkAgAAAAAAgGgSAQAAAAAAQDSJAABAMXO73Ro8eLBeeuklnTx50uk4AAAAOIUmEQAAAbZ3716nIyg9PV1fffXVOc+zceNGvfzyyzpy5Eie6ytWrNBnn31WqDkXLVqkl19+WXPnzlVoaOg5ZwQAAEBg0CQCACCApk+frjp16mjlypWO5pg6daqaNWumJk2ayO12F3qeoUOHavDgwZoxY4bn2o8//qh27drprrvuUnZ2tt9zPv3005KkESNGFDoXAAAAAs9Ya63TIQAAKO2ysrI0f/58denSRddee62+//57LV++XI0bNy72LG63W5dccok2btyoxx57TKNGjSrUPCtWrFCbNm1Uu3ZtbdiwQcHBwZ737rrrLk2aNEkTJkzQwIED8/38U089pQ0bNuS5lpaWpvnz5yssLEy9e/c+a4YWLVpo0KBB+b63ceNGvfjii1q8eLGaNWumWbNm+fHtgMD79ttv1b59e11yySXq27ev7rzzTgUFBTkdCwAA31kAAOCzjIwM26dPH1ulShU7efJkz/VPPvnEBgUF2aFDh9ojR47Yu+66y+7atcuRjBMnTrSSrCSbnJxsb7vttjN+LVy40Oschw8ftomJidYYY5ctW3bG+3v27LEVK1a0MTExdvfu3fnOcc0113hyFPbXzTffnO/cr7zyig0LC7Px8fF2+PDh9vPPP/e8J8leffXVZ3xmxIgR9lz+12fRokW2QoUKVpLNzMzMd8yyZctsixYtbHh4uI2JibG333673bdv3xnjZs6caRs3bmzDwsJsQkKCHTx4sD127Fihs6Ho/f3vf7d16tSxkZGRNjQ01NaqVcsOGzYsTy0cPXrUTpw40Xbs2NFKsldccYXdu3evg6kBAPAPTSIAAHyUlZVlu3fvnqeJ8eKLL3renzRpkjXG2NmzZzuW8dChQzYhIeGszZfY2FiblZWV7xxut9v26tXLSrJ//etfC7zXzJkzrSTbsGFDe/jw4bNm+9e//mUl2e7duxf6+1lr7RtvvGEl2a5du+Z739zvuHTp0jzXC9sk2r17tx04cKANDg62lStXLrBJ9P7779ugoCB76aWX2vHjx9snnnjCRkZG2kaNGtmMjAzPuFdffdVKsldddZWdOHGiHTJkiA0KCrKdO3e2brfb73woHldffbV99NFHbUpKip08ebLt2bOnlWSHDx+e7/j//ve/tly5cvaSSy6xR48eLea0AAAUDk0iAAB81L9/fyvJtmzZ0n755Zc2Pj7eSrKvvvqqZ8xnn33mYEJre/ToYSXZbt265fv+nXfeaSXZBx98sMA57r33Xs8qiPT0dK/3u/vuu60k26pVK3v8+PECx+3cudNGR0fbkJAQ++OPP1prc5pGUVFRdv78+T58sxw7duyw5cqVs82bN8/TeDldbpPoyiuvzHO9sE2iO+64w9apU8d+9NFH9o477si3SZSdnW0TExNttWrV7JEjRzzXn3zySSvJPvHEE9banBVa0dHRtmnTpvbEiROecQMGDLCS7LRp0/zOB2e43W5bq1Yt26BBgwLHLF++3LpcLnvfffcVYzIAAAqPg6sBAPDBY489psmTJysxMVELFixQs2bNtGjRIkVHR+vee+/VW2+9JSnnDB2nDBkyRAsWLFC1atWUkpJyxvuHDh3SzJkzFRoaqoceeuiM97OysjRo0CBNmDBBdevW1cKFCxUeHu71nuPHj9eNN96olStXqlmzZlq7du0ZY9xut/r166djx45p+PDhqlu3riTp5MmTOn78uLKysnz+jk899ZTS0tL0+uuvKywsrMBxDRo00BdffKH33nvP57kL8tBDD2n9+vVq27ZtgWNWrFih7du3a+DAgSpfvrwkKSMjQ4sWLVLHjh01efJkSdL8+fN17NgxPfjgg54nu+3fv1+bNm1So0aNPOOKyvjx49WgQQOFh4crKipKdevWzXNm1bZt22SM0euvv57nc5s3b5YxRiNHjjxjzjfeeENNmjRReHi4KlWqpMsuu0xTpkwp1LjU1FTdc889io+PV2RkpJKTk7V161a/v4e/4wojOztbJ0+eVHx8fIFjrrnmGg0YMEATJkzQjh07AnJfAACKUvDZhwAAcH6bOHGinnzySUVGRmrBggWKi4uTJDVt2lTz5s1TcnKy+vXrp4oVK6pLly7Fni8zM1ODBw/Wq6++qoiICI0ZM0bz5s3TXXfdJWOMZ9zTTz+ttLQ0DRgwQNWqVcszx7Fjx9StWzctX75cVatW1QUXXKAHH3zQp/sPHDhQERERmjp1qlq0aKG1a9fqoosu8rz/+OOP68MPP9QFF1ygzp07a82aNZKk3bt3S5J+/vlnz7XT5f7jPld6errefPNNJScn6+KLL/aaqWfPnjp58qT++c9/qkuXLnl+H/x16aWXnnXMqlWrJEnt27f3XBs7dqzuuusubd26VUuXLtXevXvzHTds2DA9/vjjSklJ0Zw5c5SdnV0khx0vXrxYgwYNUosWLfTSSy8pNDRU27Zt06+//lroOYcMGaJx48apa9euuueee+R2u/Xll1/qo48+0p133unXuNTUVF155ZVKS0vT//3f/yk6OlrPPfecrrnmGq1bt04VK1b063sUxff97bfflJ6ers2bN2vMmDH69ddfPU/rK8gDDzyg1157TZMnT9bjjz9e6HsDAFAsnF7KBABASfbee+/ZoKAgK8nOnDkz3zEpKSlWko2IiLCffvqp1/nCwsIKfYhzYmLiGfOtW7fOXnbZZVaSjYmJscuXL7ctW7a0kmyHDh3szp07rbXWrl+/3oaEhNigoCC7adOmfLPdfvvttlWrVnbr1q1+5Xr99dettda+8MIL9t///neeOV977bVCf9/vvvsuz1wfffSRlWRfe+01r7/HkuyIESPsm2++aSXZefPmWWvzbjdLT0+3hw8f9vrr9O1guQrabnbfffdZSXb//v3WWmu3bt1q27RpY91ut+cMotWrV9vk5GQbGRnp+dznn39ub7rpJmuttcOHD88zR6ANHjzYSsr3IO1cuX/2KSkpea5v2rTJ8/uaK/fPY9CgQWfMk52d7fe4AQMG2HLlytmtW7d6rm3ZssVKsi+88IJf38PXcUePHj1rHZyeMTEx0VOfoaGh9t133/WaIVeNGjVsq1atfBoLAICT2G4GAEABfvjhB910003Kzs7WsGHDdPPNN+c7rl+/fnr00UeVnp6u5ORkrVu3rsA5o6OjFRUVVehfp/vll1/UqlUrffXVV7rooou0atUqXXPNNZo1a5batm2rZcuWqXHjxnrnnXf0l7/8RZmZmXrggQdUu3btfLO9+uqr+vjjj5WUlCSbc26hrLWe1VGZmZl5rr/44ouS5Fml87e//U3333+/Z74XX3xRAwcOVOXKlVWpUiVFRUVp0aJFnl9/+ctfJEmPPPJInusdOnTIN9+GDRsk+bayR5JuvfVW1a9fXyNGjJDb7c7z3pgxY1SpUiWvv3K3EPoiPT1dkjxbzR566CGNGTNGxhhlZ2dLktLS0pSenu4Z43a7NWTIED3zzDOSlGdcUahZs6YkqU+fPvrqq6/Oeb7XXntNISEh+a6Ocblcfo1zu92aM2eO2rdvr4oVK+rXX3/Vr7/+qpiYGMXFxemLL77w+3v4Mi45OfmsdXD6NrGZM2fqgw8+0OTJk9W+fXtdf/31evTRR739NknKWRW3cePGs44DAMBpbDcDAKAAtWvXVpcuXXT8+HE99dRTXseOGjVKW7Zs0YwZM3Tdddfps88+U2Ji4hnjUlNTA5bvggsu0LPPPqs1a9Zo3LhxqlChgiSpevXqWrp0qYYMGaLnn39ePXr0kJSzPe7JJ58scL6Czh/av3+/ypcvr+DgvP/bYK2VpHy3Rq1bt05DhgxRuXLltHjxYt18881KTU1V586dPWNyzy9q3LhxnuszZ87MN8fhw4clybPt6GxcLpcef/xx3XjjjZo1a1ae9/r166c2bdp4/fzZtrSdLjIyUlLOuU7Lly9XdHS0mjdvLinnLChJioqKUmRkpOcMpokTJ6pdu3aeOjl9XFEYOHCgvv76a02fPl3NmjXT1Vdfrb///e+6/vrrC7Udb926dapVq5YqV658zuMOHDigI0eOaMGCBVqwYMEZ7+/fv9/v7+HLuBdeeOGs288SEhI8P1955ZWen++880716dNHTz31lLp06aKWLVsWOEfFihU99QsAQElGkwgAgAKEhYVp5syZSk9Pz7MyIj/GGKWkpGjHjh1auXKlOnXqpJUrVyo2NrZIM95zzz265557zrgeFBSk5557Tlu3bvWcozRv3ryzHkSdn61bt6p69epnXM9d+ZJfk6hhw4Z67rnndPnll+uyyy7z+575yW2eHD9+3OfP9OrVS40bN9bIkSPzrASLj49XdHS018+WK1fO5/vUqFFDUk6z4+GHH9bChQs97+WuILnwwgtVo0YNLVmyRKmpqXrllVe0cuXKPOOioqIUExPj8339ERERoWnTpunhhx/Wa6+9pqlTp6pHjx5q06aNFi5cqKioKE+d5zYAc2VkZJwxn7XW0xzzxtdxknTfffepV69eZ1zPbYD6+j18HVezZk1lZmZ6zZR7wHh+Bg8erBkzZmjlypVem0RpaWlF1vwDACCQ2G4GAIAXLpfL53/chYWFacGCBfrTn/6kn376SV27dvWroRFIJ06c0J133qkFCxYoKipKCxYsUFJSkt/z/Pjjjzp48KAuueSSfO8hSSEhIfl+dtCgQbriiiv8vmdBcvNv3rzZ588YYzRq1Cj99NNPeZ509swzzyguLs7rrxkzZvh8n9xVQwMHDtSNN96oqlWrSspZWbR06VLVrFlTVapUUfPmzXXixAn17dtXQ4cO9TRPDh06pFWrVumKK644a0PyXNWrV0/PP/+8du7cqb/+9a9avny5Z2tdRESEJOno0aN5PvPSSy+dMc/FF1+s9evX69ixY17v58u4uLg4lS9fXidOnFCbNm3O+NWkSRO/voev47p3737WOvDlqWRn+zP7+eef811ZCABASUOTCACAAIqNjdV7772nihUravXq1XrssceKPcNXX32lyy+/XCkpKapQoYIWL17sdZWDN0888YQkqVu3bme8l9sA87bSIpBat24tl8ulpUuX+vW5bt266YorrtDXX3/tudavXz99/PHHXn+dvgXOl2w1a9bUkiVL1LdvX8/11157Tfv27VP//v0lST169FC5cuX0+eefe7YBSjlnJGVmZnrGnS4lJUUXXXSRxo0b59f3/qPTv7+Us0Uu936//fabpJz6jYiI0HfffecZN2fOHM2ZM+eM+fr166eMjAxPjRTEl3Eul0s33XST3njjDc82xNPt2rXLr+/h67gXXnjhrHWQkJDgOU/qj8aPHy9jjDp16lTgdztw4IC+//57XX311QWOAQCgpGC7GQAAAXbxxRdr7ty5mjBhwln/AR1Ie/bs0dNPP60JEyYoKytLjRo10ty5cws8qPpsXnrpJU2fPl2JiYm65ZZbzng/94wVX7cS+Sp3+88fz8mJjY1Vly5dNHPmTI0dO/asZ+Gc7oknntC1117reZ2UlOTzyqo1a9Zo27ZtkqTt27dLkt5++23P6pHevXvLGKPx48crOTlZXbt21QMPPKCff/5ZzzzzjBo1aqSHHnpIUs4Wtn/961+6++671blzZ/Xv319ffvmlXn75ZXXu3Fl9+vQ54/4vvPCCfvrpJw0bNkx9+/bNc0aOP9q3b6+GDRuqS5cuqlGjhlJTU/Xqq68qLCzM0xAzxqhDhw6aOnWqwsLCtGvXLi1evFivvPKK+vXrl2e+66+/XnfddZeeeeYZbdiwQV27dlVISIg2btwot9vtaWr5Om7MmDFavny5rrrqKs8B61u2bNH8+fPVvHlzTZo0yefv4eu4xo0b+/R7t3btWrVt21Y33XSTZ1XT/PnztWTJEg0fPtzrYeoTJ06U2+0+4/cPAIASyanHqgEAgHPndrvtxx9/bPv27WvDwsKsJBsSEmIffvhhm56eXqg5d+7caW+55RYryQYHB9tPPvkk33Ht2rWzkuyaNWvOOmetWrVshQoV8lx7+umnrSQ7Z84ca621u3fvtpGRkZ5HjP/yyy9nzPP5559bl8tl+/btW+C99IdHtee66qqrbGH+1yf3sfcF/TrdsmXLbIsWLWx4eLiNiYmxf/7zn/N9BPvMmTNto0aNbGhoqK1SpYodPHiwPXr0aL73f/bZZ214eLiVZD/99FO/8+d64403bPfu3W1iYqINDQ210dHR9sorr7TvvfdennG7d++2HTt2tFFRUfbyyy+3//vf/+ymTZsK/H199dVXbaNGjWxYWJiNjo62DRo0sM8//3yhxh04cMDefffdNjY21pYrV842a9bMjho1ym7bts3v7+HrOF8cOnTI/uMf/7CNGze2UVFRNioqyrZs2dK+9dZbXj+3c+dOW758eXvttdf6fU8AAJxgrP3DyYQAAKDU2L17t1q3bq1t27Z5tuyMGDHCrydznW7x4sXq1q2bMjMzFRcXpylTpig5Odnz/rZt29S5c2dlZmbq559/lsvl0oEDB8562HKbNm106NAhrVu3znNtzJgxGj58uObMmaPevXtLkm6//XZ9//336t69u0aMGJHvXA899JCee+45jR07Vv/4xz8K9T1Lm5deekn333+/NmzYUOg/WxSvvXv3qkOHDtq9e7e++eabQp0JBgBAcWO7GQAApVi1atX0n//8R59++qnuueeecz4ct2PHjmratKmaNGmip5566ozmT1JSkoKCgrR//341a9ZMgwYN8ulpXMuXL/fp/lOnTj3rmGeffVaHDh3S0KFD9eWXX2rYsGFq0qRJkR/47KQPP/xQVatWLfTWQRSfQ4cOadasWRo5cqROnDihRYsW0SACAJQarCQCAAB5nDx5stgOoz4X06ZN06hRo7R582b16tVLc+fOdTpSkXj66af18MMP6+WXX9agQYOcjgMvvvvuOzVs2FDBwcHq2bOnnnvuOVWvXt3pWAAA+IwmEQAAKLWstVq7dq2OHz+u1q1bOx2nSNSuXVv9+/fX8OHDnY6Cszh06JBWrVqlyy+/XLGxsU7HAQDAbzSJAAAAAAAAoLK7eR8AAAAAAAA+o0kEAAAAAAAAmkQAAAAAAACgSQQAAAAAAABJwU4HKCqHDx9WVlaW0zFQAsXFxenAgQNOx0AJRG3AG+oD3lAfKAi1AW+oD3hDfSA/wcHBqlSpUtHNX2QzOywrK0uZmZlOx0AJY4yRlFMfPNgPp6M24A31AW+oDxSE2oA31Ae8oT7gFLabAQAAAAAAgCYRAAAAAAAAaBIBAAAAAABANIkAAAAAAAAgmkQAAAAAAAAQTSIAAAAAAABICnY6AApm09NkP1ki/bJD9teD0vFjMkm1ZRo0ky6+VCYs3OmIAAAAAACgjKBJVALZzEzZ/30gu3CWdOy3vO9t2yS7fJEUFiHT4zaZdl1lXEEOJQUAAAAAAGUFTaISxh79Te7nHpN2bc25UKWazJXXSJVipdBw6afvZb9bIx3cLztrkuyXn8h1+2CZajWcDQ4AAAAAAEo1mkQliD32m9zPPSrt2iaVq5CzUqhlB5ng0/6YLm8t63bLfrJEdt7r0s8/yj36IbkGPSJTv7FT0QEAAAAAQCnHwdUlhD1+9NQKom1ShUpy/eNpua7unLdBdIpxueS6prNcj/9HqtdIOnlC7peekF33ZfEHBwAAAAAAZQJNohLAWiv3xHHSzq1S+YpyPfSkTEL1s37OVKos1+B/So2bS1mZco8fLfvNF8WQGAAAAAAAlDU0iUoA+9mH0g/fSMEhcj04SqbqhT5/1oSEyHXPUJnLr5Kys+WeNE52x5YiTAsAAAAAAMoimkQOs78elJ01WZJkuveRqZ7k9xwmOFhmwN+lBk2lkyfl/s9o2d9+DXBSAAAAAABQltEkcpC1Vu43J0jpx6WkOjIdexR6LuMKkmvA/0nxF0iHDsj96ljZrKwApgUAAAAAAGUZTSInrV0lfbtaCgqWq9/9MkFB5zSdiYqW675HpPAI6acfZP87PUBBAQAAAABAWUeTyCHWWrnfnSFJMp16yFRLDMi8puqFcvX7W849Ppgvu2VjQOYFAAAAAABlG00ip3y3JudpZmHhMp0Kv80sP+ayljLNr5GsW+6Uf8ueOBHQ+QEAAAAAQNlDk8gB1lq5F86SJJk218lElw/4Pcyt90gVY6R9u2UXTAv4/AAAAAAAoGyhSeSEjeukrT9JIaHndFi1NyYqWq7bB0uS7LL/su0MAAAAAAB4RZPIAe73ZkuSzNXXylSoVGT3MZdeJtOiXc49Z7wm684usnsBAAAAAIDSjSZRMbM7fpZ+/E4KCpbp1LPI72d63yFFRErbN8t+urTI7wcAAAAAAEonmkTFzH72oSTJNG4uExNb5Pcz5SvJXN8n595vT5M99luR3xMAAAAAAJQ+NImKkc3KlF21XJJkWnUotvuatl2laonS8aOy898stvsCAAAAAIDSgyZRcfr2S+nY0Zynjl3SuNhua4KC5OpzjyTJfrJEds/OYrs3AAAAAAAoHYKdvHlKSorWrFmjY8eOKTo6Wrfeeqtat26tm266SZGRkZ5xsbGxGjdunINJA8O9cpkkybRoJ+MKKtZ7m7oNpMZXSmu/kPvtaQoa9HCx3h8AAAAAAJRsjjaJEhMTdcMNN6hChQpavHixpkyZotatW0uSXn/9dSejBZz99aD0/deSJNOyvSMZXDf8We5vV0trv5DdvF6mdn1HcgAAAAAAgJLH0e1m7dq1U4UKFSRJ1apVU0hIiJNxipT9fLlk3VLtejIJ1RzJYKpeKNM65ywk97w3ZK11JAcAAAAAACh5HF1JlCsjI0NTp07VrbfeKkkKCQnRHXfcoZiYGDVv3lw33HCDQkND8/1sZmamMjMzPa+NMYqIiJAxRsaYYsnvC/vFx5IkV6sOjuZyXd9H2V8slzZvkL5dLdPkSseyOCH3974k1QZKBmoD3lAf8Ib6QEGoDXhDfcAb6gMFKeqaMNbh5SRZWVl69tlnVbduXfXq1SvP9R07diglJUXx8fEaPHhwvp+fPXu25s6d63lds2ZNjR07tshz+yPzl53ae1dPKShI1d5aJld0OUfz/Pr6yzo653WF1KyjKi9Ol3FxfjkAAAAAAOc7R5tEbrdbL7zwgmrUqKHevXvnO2bv3r164IEH9MYbbygsLOyM9wtaSZSamprnupPci+fLPWeKTL1GCnroSafjyB47quxh/aWMdLnuHS7XZS2djlRsjDFKSEjQ3r172W6HPKgNeEN9wBvqAwWhNuAN9QFvqA8UJCQkRLGxsUU2v6PbzWbMmKGIiIgCG0SSFB4eLmttgUuqQkJC8j3LyFpbYv4yub/5IueHxs1LRqaoaJn23WTfmy33uzOkxs3Pu9VEJak+ULJQG/CG+oA31AcKQm3AG+oD3lAf+KOirgfHOgMZGRlauHDhGQ2iw4cPa82aNTp58qQyMjI0bdo0NW/evMAziUo6+9uv0paNkiTTqLnDaX5nOnaXwiOkXduktV84HQcAAAAAADjMsZVER48eVXZ2toYMGZLn+rhx4/Tee+/pP//5j0JDQ9W0aVP179/foZTnzq77MuepZjVqyVSOczqOh4kqd9pqoplyNb7yvFtNBAAAAAAAfud3k+jHH3/URRdddM43jouL0+zZs/N9b8SIEec8f0lh166SJJnGJWcVUS7TsbvsRwtPrSZaJTVt4XQkAAAAAADgEL+XjowaNUpDhgzRkiVLlJ6eXhSZygx7IkNav1aSZJqUwCZRVDmZtl0lSe4P5rHXFQAAAACA85jfTaIJEyaodevWeu+99zRw4EC9+uqr+vnnn4siW+m3fq2UeVKqHC9VS3I6Tb5M+2QpOETa+pO06Qen4wAAAAAAAIf43SQqX768unfvrn//+98aMmSI0tPT9dhjj+nhhx/WihUrlJ2dXRQ5SyX73RpJOVvNCno6m9NM+UoyLdtLktwfvO1wGgAAAAAA4JRzOqm4QYMGuummm9SoUSNt2bJFkydP1v3336/Vq1cHKl+pZjd8K0ky9Rs7nMQ7c20Pybik79bI7trmdBwAAAAAAOCAQj/dbNOmTXrnnXe0Zs0a1a5dW0OGDFHTpk31+eefa8KECcrKylLLli0DmbVUsQf2Sqn7pKAgqe4lTsfxysRfINO0hexXK2UXz5fp/6DTkQAAAAAAQDHzu0n0zTff6J133tGGDRvUsGFDPfbYY7rkkt+bIK1atZK1Vm+//fb53STauC7nh6Q6MuGRzobxgel8Q06T6Mv/yfboK1M5zulIAAAAAACgGPm93Wzs2LEqX768xowZo0ceeSRPgyhXo0aNtG/fvoAELLVONYlMvUYOB/GNSaojXdxQys6WXbrA6TgAAAAAAKCY+b2S6Pnnn1fVqlW9TxocrJtvvrnQoUo7a61nJZG5uHQ0iSTJ1bmX3BvXyX6yRDb5Zpno8k5HAgAAAAAAxcTvlURVq1bVwYMH5Xa781zfsWOHdu3aJUmKiIhQcnJyYBKWRr/skH77VQoNlf50kdNpfFe/sXRhTenkCdnl7zudBgAAAAAAFKNCPd3s//7v/3To0KE81w4fPqyUlJSAhCrtcp9qptqXyISEOBvGD8YYmWtvkCTZDxfKnjjhcCIAAAAAAFBcCtUkyk/SFv2FAAAgAElEQVRiYqI2b94cqOlKNc9Ws3oNHU7iP9OstVQ5Xjr2m+xny5yOAwAAAAAAiknAmkQ//PCDoqOjAzVdqWWzs6WfvpdUeg6tPp0JCpLp1EOSZJe+I+vOdjgRAAAAAAAoDj4dXL18+XKtWLHC8zojI0MvvPCCQk5tpUpLS9P27dvVt2/foklZmuzYIqWnSZFROef7lEKmVQfZd6ZLB/ZK69ZIjZs7HQkAAAAAABQxn5pE8fHxql+/vuf1Tz/9pNq1aysqKkqSFBkZqTp16qhu3bpFk7IUsZs35PxQ5xIZV5CzYQrJhIXLXHWt7OK35V72XwXRJAIAAAAAoMzzqUlUv379PE2i999/X8nJyYqNjS2yYKWV3ZLTJDK16jmc5NyYtl1lly6QfvxOdtdWmeqlc1UUAAAAAADwTaHOJOrfv7/KlSsX6CylnrVW2rJRkmRqXexwmnNjKsfJNGkhSbLL3nU4DQAAAAAAKGqFahK1bt1aYWFhgc5S+h06IP16SAoKkpJqO53mnJkO10uS7KoVskePOJwGAAAAAAAUpYA93cztduvbb78N1HSlkuc8ogv/JBNaBppotS6WkupIWZmyKz5wOg0AAAAAAChCAWsSZWRkaPTo0YGarnQqI1vNchljfl9NtPx92axMhxMBAAAAAICi4tPB1YcOHVJ4eLgiIyOVmpqa75i0tLSABiuN7KkmkUr5odWnM5e1lJ2bIv16SHbNpzJXtnU6EgAAAAAAKAI+NYmGDBmixMRE/fOf/9SgQYOKOlOpZDPSpV1bJZWdlUSSZIJDZNp0kV3wpuyyd2Wbt5ExxulYAAAAAAAgwHxqEg0bNkzR0dGe188++6wiIiKKLFSptG2T5HZLMbEyMbFOpwkoc/W1su/NlrZvlrZskGrXdzoSAAAAAAAIMJ+aRHXq1MnzOjY2VpGRkUUSqLSynvOIys5Ws1ymXAWZK9vIfrJE7mX/VRBNIgAAAAAAyhy/D65+9NFHFR4eXhRZSrXfzyMqO1vNTmfad8v54esvZA/udzYMAAAAAAAIOL+bRJdeeqlcroA9FK1MsNZKP/8oqWydR3Q6Uy1RqtdIsm7Z5YucjgMAAAAAAALM727PwoUL9cknn5xx/ZtvvtG8efMCEqrUSd0nHT8qBYdI1ZOcTlNkXO2SJUn2kyWyJ044nAYAAAAAAASS302iDz74QHFxcWdcr1ixoj788MOAhCpt7LbNOT9UT5IJDnE2TFFq2EyKrSIdPyq7eoXTaQAAAAAAQAD53SQ6fPiwEhISzrgeExOjI0eOBCRUqbM9p0lkEms5HKRoGVeQTNsukiT74bs52+wAAAAAAECZ4NPTzU4XGRmpvXv3qmLFinmu79mzx+8nnqWkpGjNmjU6duyYoqOjdeutt6p169basmWLJk+erD179igpKUn33nuv4uPj/Y1abOyOLTk/JNZ2NkgxMK06yr7zlrR7u/TTD9JFDZyOBAAAAAAAAsDvlUQtWrTQlClTtH//70+42rdvn1JSUtS8eXO/5kpMTNTo0aP1xhtv6Prrr9eUKVN07NgxPfXUU2rfvr0mTpyopKQkjRkzRm6329+oxcJaK23PaRKV9ZVEkmSiomVatJUkuT961+E0AAAAAAAgUPxeSdSnTx+NHj1aDzzwgKpVqya32609e/YoMTFRt912m19ztWvXzvNztWrVFBISotWrVys6Olrt2rWTMUZhYWEKCgrS+vXr1aBBCVy1krpPSjuWc2j1BTWcTlMsTNtk2RUfSN+skj24X6ZyyV3lBQAAAAAAfON3kyg8PFwjR47UF198oY0bN8paq+7du6tly5YKDvZ7OklSRkaGpk6dqltvvVXbtm1TnTp1ZIzR119/rYYNG2rPnj3atm1bvk2izMxMZWZmel4bYxQRESFjjIwxhcrjD7v990OrXSGhRX6/ksBUT5St10h2w7eyyxfJ1buf05F8llsTxVEbKF2oDXhDfcAb6gMFoTbgDfUBb6gPFKSoa6JQXR2Xy6WWLVuqZcuW5xwgKytLzz//vJo3b642bdrolVdeUcWKFZWVlaWNGzeqT58+Wrp0qdLS0vL9/Pz58zV37lzP65o1a2rs2LGKjY0952y++DV1r45KiqrfUDFVqxbLPUuC9N63K/WJh2RWLlOVux6UKzzc6Uh+ye/wdUCiNuAd9QFvqA8UhNqAN9QHvKE+UNwK1SQ6efKk9u7dq7S0tDPOCqpfv77P87jdbr344ouqU6eOevXqJSlnpZIkLVmyRB06dJCU80S1WrXyP++nZ8+eSk5O9rzO7aqlpqbmWWFUVLLXfytJSo+tqj179hT5/UoKe2EtKbaK3Kn7tOedmXJdfa3TkXxijFFCQoL27t3L09mQB7UBb6gPeEN9oCDUBryhPuAN9YGChISEFOmiGL+bRCtWrNDEiRMLbMDMmjXL57lmzJihiIgI9e7d23MtKSlJH3/8scqVK6f4+HhlZGRoy5YtuvHGG/OdIyQkRCEhIWdct9YW+V8ma63s9twnm9U6v/7yGpdM2y6yc1Lk/vBdqXXHUrUUsjjqA6UTtQFvqA94Q32gINQGvKE+4A31gT8q6nrwu0n01ltvqW/fvmrdurWio6MLfeOMjAwtXLhQL774Yp7rl19+uaZNm6bmzZsrKytL06dPV3x8vOrVq1foexWZ8/DQ6tOZVh1l33lL2r1d+ul76aJLnY4EAAAAAAAKye8mkdvtVpMmTc6pQSRJR48eVXZ2toYMGZLn+rhx4/TII49o0qRJmjNnjpKSkjRs2DC5XK5zul9RsNt+P7TaBJ+5mqmsM1HRMi3ayq74QO6PFiqIJhEAAAAAAKWW302iTp06ac6cObrvvvvO6cZxcXGaPXt2vu/FxsZq9OjR5zR/sTj1ZDOTmP95SecD0y5ZdsUH0jerZA/ul6kc73QkAAAAAABQCH4vz6lfv75SU1P15JNP6ttvv9X69evz/Dqf2J0/5/xQ4zxuEl1QQ6rXSLJu2Y/fdzoOAAAAAAAoJL9XEo0aNcrz83fffXfG+/4cXF3q7dwqSTIX/snhIM5ytUuWe8O3sp8ske12q0xYmNORAAAAAACAn/xuEr388stFkaPUsb8dlo4ekYzrvDy0Oo+GzaTYKlLqPtnVK2Su6uR0IgAAAAAA4Ce/m0RxcXFFkaP02bUt579Vqp73K2eMK0imbVfZOVNkP3xXtnVHGWOcjgUAAAAAAPzg95lEbrdbc+bM0eDBg3XbbbcpNTVVkjRlyhR9+umnAQ9YUtmd2yRJplqSozlKCtO6gxQaJu3eLv30vdNxAAAAAACAn/xuEs2ZM0efffaZunTpkuex9A0bNtSKFSsCGq5Ey11JVD3JyRQlhomMlmnRVpLk/vBdh9MAAAAAAAB/+d0kWr58ue6//35dd911Cg7+fbdaUlKSfv7554CGK8nsqSaRubCms0FKENMuOeeHtatlD+53NgwAAAAAAPCL302i9PR0hYaGnnH98OHDys7ODkioks5mZUl7dua8YCWRh7mghlSvkWTdsh+/73QcAAAAAADgB7+bRE2bNtXUqVOVlpbmuXbo0CFNmjRJzZo1C2i4EmvvLik7S4qIlGI4yPt0rlOriewnS2RPnHA4DQAAAAAA8JXfTaI77rhD+/fv14ABA5SWlqann35a999/vzIzM3X77bcXRcYSJ3ermaol8RSvP2rYTIqtIqUdk1213Ok0AAAAAADAR8FnH5JXhQoVNHbsWH366aeeM4i6du2q1q1b57sNrUzynEeU5GiMksi4gmTadpWdM0X2o4WyV3WikQYAAAAAQCngd5NIkkJDQ9WuXTu1a9cu0HlKBbt7W84PnEeUL9O6g+w706Xd26WfvpcuutTpSAAAAAAA4Cx8bhKlpqZq4cKF+uGHH5SamqqMjAyFh4crISFBTZs21XXXXafo6OiizFpy7NwmSTLVkhyNUVKZyGiZFm1lV3wg94fvKogmEQAAAAAAJZ5PTaKdO3fqscceU1xcnC677DLFxcUpODhYJ06c0L59+/TZZ59p2bJlGjlypKpWrVrUmR1ljx6RjhySjJGqJTodp8Qy7ZJlV3wgrV0te3C/TOV4pyMBAAAAAAAvfGoSzZgxQ5deeqkeeuihfN+/7bbb9PTTT+utt94qcEyZkXtodVyCTHiEo1FKMnNBDaleI2nDt7Ifvy/Tu5/TkQAAAAAAgBc+Pd1s48aN6tGjR8GTuFzq1auXfvjhh4AFK6k85xGxiuisXO2SJUn2kyWyJ044nAYAAAAAAHjjU5Po+PHjZ91GVq1aNR0/fjwgoUq0X3ZKkgxNorNr2EyKrSKlHZNdtdzpNAAAAAAAwAufD66eNGmSgoMLHp6VlRWQQCWd/WVHzg9VL3Q2SClgXEEybbvKzpki++G7sld1kjHG6VgAAAAAACAfPjWJ6tWrp8OHD591XP369c85UElmrf19JdEFNRxOUzqY1h1k/ztD+mWHtH6tdEkTpyMBAAAAAIB8+NQkGjlyZBHHKCV+PSSlH5dcLqlKNafTlAomMjqnUfThu3IvWaAgmkQAAAAAAJRIPp1JhFNyt5rFV5UJCXE2Syli2neTjEta/43s7u1OxwEAAAAAAPmgSeQHu+dUk4itZn4xcQlS0yslSXbpAofTAAAAAACA/NAk8gfnERWaq2MPSZJdtUL2yNnPtwIAAAAAAMWLJpEfeLJZ4ZlaF0u1LpaysmSXv+90HAAAAAAA8Ac0iXzEk83Onatjd0mSXf6+7IkTDqcBAAAAAACno0nkK55sdu4aXylVjpeOHZX94mOn0wAAAAAAgNPQJPIVTzY7ZyYoSKbD9ZIku/QdWbfb4UQAAAAAACAXTSIf8WSzwDCtO0gRUdK+3dJ3XzkdBwAAAAAAnBLs7wfefvtttW/fXhUqVDjnm1trtWzZMqWkpKh///5q3769JOmmm25SZGSkZ1xsbKzGjRt3zvc7J7nnEXFo9Tkx4ZEyV3eSXTxf7qULFNTocqcjAQAAAAAAFaJJ9MUXX2ju3Llq1qyZOnXqpAYNGhT65qNHj1a5cuVUvXr1M957/fXXCz1vUfA82YyVROfMtEuWXfZf6cfvZLdvlkms7XQkAAAAAADOe343iZ555hlt3rxZy5Yt0zPPPKOYmBh16NBBbdq0UXR0tF9z9e/fXwkJCRo5cqS/MTwyMzOVmZnpeW2MUUREhIwxMsYUet7Tnf5kM1e1xIDNe74yleNlm10lu2q57Afz5Bo4rPjuferPjj9D/BG1AW+oD3hDfaAg1Aa8oT7gDfWBghR1TRhrrS3sh9PT0/Xpp59q2bJl2r17t1q2bKkePXroggsu8GuekSNH6qqrrvJsN7vtttsUHBysmJgYNW/eXDfccINCQ0Pz/ezs2bM1d+5cz+uaNWtq7Nixhf1K+cpK3a89d3SRXEGq/vYnMiH5Z4HvTm7brH2DbpGMUcKr8xRSjRVaAAAAAAA4ye+VRKcLCQmRMUYZGRmSpMOHD2vo0KFq3769br/9drlchTsXe/r06crKytKOHTuUkpKiAwcOaPDgwfmO7dmzp5KTkz2vc7tqqampeVYYnQv3hm9zfoiror2pBwMy53kvLEqm4eWy677UvjdfVdDt9xXLbY0xSkhI0N69e3UO/VGUQdQGvKE+4A31gYJQG/CG+oA31AcKEhISotjY2CKbv1BNorS0NC1ZskSLFi3SyZMn1bFjR3Xt2lUVKlTQgQMHNGbMGM2bN0833nhj4YMFB+tPf/qTBg0apAceeEB33323wsLCzhgXEhKikHweSW+tDdhfJrsnZ6uZEqrzFzSAzHW9ZNd9KfvZh3J3u0WmYuViu3cg6wNlC7UBb6gPeEN9oCDUBryhPuAN9YE/Kup68LtJ9Oabb2rZsmUKCQlR165d1alTpzxPIouLi9OAAQP00ksvnVOTKFd4eListc7uxdy7W5JkqlRzLkMZZGrXl2rXlzavl132X5nef3E6EgAAAAAA561CPd3s1ltvVdu2bQs8J+jCCy/UsWPHChXo8OHD2rJlixo2bCi3261p06apefPmBd6rONi9u3J+SKBJFGiu63rJ/dJ62eUfyF53o0yUf4efAwAAAACAwPC7SfTiiy+e9ayhsLAwPfbYY2eda+jQodq3b58yMjK0efNmTZs2Td26ddP333+v//znPwoNDVXTpk3Vv39/f2MGVu5KoqrVnc1RFl3aTKqWKO3eLrv8fZmuNzmdCAAAAACA85LfTSKXy6WVK1eqWbNmec4IWrt2rU6ePKkrrrhCISEhqlOnzlnnKugpZL169fI3VpGxJzKkQwdyXlShSRRoxhiZ63rLTvqX7IfvynboLpPP2VMAAAAAAKBoFerxY5MmTdLRo0fzXCtfvrzmzZsXkFAlyr5fcv4bXU6mXHlns5RRpllrqXK8dPSI7GfLnI4DAAAAAMB5qXDPqM9HpUqV9MsvvwRquhLDcx4Rh1YXGRMUJHPtDZIku3i+bFaWw4kAAAAAADj/BKRJlJ6ernfffVcJCQmBmK5kyT2PKIGtZkXJtGovlasgHdwv++UnTscBAAAAAOC849OZRHPmzNHcuXPzXBs0aFCe19HR0frb3/4WuGQlBU82KxYmNEymw/Wy86fJLpor2/xqGVeQ07EAAAAAADhv+NQkatOmjS655BJJkrVWzzzzjO655x5VrFhRkhQZGamqVavmOci6rMjdbsZKoqJn2nSRXTxf2rNTds1KmSuudjoSAAAAAADnDZ+aRHFxcYqLi/O8drlcqlu3rmJjY4ssWElg3e7fD65mJVGRM5FRMh27y74zXXbhLNlmrVhNBAAAAABAMSnUmUTPPvusYmJiAp2l5Dl8UDp5QgoKkmLL4HlLJZBplyxFRnlWEwEAAAAAgOJRqCZRbGysXK6APRit5Np36jyiuKoywT4tusI5yl1NJClnNZE72+FEAAAAAACcHwLW6cnMzNRbb70VqOlKBLsn58lmbDUrXqZdt99XE331mdNxAAAAAAA4LwS0SfTOO+8EarqSYR+HVjshz2qid2eymggAAAAAgGLg0x6q1atXq2LFiqpbt65WrFiR75gTJ04ENFhJYPeyksgppl032aXveFYTmcuvcjoSAAAAAABlmk9NosmTJyspKUnDhw/X+PHjVblyZRljzhhX5p52toeVRE75/Ulnb+WsJrqsJU86AwAAAACgCPnUJHr55ZcVFPT7P9DHjRunyMjIIgtVEtiMNOnXgzkvWEnkCFYTAQAAAABQfHw6kygkJOT8eJrZ6fb9kvPfchVkoso5m+U8ZSKjZDpwNhEAAAAAAMXB787PjBkzyvwqIum084iqsIrISab9aU86W/U/p+MAAAAAAFBm+d0kOm9WFO09dR5RVc4jcpKJjJLp3FuSZN+ZLpuZ6XAiAAAAAADKpoB1fNxutzZv3hyo6ZzHk81KDNMuWaoYIx3cL7vifafjAAAAAABQJgWsSZSRkaFHHnkkUNM5zuauJKrCSiKnmbAwmev7SJLse7Nl0447nAgAAAAAgLLnPNk75h/rdv9+cHVVVhKVBKZleymhunTsqOyS+U7HAQAAAACgzAn2ZdB999131jHW2nMOU2IcOiBlnpSCgqXKVZxOA0kmKEiunn+We8LTskvfkW3TRaZijNOxAAAAAAAoM3xqEh04cEDdu3dXSEhIUecpGU5tNVN8VZmgIGez4HdNrpT+dJH084+y782Sue1epxMBAAAAAFBm+NQkkqQePXooMjKyKLOUGJZDq0skY4xcN9wh97iHZf+3WLZDd5kqFzgdCwAAAACAMsGnM4kaNWqk4GCf+0mlX+6h1QkcWl3SmIsaSJc2k9xu2QVvOh0HAAAAAIAyw6cm0cMPP6zQ0NCizlJisJKoZHPd8GfJGNk1n8pu3eR0HAAAAAAAygSebpafU00iVhKVTKZ6TZkr20iS3LMnl61D0wEAAAAAcEjAmkQZGRl69NFHAzWdY2x6mnTkUM4LVhKVWKZHXyk0VNq8XvbLT5yOAwAAAABAqRewJpHb7damTWVg60/uVrPyFWUio53NggKZmDiZ63pLkuzc12VPZDicCAAAAACA0s2n06jfeOMNJSQk6Nprr9X48ePzHZOVleX3za21WrZsmVJSUtS/f3+1b99ekrRlyxZNnjxZe/bsUVJSku69917Fx8f7PX9h2FOHVoutZiWe6dRT9tNl0sH9sovm5qwuAgAAAAAAheLTSqJ9+/YpNTVVkrRixYp8G0LBwcG65ppr/Lr56NGjtWHDBlWv/ntD5tixY3rqqafUvn17TZw4UUlJSRozZozcbrdfcxea5zwitpqVdCY0TK6b7pQk2cXzZQ/sdTgRAAAAAACll08rif7xj3/keT1gwABFRkae88379++vhIQEjRw50nNt9erVio6OVrt27WSMUVhYmIKCgrR+/Xo1aNDgnO95Np6VRFVoEpUKTVpIFzeUNq6Te26Kgu4d7nQiAAAAAABKJZ+aRKe76KKLFBQUFJCbJyQknHFtx44dqlOnjowx+vrrr9WwYUPt2bNH27Zty7dJlJmZqczMTM9rY4wiIiJkjJExxv9Q+06tJKp6YeE+j2JljJG59W5lP36/9PXnshvXyVWvkdfxp/8XyEVtwBvqA95QHygItQFvqA94Q32gIEVdE343iUaNGlUUOTwyMjJUsWJFZWVlaePGjerTp4+WLl2qtLS0fMfPnz9fc+fO9byuWbOmxo4dq9jYWL/vbbOztWv/HklSlYZNFFy1auG+BIpX1ao63KW3ji2cLdfcKUp46S2ZIO+lnV+DEpCoDXhHfcAb6gMFoTbgDfUBb6gPFDe/m0SSdOLECW3atEm//fabKlasqLp16yo4uFBTnSE8PFyStGTJEnXo0EGSdPjwYdWqVSvf8T179lRycrLndW5XLTU1Nc8KI1/YA3ulzJNScIj2uyWzZ09hvgIcYDv0kJYvUtb2n/XLjClyte+W7zhjjBISErR3715Za4s5JUoyagPeUB/whvpAQagNeEN9wBvqAwUJCQkp1KIYX/nd2fn+++/1/PPP69ixY55r5cuXV79+/dSqVatzDpSUlKSPP/5Y5cqVU3x8vDIyMrRlyxbdeOON+Y4PCQlRSEjIGdettX7/ZbJ7Tp1HFF9VMi7+MpYmUdEy3fvKTp8g9zvTpctayVSoVODwwtQHzg/UBryhPuAN9YGCUBvwhvqAN9QH/qio68Gnp5udbvLkyWrRooUmTJig6dOn64UXXlDbtm01fvx4rVmz5pwDXX755dq1a5dCQ0OVlZWl6dOnKz4+XvXq1Tvnuc/Gc2h1QnXvA1Eimas7STVqSWnHZWdOdDoOAAAAAAClit8riQ4fPqwePXooJiZGklS1alX16dNHkZGRmj17tpo1a+bzXEOHDtW+ffuUkZGhzZs3a9q0abrlllv0yCOPaNKkSZozZ46SkpI0bNgwuVx+97P8t/fUodUJPNmsNDKuILnuuE/upx6SXfOp7JVtZBpd4XQsAAAAAABKBb+bRM2bN9eyZct0yy235LnesmVLzZo1y6+5xo4dW+B7o0eP9jfaOWMlUelnatSS6dhddvF8uae/IlfdBjIRkU7HAgAAAACgxPN7ec5f/vIX/fTTT3r99dd18uRJz/WsrCy53e6Ahit2+1hJVBaYbn2kuATpcKrs/GlOxwEAAAAAoFTwaSXRzTfffMa1H374QYsWLQp4IKfYtOPSkcM5L6rQJCrNTFiYXH3/Kvfz/5Rd/r5s82tkal3sdCwAAAAAAEo0n5pEI0aMKOoczju1ikgVKslERjmbBefM1G8s06Kd7OcfyT31Zbkee14m+Myn4AEAAAAAgBw+NYnq169f1DkcZ08dWs0qorLD3HSn7PdfSb/skP1gnkzyLWf/EAAAAAAA5ym/D65+/PHH872elZWlgwcPavz48eccyhGnDq02HFpdZpjo8jI3D5Cd9C/Z92bLXtZK5oIaTscCAAAAAKBE8vvg6vXr16tKlSqqX7++6tevr1q1amnbtm0KDg5Wnz59iiJjsfCsJOLQ6jLFXHG1dGkzKStL7knPyWZlOh0JAAAAAIASye+VRJJ04403qnLlyp7XnTt31ujRoxUeHh6wYMWOJ5uVScYYuf48SO6Rg6UdW+T+7wxp0FCnYwEAAAAAUOL4vZIoMjJSR44cyXMtNjZWd999t2bOnBmwYMXJurOlfb/kvGC7WZljKlWW6/ZBkiS7aK5OfP+Nw4kAAAAAACh5/G4SXXrppZo6daoyMjLOeO/AgQMBCVXsDh6QsjKl4BCpcpzTaVAEzGWtZFq0k6zVwX/9UzbtuNORAAAAAAAoUfzebnbHHXdoxIgRuu+++9SkSRPFxMToyJEj+vzzz3XZZZcVRcaid+rQalW5QMYV5GwWFBlz692ym75X9v49MjNek+vOB5yOBAAAAABAieH3SqLKlSvr2WefVZcuXbR//3598cUX+uWXX9SzZ08NHDiwKDIWOc+h1VU4j6gsMxGRCur/d8nlkv38I7m//NTpSAAAAAAAlBiFOrg6IiJCN9xwg2644YZA53HG3txDqzmPqKwzdS5R+Rv76bdZU2TfHC9b62KZmFinYwEAAAAA4LhCNYk2bdqkb775RkeOHFGlSpXUrFkzJSUlBTha8bG52814stl5oXyfu/Xbqv9J2zbLPeV5uf4+im2GAAAAAIDznt/bzRYuXKhHH31UK1eu1Pbt27V06VINHTpUEydOLIp8xWNf7koimkTnAxMcrKABD0mhYfr/9u48Psrq7vv491yZTBaGEJIYI2tQqYASUURaZRNQ8CW1KlatK1Zc0VrbW+UutUJ5VCwqBdQ+ersEeVDZqXUDQSiKCzeiCAKV1bAkhiwQEhKyzHn+GBiJJGMSZjLJ5PN+veY113rmd4Ufh8mPc51L/1kv+6/Z4Q4JAAAAAICwq/dIon/+858aM2aMBgwY4N/21Vdf6eGDNzMAACAASURBVNlnn9VJJ52kK664IqgBhpo9VCIdKPStMCdRi2HSOsjcdI/sy1Nk35kte3p3mTPPCXdYAAAAAACETb1HEkVFRen000+vtq1Xr1667bbb9MEHHwQtsEZzZBSR2rSViW8V3ljQqJyfXyQzYJhkrbwvPS1bkBfukAAAAAAACJt6F4kuvfRSzZw5U1VVVdW2n3HGGcrPzw9aYI3F/2QzJq1ukcx1t0sdu0jFRfK+8KRsRUW4QwIAAAAAICzqVCTauHGj/9W1a1cdPnxYjzzyiL788kv/9vXr18taG+p4g+/IpNWGW81aJBPtlnPXWCm+lbT9P7JvvNA88xgAAAAAgBNUpzmJJkyYUOP2SZMmBTWYcPhhJBFFopbKpJ4i5/YH5Z32V9mPlkidTpMZdGm4wwIAAAAAoFHVqUg0e3YEP/3p6Egibjdr0cxZ58pceZPsghmyb/6PbLuOMj87K9xhAQAAAADQaOo9J9FR+/bt09q1a7V27Vrl5TXPCX+tt0rKzfatMJKoxTPDr5Lp01+qqpT3+Sdkc/eGOyQAAAAAABpNnUYSHauiokIvvPCCPvroI/82Y4wGDhyo22+/XS5XvZsMn/x9UmWF5IqWkk8KdzQIM2OMdMvvZPflSDu3yDt9opyxk2VaecIdGgAAAAAAIVfvkURvvPGGduzYoQkTJmjmzJl67bXX9Kc//UmbNm3Sm2++GYoYQ+fIrWY6uZ2MExXeWNAkmJgYOff+WUpKkXL2yPt/J8lW8sQzAAAAAEDkq3eRaNWqVbrjjjvUrVs3ud1uxcTEKCMjQ3fddZdWrlwZihhDxj9pNU82wzFMm7Zy7ntEiomTNn8t++o0Wa833GEBAAAAABBS9S4SlZSUKCkp6bjtSUlJKi0tDUpQjYZJq1EL06GLnLvHSlFRsqv/LTt/RrhDAgAAAAAgpOpdJOrWrZsWLlwo7zEjKyorKzVnzhx169YtqMGFmn8kEZNWowbmzHNkbvmdJMkuWSjvkoVhjggAAAAAgNCp9yzTo0aN0oQJE/T1118rPT1dkrRt2zZVVVVp/PjxQQ4vxL73FYkYSYTaOL+4SN79BbILZsjOfVXeuFZy+l8S7rAAAAAAAAi6eheJOnTooMmTJ+v999/Xjh07JEkDBw7U8OHDlZiYGPQAQ8WWFEsHCn0rjCRCAGb4VVJJkezihbIzn5PXHSOn78BwhwUAAAAAQFA16Hn1iYmJuu6664IdSzXXXHON4uPj/espKSl66qmngvcBe7N870kpMnHxgY9Fi2aMkUaOkspKZf/9vuwrU2TdMTLn/DzcoQEAAAAAEDT1LhI9//zzkiS3263Ro0cHPaBjZWZmhqxtu+c730K7ziH7DEQOY4x0/V3S4TLZz1bI+8KTcu54SObcX4Q7NAAAAAAAgqLeE1f/+9//VkJCgpKTk0MRT+PZ6ysSmfadwhwImgvjODKj7pc5f6BUVSXvi3+T/eKTcIcFAAAAAEBQ1HskUUpKigYNGqQOHUI72XN0dLRuueUWJSUlqW/fvrrqqqvkdruPO66iokIVFRX+dWOM4uLiZIzxjf6ohd3ju93MtE8PeBwiy9E/64b+mRuXS/a2B+Q1RvbzFfK++Dc5o+6Xc8HgYIaJMDjR3EBkIz8QCPmB2pAbCIT8QCDkB2oT6pww1lpbnxPefvttbdiwQQ8++KCioqJCFZckqbKyUllZWXr11VeVmpqq++6777hj5syZo3nz5vnXu3TpoieffDJgu9Za7b1+qLxFB3Ty1P8n9+ndgh47IputqlLB1L/q0LJ3JEltfnu/EkbeFOaoAAAAAABouHoXiTZu3Kj58+ervLxcI0eOPG50T48ePYIaoCTl5OTo97//vWbMmKGYmJhq+2obSZSXl1dt+7HsgUJV/fFmyRhFPTdXxh1T43GIPMYYpaWlKScnR/VM/eNYr1feua/KfrDI1/YlV8i5+lYZp953caIJCGZuIPKQHwiE/EBtyA0EQn4gEPIDtYmOjlZKSkrI2q/37WYTJkzwLz/xxBPH7Z89e/aJRVSD2NhYWWtrHFYVHR2t6Ojo47Zba2v9y+SftPqkNCnazV+6FihQftSZMXKu+a28bRJl52XKLlkkb9EBmVvuk3E16MGBaAKCkhuIWOQHAiE/UBtyA4GQHwiE/MCPhTof6v2bbCiKQD9WWFiobdu2KSMjQ16vVzNnzlTfvn1rnJOoIXiyGYLJGXaVvK3byM6YLvvZctniIjl3PSwTExvu0AAAAAAAqLMmeV+M1+vVO++8ozvvvFP333+/3G637rrrruB9wN6jk1bzZDMEh3PBEDn3jJPcbmnDF/I++bBsfm64wwIAAAAAoM7qPJJo586dWr16tYqKitSlSxf1798/aCN7fiw5OVmPPvpoSNqWjhlJ1J6RRAgec3YfOQ9MlPf5x6VdO+T9P3+Qc+dDMt0ywh0aAAAAAAA/qU4jiVatWqWxY8fqk08+0Y4dO/Taa6/pwQcfVF5eXqjjCzprrbTnyEgibjdDkJnTu8sZ94zU6TSpuEjeKX+Rd+lb3EcMAAAAAGjy6lQkevvtt3XllVfq73//ux577DE999xzSk5OVmZmZojDC4GCfdLhUinKJZ18SrijQQQyySfJeXiSzM8HSV6v7OyXZF/9u2z54XCHBgAAAABArepUJNq9e7cGDhzoX/d4PLr55pv11Vdfyev1hiy4kDh6q1laexnX8U9FA4LBuGNkfvuAzLW3SY4j++lyeZ8cK5uzJ9yhAQAAAABQozoVicrLy5WQkFBtW6dOneT1elVUVBSSwELFHp20uh2TViO0jDFyhv5Kzu8nSJ7WUtY2eSf+Xt6V73P7GQAAAACgyanzxNVz5syRy1X9cMdxNGfOHMXHx/u33XjjjcGLLhSOzEckikRoJKb72XIemSpv5lRp0zrZmc/Lrv9Czs33yrRuE+7wAAAAAACQVMciUUpKiv73f//3uO1t2rTRunXr/OvGmCZfJLK7tkuSTAcmrUbjMUkpcn4/QXbpP2UXzpS++lzeHd/KGfU7mbN6hzs8AAAAAADqViR67rnnQh1Ho7CHy6S9u3wr6V3DGwxaHOM4MpdcKdu9l7wvPS3tzZJ36gSZfhfLjLxFxpPw040AAAAAABAidZqTKGJkbZesV2qTJJOYHO5o0EKZjl3kjHtaZsgvJUn24w/kfeQeeT9ZxlxFAAAAAICwaVFFIvvdFt9C+unhDQQtnnHHyLnudjkPPym17ywVF8m+OlXep8bJZu8Od3gAAAAAgBaoRRWJtGOrJMlQJEITYU7vLufPU2RG3iK53dK3G+Sd8Dt5F8yQPVQS7vAAAAAAAC1IiyoS2e+OFomYjwhNh3G55AwfKWfCc1JGH6mqUva9+fKOu0Pepf+UragId4gAAAAAgBagxRSJ7KFi6fs9vpXOFInQ9JiUk+Xc+2c5Y8ZJaR2k4oOys1+W95G75f1shazXG+4QAQAAAAARrE5PN4sI323zvSenyrTmKVJomowxUq++cnqeJ/vJMtm3Xpfyc2VffkZ2yUI5I66Vev1cxmkx9V0AAAAAQCNpMUUiu5NbzdB8mKgomf6XyJ4/UHbpP2UXL5B27ZD3H5OktPYyw0fK9B0o44oOd6gAAAAAgAjRYoYj8GQzNEcmJkbOZdfIeexFmcuukeJbSTl7ZDOnyfunO31zFpWVhjtMAAAAAEAEaDFFIjGSCM2YaZ0g54ob5Ux6WebqUVKbtlJhnm/Ooodulff1F2T3ZIU7TAAAAABAM9YibjezBw9I+bm+lU6nhTcY4ASYuHiZYVfJDh4h++mHsosXSbl7ZZe/I7v8HalrD5mBl8qce4FMNLeiAQAAAADqrkUUiY6OIlJae5n4VuGNBQgCE+2WGTBctt8l0uZ18q54T1q3WtqyUXbLRlnP/8ic10+m70DptG6+CbEBAAAAAAigRRSJ7NaNkrjVDJHHOI7U4xxF9ThHtjBf9qMlsh8tlvYXyK54V3bFu74n+vUd6Hu16xTukAEAAAAATVTLKBJ9vca30OOc8AYChJBpmyxz+W9kL7tG2vSV7OcrZb/8TMrPlX13ruy7c6VTOsr0Ol/m7L5Sl64yTlS4wwYAAAAANBERXySyBXnS7h2SMTJnnRvucICQM1FR0lm9Zc7qLXv4sOy6z2VXr5Q2fCFl75LN3iX73nypdRuZjD4yGX2kbj1l4j3hDh0AAAAAEEaRXyTacGQU0alnyLRuE95ggEZmYmJkzh8gnT9AtqRYdsMX0rrVvveDB2RXLZVdtVQyjpR+uky3DJnuZ0und5eJdoc7fAAAAABAI4r8ItGRW81Mz/PCHAkQXqaVxzeRdd+BspUVvkmu162W/WatlLNH2vGt7I5vZd+bJ0W7fUWj07rLnNbNN/k1RVYAAAAAiGgRXSSyFeXSpnWSKBIBxzKuaKn72b5RQ/Ldlmk3r5M2fS27aZ10oOCHJ6UdPSn1FJlTz5A6nirT+TSpYxduUQMAAACACBLRRSL9Z71UflhKTJY6dgl3NECTZZJSZC4YIl0wRNZa6fs9stv+I23bJLtts7Q3S8rNls3Nlj5b8UPh6KQ0X9GoXUcprYPMKR2ltPYy7phwXg4AAAAAoAEiukj0w61mvWWMCXM0QPNgjPEVfNI6SBcOkSTZkmJp+39kv9sqm7VNytou5edK+3KkfTmya33nWl8DUtJJviepndLB955yspRystQ2RcYV0d0OAAAAADRbEf3bml1/pEiUwa1mwIkwrTxSz94yPXv7t9niImnXDtld26Xs3bI5u6Xs3VLJQV8BKT/XN0G29MPII+NIbZOllFSZ5FQp+WTfcttkqU2S1Kat1Ko1RV0AAAAACIOILRLZ/6yX8r6XXC6p29nhDgeIOMaTUG1eI0m+W9UOHjimaLTL956fK+XlSpUVUsE+qWCfrL754bxjG3a5pIS2UqKvaGTaJEmt20ie1r4CkidBauWRWrWWPAlSbBxFJQAAAAAIgiZbJNq2bZtefvllZWdnKz09XXfffbdSU1PrfL73n69Lksygy2Ri40IVJoBjGGOkhEQpIVHmjLOq7bNer6+AlPe9bH6ur4ibnyubl+ubKPtAgVR8UKqs9BeSpB8VkGpYV5TrmKJRaymulUxsvBQbJ8XF+d5j46TYeF9fcHRfbJwUE+t7kpvbLbljfUUuAAAAAGihmmSRqLi4WI899phuuOEGDRw4ULNmzdKkSZP01FNPyXGcujVyuNT32O6RN4c2WAB1YhzHdztZm7Yyp3Wr8RhbUSEVFUr7C6QDhbIHCn3Fo4NFsiVFviJSycEf3ivKpapKqWi/73W0nVpi+KkS0G5jJHfMD4Wj6CPLMUfej2w3US5fccr1E+9Hl49ZNy6X5Di+l3F+WK62LaqW7TUdHyVFHdlv5HuXObJsjiwffcm/zugrAAAAAD/WJItEq1evlsfj0eDBg2WMUUxMjKKiorRx40addVb10QkVFRWqqKjwrxtjFBcXJ8V7FHXTvTLR7sYOH03Y0V+M+QW5aTJut2+C65ST63S8PVzmKxaVFPvmSCo5KB0qkcpKpbJS2bJD/mWVHZIt/WFZZaW+px+Wl0vWe6RBKx0u870Cfe4JXGOTGqvkLyLph0KSji8oVd93ZLtzzLH+NpzqxanjPq/GIGqIqQHn1XhuHduqsT84flt2tEuVlVX64U+xLjHU8pl1uW40KzkulyorK8MdBpogcgOBkB8IhPxATZxOXaQ/PxWy9ptkkSgrK0tdu3aVMUZr165VRkaGsrOztXPnzuOKRAsXLtS8efP86126dNGTTz6ptqPuVUy3no0dOpqJtLS0cIeAJsJaK1VWyB4+LFt+WPZwme+9/LC8h8tq3K7KStnKiiPvxy5XHFmuOLK98pjlI9srKiVvle/2O2+V5PVWX646drvvvfpy1ZFjjhx3YhcvyTaxylXTxVc0BFLx04eghSI3EAj5gUDID9QoNjakzTfJIlFZWZkSExNVWVmpzZs36/rrr9cHH3ygQ4cOHXfslVdeqREjRvjXj44QOZhyigqysxstZjQPxhilpaUpJyeH+WdQTbXccEdJ7njJE4Y4almuibXWNwrK6/UVfI59HS3+HLft2HX9sN17zH6phvbk+yx7tKj0U20dF21NF/CTh9TpvBpPq+t5NR1Xc5NJyUkqyC9QtZ/TT4bVwPjRrBhJSUlJKigooOaKasgNBEJ+IBDyA7WJ8iSEtP0mWSSKPVIZW7JkiYYOHSpJKiws1GmnnXbcsdHR0YqOjj5uu7WWIgBqRX6gNs0uN4zjm5MIIWWMUewpp8hkZ9eaH9ww1nIdmx8UAXEscgOBkB8IhPxArWqofwRTk/zNIj09XVu3blV5eblSU1NVVlambdu2qUuXLuEODQAAAAAAICI1ySJRnz59tHv3brndblVWVmrWrFlKTU1V9+7dwx0aAAAAAABARGqSRaJWrVpp3Lhx+vjjj3X77bdr9+7dGjt2rBynSYYLAAAAAADQ7DXJOYkk6dRTT9Xjjz8e7jAAAAAAAABaBIbmAAAAAAAAgCIRAAAAAAAAKBIBAAAAAABATXhOohPlckXspSEIyA/UhtxAIOQHAiE/UBtyA4GQHwiE/MCPhTonjLXWhvQTGllFRYWio6PDHQYAAAAAAEBIhKr2EXG3m1VUVGjq1KkqLS0NdyhogkpLS/Xwww+THzgOuYFAyA8EQn6gNuQGAiE/EAj5gdqUlpZq6tSpqqioCEn7EVckkqRVq1YpwgZIIUistdqxYwf5geOQGwiE/EAg5AdqQ24gEPIDgZAfqI21VqtWrQpZ+xFZJAIAAAAAAED9UCQCAAAAAACAosaPHz8+3EEEm+M4OvPMMxUVFRXuUNAEkR+oDbmBQMgPBEJ+oDbkBgIhPxAI+YHahDI3Iu7pZgAAAAAAAKg/bjcDAAAAAAAARSIAAAAAAABQJAIAAAAAAIAoEgEAAAAAAECSK9wB1MW2bdv08ssvKzs7W+np6br77ruVmpoa8JyVK1dq3rx5OnTokM4991yNHj1abrdbkvTWW29p8eLFqqioUP/+/XXDDTfIcaiXNUfBzI0VK1bohRdeUExMjP/YK664QldccUWoLwMh0pD8KC8vV2ZmppYuXapp06YpLS3Nvy9Qv4LmJ5j58c0332jixImKjY31H9uvXz+NHj06pNeA0Khvbuzfv18vvfSSvv32Wx0+fFhpaWm6++67lZ6eLom+I9IEMz/oOyJLfXNj165dmjJligoKCmStVbdu3XTffffJ4/FIou+INMHMD/qOyNOQ76VHPf/888rNzdWxD64/oZqHbeIOHjxob731Vrt06VJbUVFhMzMz7QMPPGCrqqpqPeebb76xN910k92+fbstLS21f/7zn+2LL75orbV2+fLl9s4777Tff/+93b9/v73vvvvswoULG+tyEEShyI1nn322scJHiDUkPw4ePGjHjBlj58+fb3/961/b7Oxs/75AuYPmJ9j5sWHDBvvoo482QuQItYbkRl5env38889tWVmZraiosFOnTrWTJ0+21tJ3RJpg5wd9R+Ro6L8r+fn51lprS0tL7X//93/bt956y1pL3xFpgp0f9B2RpSH5cdSKFSvs6NGjq+XDidY8mvzwmdWrV8vj8Wjw4MFyuVyKiYlRVFSUNm7cWOs5y5Yt0y9+8Qt16dJF1lp1795dK1as0OHDh/Xhhx9q6NChSk1NVUVFhTIyMrR48eJGvCIES7BzA5GlIfnhdrv117/+VVddddVx+8idyBLs/EDkaEhuJCcn6/zzz1dMTIxcLpfS0tL8/9tP3xFZgp0fiBwNyQ2Px6OkpCRZa5Wbm6uSkhKddtppkug7Ik2w8wORpSH5IUl79uzR3Llzdcstt1TbfqI1jyZfJMrKylLXrl1ljNHatWuVkZGhdu3aaefOnQHPOeOMMyRJixcv1uWXX674+Hjt2bOn2r6lS5fq+uuvV35+voqLixvjchBEwc4Nx3H0ySef6LbbbtOECRO0Zs2aRroShEJD8sPtdispKanW9mrLHTQ/wc4PY4y2bNmiW2+9VePGjdOHH34YosgRag3JjWPt3r1by5Yt89+qTN8RWYKdH/QdkaOhubFgwQLdcMMN+q//+i/17t1bPXr08LdH3xE5gp0f9B2RpSH5UV5err///e8aPXq02rZte1x7J1LzaPJForKyMiUmJqqyslKbN29Wjx495DiODh069JPn7N27Vx6PRx6PR47jqLS01L9vw4YN+tnPfub/n5xA7aFpCnZuDBgwQLNmzdL06dM1dOhQTZ8+XV9//XUjXhGCqSH5UZf2asodND/Bzo8ePXpo1qxZ+sc//qGrr75ab775Jl/YmqkTyY2CggI99dRTGjNmjDp16lStPfqOyBDs/KDviBwNzY2rrrpKs2bN0pQpU7Rp0ybNmDGjWnv0HZEh2PlB3xFZGpIfr7zyinr16qVevXrV2l5Dax5Nvkh0dDKuJUuWaOjQoZKkwsJCxcfH13pOXFycJN8wq8GDB6uqqkpFRUWKj49XXFycqqqq9NVXX+ncc89VQUGBJAVsD01TsHPjqPj4eF144YW69NJL6WybsYbkRyB1yR00H8HOj2PbPeecc3Tttddq2bJlJxwnGl9Dc6OoqEiTJk3SrbfeqoyMDP92+o7IEuz8OLZd+o7m7UT+XTHGqH379rrxxhu1fPlySfQdkSbY+XFsu/QdzV9D8uOjjz7SkiVLNGrUKE2aNEmbN2/WqFGj9PHHH59wzaPJF4nS09O1detWlZeXKzU1VWVlZdq2bZu6dOlS6zmdO3fWkiVLlJGRIcdxtHnzZjmOo/bt26tz586aO3euBg0aJEnauHGjTjrpJP9TBNB8BDs3fiw2NlbGmFBeAkKoIfkRSH1yB01fsPPjx+g/mq+G5IbX69XkyZM1YsQInX322dX20XdElmDnx4/RdzRfwfh3JSoqyv+UXfqOyBLs/Pgx+o7mrSH5MWvWLM2YMUOZmZkaO3asunXrpszMTPXr1++Eax5NvkjUp08f7d69W263W5WVlZo1a5ZSU1PVvXt3Sb5q25133umvjknS4MGD9fXXXysuLk7FxcV6/fXXddFFF8ntdmvw4MHatGmTXC6X8vPzNX/+fF1yySXhujycgGDnxvr165WVlSVrrXbt2qXFixerX79+4bo8nKCG5EcggXIHzU+w82Pr1q3aunWrvF6v8vLytGjRIvqPZqohubF+/XoVFBRowIABx7VH3xFZgp0f9B2RoyG5sWnTJm3cuNE/Smj27Nnq37+/JPqOSBPs/KDviCyh+L3lRGoergZdRSNq1aqVxo0bp5deeklz585Venq6xo4dK8f5ob7l9XqrndO9e3fdfffdmj59ukpKStS7d2/dfPPNkqQBAwbowIEDmjhxosrLyzVw4ED98pe/bNRrQnAEOzcKCwv16quvKj8/X0lJSfrNb36j3r17N+o1IXgakh8FBQX6wx/+4F9/+OGHZYzR+PHjA+YOmp9g50dxcbHeeOMN5eTkKCEhQZdccomGDRvWaNeD4GlIbhQWFio/P1+jRo3yb0tJSdFTTz1F3xFhgp0f9B2RoyG5UVVVpZkzZyonJ0cej0f9+/fX1VdfLSnwd1Y0P8HOD/qOyNKQ/AjkRGsexlpr63UFAAAAAAAAiDhN/nYzAAAAAAAAhB5FIgAAAAAAAFAkAgAAAAAAAEUiAAAAAAAAiCIRAAAAAAAARJEIAAAAAAAAokgEAAAAAAAAUSQCAAAAAACAJFe4AwAAAAhkzJgx2rdvnyQpJiZGHTp00GWXXaZ+/fpVOy43N1cLFizQunXrNHToUI0cOTIc4YbNH/7wB0VHR6t///4aPny4XC6+5gEAgPox1lob7iAAAABqM2bMGPXu3VvXXHONioqKtHLlSi1cuFBPPPGETj31VEnShg0b9OSTT6pnz566/PLLdfrpp7e4IklRUZHWrl2r2bNnKykpSePGjVN8fHy4wwIAAM0It5sBAIAmz+VyyePxqF27drruuuuUlpam1atXS/KNIPrb3/6m4cOH66GHHlK3bt1aXIFIkhISEjRo0CBNnjxZlZWVmjZtWrhDAgAAzUzL+wYFAACavVatWungwYOSpBkzZqhz5876zW9+U+Ox69ev1+uvv649e/bIWqu0tDTdf//96tChg3Jzc3Xvvffq2WefVWpqqiT5tz366KM688wzJUllZWV6/fXX9emnn6qkpERt27bVsGHDdPnll/s/Z9WqVVqwYIGys7MVGxurU089VQ899JDcbrck6eOPP9b8+fOVl5en9PR0jR49Wp07d/af/8knn2j+/PnKzc2V4zhq3769xo4dq4SEBH9cL7/8srZs2aKysjIlJyfruuuu04UXXljtej0ej+6//3798Y9/1Jo1a3TeeecF6acOAAAiHUUiAADQrBQXF2v37t266KKLVFBQoDVr1uihhx6S4xw/QNrr9WrKlCm64IIL9Mc//lGStHfv3nrdhuX1evX444+rtLRUDzzwgNLS0pSTk6O8vDz/MUuXLlVmZqZuuukm9e7dW2VlZfryyy919K7+5cuX66WXXtLtt9+unj17asWKFRo/frymT58uj8ej/fv3a/r06Ro5cqQGDx6sqqoqfffdd9VGRGVmZurQoUMaP368EhIStG/fPkVFRdUYc7t27dSnTx8tXryYIhEAAKgzikQAAKBZ8Hq9ysrK0muvvaY2bdqof//+WrNmjRzHUc+ePWs8p7y8XIcOHVJMTIwSExPlcrmUkpJSr8/97LPPtGXLFk2ZMkVpaWmSpKSkJP/+iooKzZo1S7/61a80bNgw//YOHTr4l+fMmaNhw4Zp0KBBkqSRI0fq7bff1tq1azVgwAAVFxerqqpK8fHxSkxMlOM4Oumkk6rFsX//fsXGxioxMVEJCQlKTEwMGHdGRoYyMzPl9XprLKABAAD8GEUiAADQ5C1eoMnLgAAABLdJREFUvFjvvfeeHMfR2Wefrb/85S+KjY1VXl6eWrdu7b+l68diY2N100036Y033tBHH32kIUOGaOjQoUpOTq7zZ2/atEnt2rXzF4h+LCsrSyUlJerdu3eN+w8cOKD8/HwtWbJEH374oX97aWmpfzRShw4ddPnll2vmzJl67733dPHFF+uiiy5S69at/cffcMMNmjZtmu655x5deOGFGjZsmH/i7pqkpKSovLxcBw8eVJs2bep8vQAAoOWiSAQAAJq8AQMG6Ne//rUSEhKq3YLlOI68Xm/Acy+77DL1799fK1eu1LJly7Ro0SKNGTNG/fr1q3GEzaFDh6qtW2sV6GGwR/f9VByjRo1SRkZGtW2tWrXyL994440aPny4VqxYoffff1/z58/Xgw8+qLPOOkuSdOaZZ+q5557Tp59+qhUrVmjs2LEaNmyYbrvttho/72g8jCICAAB1xbcGAADQ5MXFxSkpKem4p5adfPLJOnjwoIqLiwOen5CQoBEjRuiZZ55Rnz59tGjRIn+7ku/x8ZKv4PPmm29WO7dr167Kzs5WTk5OjW137NhRMTEx+vLLL2vc36ZNGyUnJ2vnzp1KTU2t9jq2SCT5Rv9cffXVmjZtmjp27Kh33nmn2n6Xy6X+/fvrkUce0W9/+1stXry41mvPyclRfHy8PB5PwJ8NAADAURSJAABAs3XmmWfK5XJpzZo1Ne7/7LPP9Morr2jz5s0qLCxUVlaW9u3b53+SWatWrZSSkqJ3331X69ev1zPPPHNcIerCCy9U586d9fTTT2vTpk3av3+/tm/frrVr10qSYmJidO2112rRokX64IMPlJeXp5ycHK1cuVLl5eWSpGuvvVZLly7VokWLlJeXp507d2revHn66quvJEnvvfee3nzzTW3ZskX79+/Xjh07VFhYqJNPPlmSlJeXp8mTJ2vt2rXat2+fcnNztX37dnk8nlon4f7iiy+UkZEhY8yJ/6ABAECLwO1mAACg2fJ4PBo4cKAWLFigCy644Li5idq1a6dPP/1UU6dO1f79++V2u9W9e3eNHj3af8wdd9yhF198UV9++aUuvvhiDRkyRKtXr/bvd7lc+stf/qJZs2bp6aefVklJiRITEzVkyBCde+65kqQRI0bI4/HoX//6l1555RXFxsaqQ4cO6tu3ryRp0KBBio6O1rx58zR37lylpKTonHPO8U9unZ6errfeekvLly9XUVGRYmNj1atXL1177bWSfMWstm3baubMmcrLy5PX61WnTp1qfarb5s2btWHDBk2cODG4P3AAABDRjA10kz0AAEATV1xcrAcffFBdu3bV7373u+NGArU0e/fu1cSJE3XeeefVOl8RAABATSgSAQCAZi87O1tPPPGEYmNj9ctf/lIZGRkt7oleO3fu1OrVq/X222/r/PPP1z333MOk1QAAoF4oEgEAgIhQXl6uJUuW6PPPP1ffvn01YsSIcIfUqB599FElJCRoyJAh6tWrV7jDAQAAzRBFIgAAAAAAAPB0MwAAAAAAAFAkAgAAAAAAgCgSAQAAAAAAQBSJAAAAAAAAIIpEAAAAAAAAEEUiAAAAAAAAiCIRAAAAAAAARJEIAAAAAAAAkv4/ktdBAwBEwv4AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x108f1ac88>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import scipy.stats\n",
    "x = np.linspace(0, 0.4, 400)\n",
    "y = scipy.stats.beta.pdf(x, 3, 100)\n",
    "plt.figure(figsize=(14, 3))\n",
    "plt.plot(x, y)\n",
    "plt.xlim(0, 0.4)\n",
    "plt.title('ベータ分布(N=100, success=3)')\n",
    "plt.ylabel('Probability Dencity')\n",
    "plt.xlabel('P(success)')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## statsmodelsで95%信頼区間を求める\n",
    "\n",
    "以降 N=100, Success=3 のケースとする"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'0.8.0'"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import statsmodels.stats.proportion\n",
    "statsmodels.__version__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "confidence interval for a binomial proportion\n",
      "\n",
      "    Parameters\n",
      "    ----------\n",
      "    count : int or array\n",
      "        number of successes\n",
      "    nobs : int\n",
      "        total number of trials\n",
      "    alpha : float in (0, 1)\n",
      "        significance level, default 0.05\n",
      "    method : string in ['normal']\n",
      "        method to use for confidence interval,\n",
      "        currently available methods :\n",
      "\n",
      "         - `normal` : asymptotic normal approximation\n",
      "         - `agresti_coull` : Agresti-Coull interval\n",
      "         - `beta` : Clopper-Pearson interval based on Beta distribution\n",
      "         - `wilson` : Wilson Score interval\n",
      "         - `jeffreys` : Jeffreys Bayesian Interval\n",
      "         - `binom_test` : experimental, inversion of binom_test\n",
      "\n",
      "    Returns\n",
      "    -------\n",
      "    ci_low, ci_upp : float\n",
      "        lower and upper confidence level with coverage (approximately) 1-alpha.\n",
      "        Note: Beta has coverage\n",
      "        coverage is only 1-alpha on average for some other methods.)\n",
      "\n",
      "    Notes\n",
      "    -----\n",
      "    Beta, the Clopper-Pearson interval has coverage at least 1-alpha, but is\n",
      "    in general conservative. Most of the other methods have average coverage\n",
      "    equal to 1-alpha, but will have smaller coverage in some cases.\n",
      "\n",
      "    Method \"binom_test\" directly inverts the binomial test in scipy.stats.\n",
      "    which has discrete steps.\n",
      "\n",
      "    TODO: binom_test intervals raise an exception in small samples if one\n",
      "       interval bound is close to zero or one.\n",
      "\n",
      "    References\n",
      "    ----------\n",
      "    http://en.wikipedia.org/wiki/Binomial_proportion_confidence_interval\n",
      "\n",
      "    Brown, Lawrence D.; Cai, T. Tony; DasGupta, Anirban (2001). \"Interval\n",
      "        Estimation for a Binomial Proportion\",\n",
      "        Statistical Science 16 (2): 101–133. doi:10.1214/ss/1009213286.\n",
      "        TODO: Is this the correct one ?\n",
      "\n",
      "    \n"
     ]
    }
   ],
   "source": [
    "print(statsmodels.stats.proportion.proportion_confint.__doc__)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(-0.0034344809563718326, 0.06343448095637183)"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# 正規分布近似\n",
    "statsmodels.stats.proportion.proportion_confint(3, 100, alpha=0.05, method='normal')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.006229971538306395, 0.0851760529742801)"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Clopper Pearson (EXACT法)\n",
    "statsmodels.stats.proportion.proportion_confint(3, 100, alpha=0.05, method='beta')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(2.168404344971009e-19, 0.003826758485555125)"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Wilson Score\n",
    "statsmodels.stats.proportion.proportion_confint(0, 1000, alpha=0.05, method='wilson')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "正規分布近似は負の区間が含まれるため良くない、Clopper-Pearson法とWilson Scoreではずいぶんと差がある  \n",
    "他の設定で試してもClopper-Pearson法は広くでる"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.0, 0.00034689487247754246)"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "statsmodels.stats.proportion.proportion_confint(0, 369*30, alpha=0.05, method='wilson')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "import scipy.special"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-7.966141830024816"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "scipy.special.logit(0.00034689487247754246)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [],
   "source": [
    "def check_relative_error(cv, imp):\n",
    "    mean = cv/imp\n",
    "    lcb, ucb = statsmodels.stats.proportion.proportion_confint(cv, imp, alpha=0.05, method='wilson')\n",
    "    return (mean - lcb)/mean, (ucb - mean)/mean"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(0.8212378690492708, 3.041500267952385)\n",
      "(0.8232567935859353, 4.448619617870532)\n",
      "(0.8234536293739223, 4.642558597957938)\n",
      "(0.8234732639887764, 4.662688974013383)\n",
      "(0.8234752269621691, 4.664709659040968)\n",
      "(0.8234754232546282, 4.664911804311446)\n"
     ]
    }
   ],
   "source": [
    "print(check_relative_error(1, 10))\n",
    "print(check_relative_error(1, 100))\n",
    "print(check_relative_error(1, 1000))\n",
    "print(check_relative_error(1, 10000))\n",
    "print(check_relative_error(1, 100000))\n",
    "print(check_relative_error(1, 1000000))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}