Ch05-ex05.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## インポート"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "%matplotlib inline\n",
    "\n",
    "np.random.seed(1234)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## データの読み込み"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>PersonID</th>\n",
       "      <th>A</th>\n",
       "      <th>Score</th>\n",
       "      <th>M</th>\n",
       "      <th>Y</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>69</td>\n",
       "      <td>43</td>\n",
       "      <td>38</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>2</td>\n",
       "      <td>1</td>\n",
       "      <td>145</td>\n",
       "      <td>56</td>\n",
       "      <td>40</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>3</td>\n",
       "      <td>0</td>\n",
       "      <td>125</td>\n",
       "      <td>32</td>\n",
       "      <td>24</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>4</td>\n",
       "      <td>1</td>\n",
       "      <td>86</td>\n",
       "      <td>45</td>\n",
       "      <td>33</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>5</td>\n",
       "      <td>1</td>\n",
       "      <td>158</td>\n",
       "      <td>33</td>\n",
       "      <td>23</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   PersonID  A  Score   M   Y\n",
       "0         1  0     69  43  38\n",
       "1         2  1    145  56  40\n",
       "2         3  0    125  32  24\n",
       "3         4  1     86  45  33\n",
       "4         5  1    158  33  23"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df = pd.read_csv('chap05/input/data-attendance-2.txt')\n",
    "df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## モデル式の記述\n",
    "\n",
    "$\\lambda[n] = \\exp(\\beta_1 + \\beta_2 A[n] + \\beta_3 Score[n]) $\n",
    "\n",
    "$M[n]$ ~ $ Poirsson( \\lambda[n] )$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Stanでの実装"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "INFO:pystan:COMPILING THE C++ CODE FOR MODEL anon_model_150bf7d9aaeff585eb332c7a581df661 NOW.\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "Inference for Stan model: anon_model_150bf7d9aaeff585eb332c7a581df661.\n",
       "4 chains, each with iter=2000; warmup=1000; thin=1; \n",
       "post-warmup draws per chain=1000, total post-warmup draws=4000.\n",
       "\n",
       "             mean se_mean     sd   2.5%    25%    50%    75%  97.5%  n_eff   Rhat\n",
       "b1           3.58  2.5e-3    0.1   3.39   3.51   3.58   3.64   3.77   1449    1.0\n",
       "b2           0.26  9.8e-4   0.04   0.18   0.24   0.26   0.29   0.35   1853    1.0\n",
       "b3           0.29  3.9e-3   0.15  -0.01   0.19   0.29   0.39   0.57   1480    1.0\n",
       "lambda[0]   39.62    0.05   1.96  35.94  38.28  39.57  40.86  43.67   1534    1.0\n",
       "lambda[1]   57.53    0.04   1.99  53.71  56.21  57.54  58.83  61.59   2005    1.0\n",
       "lambda[2]   42.91    0.02   1.22  40.61  42.09  42.92  43.71  45.34   2748    1.0\n",
       "lambda[3]   52.85    0.05   2.18   48.7  51.34  52.83  54.27  57.36   1926    1.0\n",
       "lambda[4]   58.63    0.06   2.35  54.05  57.06  58.62  60.15   63.4   1809    1.0\n",
       "lambda[5]   43.41    0.02   1.27  40.99  42.54   43.4  44.28  45.95   2799    1.0\n",
       "lambda[6]   42.06    0.03   1.25  39.67  41.22  42.03  42.88  44.52   2066    1.0\n",
       "lambda[7]    57.7    0.05   2.04  53.77  56.35  57.72  59.03  61.85   1969    1.0\n",
       "lambda[8]   44.24    0.03   1.46  41.45  43.25  44.21  45.24  47.13   2582    1.0\n",
       "lambda[9]   44.17    0.03   1.44  41.42   43.2  44.14  45.15  47.04   2605    1.0\n",
       "lambda[10]   57.2    0.04    1.9  53.55  55.94   57.2  58.42  61.11   2082    1.0\n",
       "lambda[11]  43.66    0.03   1.32  41.14  42.77  43.64  44.56  46.26   2759    1.0\n",
       "lambda[12]  55.33    0.03    1.7  51.98  54.17  55.32  56.47   58.7   2363    1.0\n",
       "lambda[13]  42.06    0.03   1.25  39.67  41.22  42.03  42.88  44.52   2066    1.0\n",
       "lambda[14]  41.94    0.03   1.27   39.5   41.1  41.91  42.77  44.45   2022    1.0\n",
       "lambda[15]  53.84    0.04   1.91  50.12  52.54  53.81  55.14  57.68   2116    1.0\n",
       "lambda[16]  42.24    0.03   1.23  39.88  41.42  42.22  43.06  44.66   2132    1.0\n",
       "lambda[17]  40.35    0.04   1.69  37.17  39.24  40.31  41.42  43.82   1627    1.0\n",
       "lambda[18]  46.84    0.06   2.53  41.94   45.1  46.84  48.51  51.94   1848    1.0\n",
       "lambda[19]  58.46    0.05   2.29  54.02  56.95  58.45  59.91   63.1   1833    1.0\n",
       "lambda[20]  52.85    0.05   2.18   48.7  51.34  52.83  54.27  57.36   1926    1.0\n",
       "lambda[21]  54.54    0.04   1.77  51.06  53.33  54.51  55.73  58.05   2262    1.0\n",
       "lambda[22]  54.38    0.04    1.8  50.88  53.16  54.35  55.61  57.96   2231    1.0\n",
       "lambda[23]  57.36    0.04   1.94  53.65  56.07  57.37  58.64  61.33   2043    1.0\n",
       "lambda[24]  55.25    0.04    1.7   51.9  54.09  55.24   56.4  58.63   2359    1.0\n",
       "lambda[25]  56.13    0.04   1.72  52.74  54.99  56.11  57.26  59.61   2318    1.0\n",
       "lambda[26]  42.54    0.02   1.21  40.21  41.72  42.53  43.35  44.96   2585    1.0\n",
       "lambda[27]  42.54    0.02   1.21  40.21  41.72  42.53  43.35  44.96   2585    1.0\n",
       "lambda[28]  45.22    0.04   1.81  41.76  43.98   45.2  46.42  48.83   2268    1.0\n",
       "lambda[29]  42.18    0.03   1.24  39.82  41.35  42.16   43.0  44.61   2110    1.0\n",
       "lambda[30]  59.23    0.06   2.58  54.22   57.5  59.22  60.89  64.45   1737    1.0\n",
       "lambda[31]  51.96    0.06   2.48  47.21  50.26  51.94  53.56  57.17   1803    1.0\n",
       "lambda[32]  56.95    0.04   1.85  53.41  55.72  56.95  58.15  60.77   2142    1.0\n",
       "lambda[33]  55.81    0.04    1.7  52.46  54.66  55.79  56.95  59.21   2355    1.0\n",
       "lambda[34]  42.97    0.02   1.22  40.64  42.15  42.98  43.78   45.4   2764    1.0\n",
       "lambda[35]  43.29    0.02   1.25   40.9  42.44  43.28  44.14   45.8   2804    1.0\n",
       "lambda[36]  54.31    0.04   1.81  50.76  53.07  54.26  55.53  57.92   2215    1.0\n",
       "lambda[37]   42.0    0.03   1.26  39.58  41.16  41.97  42.82  44.48   2044    1.0\n",
       "lambda[38]  40.76    0.04   1.55  37.82  39.73  40.71  41.76  43.95   1708    1.0\n",
       "lambda[39]   40.7    0.04   1.57  37.73  39.66  40.66  41.71  43.94   1694    1.0\n",
       "lambda[40]  43.16    0.02   1.24  40.82  42.33  43.16  43.98  45.63   2797    1.0\n",
       "lambda[41]  43.41    0.02   1.27  40.99  42.54   43.4  44.28  45.95   2799    1.0\n",
       "lambda[42]  42.73    0.02   1.21  40.42   41.9  42.72  43.53  45.15   2681    1.0\n",
       "lambda[43]  44.37    0.03    1.5   41.5  43.35  44.34   45.4  47.33   2536    1.0\n",
       "lambda[44]  46.63    0.06   2.43  41.96  44.97  46.62  48.23  51.52   1870    1.0\n",
       "lambda[45]   41.7    0.03   1.31  39.17  40.83  41.67  42.56  44.32   1941    1.0\n",
       "lambda[46]  38.63    0.06   2.35  34.23  37.01  38.55  40.12  43.52   1470    1.0\n",
       "lambda[47]   44.5    0.03   1.54  41.54  43.45  44.47  45.54  47.58   2491    1.0\n",
       "lambda[48]  41.34    0.03   1.39  38.67  40.42  41.31  42.25  44.16   1835    1.0\n",
       "lambda[49]  53.84    0.04   1.91  50.12  52.54  53.81  55.14  57.68   2116    1.0\n",
       "M_pred[0]   39.68    0.11   6.59   28.0   35.0   39.0   44.0   54.0   3573    1.0\n",
       "M_pred[1]   57.33    0.12   7.89   42.0   52.0   57.0   62.0   74.0   4000    1.0\n",
       "M_pred[2]   43.13    0.11   6.67   31.0   38.0   43.0   48.0   57.0   4000    1.0\n",
       "M_pred[3]   52.97    0.12    7.5   39.0   48.0   53.0   58.0   68.0   4000    1.0\n",
       "M_pred[4]   58.57    0.13   7.92   44.0   53.0   58.0   64.0   75.0   3613    1.0\n",
       "M_pred[5]   43.34    0.11   6.72   31.0   39.0   43.0   48.0   57.0   3889    1.0\n",
       "M_pred[6]   41.93    0.11   6.63   29.0   37.0   42.0   46.0   55.0   3742    1.0\n",
       "M_pred[7]   57.61    0.13   7.91   42.4   52.0   58.0   63.0   74.0   3967    1.0\n",
       "M_pred[8]   44.33    0.11   6.82   31.0   40.0   44.0   49.0   58.0   3970    1.0\n",
       "M_pred[9]   44.29    0.11   6.79   32.0   40.0   44.0   49.0   58.0   3865    1.0\n",
       "M_pred[10]  57.24    0.13   7.74   43.0   52.0   57.0   62.0   73.0   3772    1.0\n",
       "M_pred[11]  43.67    0.11   6.82   31.0   39.0   44.0   48.0   58.0   4000    1.0\n",
       "M_pred[12]  55.32    0.12   7.63   41.0   50.0   55.0   60.0   71.0   3958    1.0\n",
       "M_pred[13]  41.96    0.11   6.62   30.0   37.0   42.0   46.0   56.0   3769    1.0\n",
       "M_pred[14]  41.92    0.11   6.48   30.0   37.0   42.0   46.0   55.0   3645    1.0\n",
       "M_pred[15]  53.92    0.12   7.64   40.0   49.0   54.0   59.0   69.0   3922    1.0\n",
       "M_pred[16]  42.02     0.1   6.51   30.0   38.0   42.0   46.0   56.0   3960    1.0\n",
       "M_pred[17]  40.48    0.11   6.55   28.0   36.0   40.0   45.0   54.0   3519    1.0\n",
       "M_pred[18]  46.95    0.12   7.27   33.0   42.0   47.0   52.0   61.0   3559    1.0\n",
       "M_pred[19]  58.48    0.13   8.07   43.0   53.0   58.0   64.0   75.0   3713    1.0\n",
       "M_pred[20]  52.74    0.13   7.57   38.0   48.0   53.0   58.0   68.0   3406    1.0\n",
       "M_pred[21]  54.53    0.12   7.52   40.0   50.0   54.0   60.0   70.0   4000    1.0\n",
       "M_pred[22]  54.45    0.12   7.53   40.0   50.0   54.0   59.0   69.0   4000    1.0\n",
       "M_pred[23]  57.19    0.12   7.73   43.0   52.0   57.0   62.0   73.0   4000    1.0\n",
       "M_pred[24]  55.25    0.13   7.69   40.0   50.0   55.0   60.0   71.0   3616    1.0\n",
       "M_pred[25]  55.99    0.12   7.61   42.0   51.0   56.0   61.0   72.0   3862    1.0\n",
       "M_pred[26]  42.53    0.11   6.68   30.0   38.0   42.0   47.0   56.0   3787    1.0\n",
       "M_pred[27]  42.54    0.11   6.63   30.0   38.0   42.0   47.0   56.0   3939    1.0\n",
       "M_pred[28]  45.33    0.12   7.03   32.0   41.0   45.0   50.0   60.0   3732    1.0\n",
       "M_pred[29]  41.95    0.11   6.65   30.0   37.0   42.0   46.0   56.0   4000    1.0\n",
       "M_pred[30]  59.28    0.13   7.94   44.0   54.0   59.0   65.0   75.0   3818    1.0\n",
       "M_pred[31]  52.04    0.13   7.71   37.0   47.0   52.0   57.0   67.0   3632    1.0\n",
       "M_pred[32]  56.84    0.12   7.59   43.0   52.0   57.0   62.0   72.0   3711    1.0\n",
       "M_pred[33]  55.71    0.12   7.69   41.0   50.0   56.0   61.0   72.0   4000    1.0\n",
       "M_pred[34]  43.06    0.11   6.64   31.0   39.0   43.0   47.0  56.59   3770    1.0\n",
       "M_pred[35]  43.15    0.11   6.64   31.0   39.0   43.0   47.0   57.0   3655    1.0\n",
       "M_pred[36]  54.41    0.12   7.53   40.0   49.0   54.0   59.0   69.0   3711    1.0\n",
       "M_pred[37]  41.95    0.11   6.62   30.0   37.0   42.0   46.0   56.0   3871    1.0\n",
       "M_pred[38]  40.72    0.11   6.54   28.0   36.0   41.0   45.0   54.0   3815    1.0\n",
       "M_pred[39]  40.63     0.1   6.49   28.0   36.0   41.0   45.0   54.0   4000    1.0\n",
       "M_pred[40]   43.1     0.1   6.61   31.0   39.0   43.0   48.0   57.0   4000    1.0\n",
       "M_pred[41]  43.34    0.11   6.69   31.0   39.0   43.0   48.0   57.0   4000    1.0\n",
       "M_pred[42]  42.68    0.11   6.58   30.0   38.0   43.0   47.0   56.0   3831    1.0\n",
       "M_pred[43]  44.43    0.11   6.83   31.0   40.0   44.0   49.0   58.0   3879    1.0\n",
       "M_pred[44]   46.6    0.12   7.17   33.0   42.0   46.0   51.0   61.0   3463    1.0\n",
       "M_pred[45]  41.61     0.1   6.48   29.0   37.0   41.5   46.0   55.0   3892    1.0\n",
       "M_pred[46]  38.77    0.12   6.53   27.0   34.0   38.0   43.0   52.0   2862    1.0\n",
       "M_pred[47]  44.49    0.11   6.86   31.0   40.0   44.0   49.0  58.59   3840    1.0\n",
       "M_pred[48]  41.46    0.11   6.56   29.0   37.0   41.0   46.0   54.0   3800    1.0\n",
       "M_pred[49]  53.92    0.12   7.49   40.0   49.0   54.0   59.0   69.0   3732    1.0\n",
       "lp__       6896.5    0.04   1.28 6893.1 6895.9 6896.8 6897.4 6897.9   1289    1.0\n",
       "\n",
       "Samples were drawn using NUTS at Thu May  3 14:09:37 2018.\n",
       "For each parameter, n_eff is a crude measure of effective sample size,\n",
       "and Rhat is the potential scale reduction factor on split chains (at \n",
       "convergence, Rhat=1)."
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model_code=\"\"\"\n",
    "data {\n",
    "  int N;\n",
    "  int<lower=0, upper=1> A[N];\n",
    "  real<lower=0, upper=1> Score[N];\n",
    "  int<lower=0> M[N];\n",
    "}\n",
    "parameters {\n",
    "  real b1;\n",
    "  real b2;\n",
    "  real b3;\n",
    "}\n",
    "transformed parameters {\n",
    "  real lambda[N];\n",
    "  for (n in 1:N)\n",
    "    lambda[n] = exp(b1 + b2*A[n] + b3*Score[n]);\n",
    "}\n",
    "model {\n",
    "  for (n in 1:N)\n",
    "    M[n] ~ poisson(lambda[n]);\n",
    "}\n",
    "generated quantities {\n",
    "  int M_pred[N];\n",
    "  for (n in 1:N)\n",
    "    M_pred[n] = poisson_rng(lambda[n]);\n",
    "}\n",
    "\"\"\"\n",
    "\n",
    "from pystan import StanModel\n",
    "data = dict(N=df.shape[0], A=df.A, Score=df.Score/200, M=df.M)\n",
    "sm = StanModel(model_code=model_code)\n",
    "\n",
    "fit = sm.sampling(data=data, n_jobs=1, seed=123)\n",
    "fit"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 実測値と予測値のプロット\n",
    "* 80% 予測区間"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "ms = fit.extract()\n",
    "\n",
    "qua = np.percentile(ms['M_pred'], q=[20, 50, 80], axis=0)\n",
    "d_qua = pd.DataFrame(qua.T, \n",
    "                     columns=['p'+str(p) for p in [20, 50, 80]])\n",
    "df_pred = df.join(d_qua)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQgAAAEFCAYAAAAfaHkhAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztnXmcXFd157/v1V5dLXVL6kW7vMhXyLYsy/ImeUXYBgIY\nwgSGgZCECQFiGIZkJjPJwGdCJskkJB9DiIfshM1kCAZ/sEkcZLxhMdhYtmVZVutK1mLtvamXql5q\nezV/vKpSVXdt3V2vazvfz0cfVb16795T1VW/c+45995npFIpBEEQCmHW2gBBEOoXEQhBEIoiAiEI\nQlFEIARBKIoIhCAIRXHX2oBSDA6GG67E0tkZZGRkstZmzItGtb1R7Yb6sL2rq90o9ppEEFXG7XbV\n2oR506i2N6rdUP+2i0AIglAUEQhBEIoiAiEIQlFEIARBKIoIhCAIRRGBEAShKC0hELtfOMXuF07V\n2gxBaDiaXiCSlsVoJMrAyBSnByIkLavWJgmC40Qikaq009QCMRKO8vgLpzk1EKF/ZJK9eoDHXzjN\nSDhalfYffPDr3HvvPUSjs9vbs+cn/Pqvf5iPfezXeOSRh6vSnyBUQn//eZ566nEGBvoX3FbTCkTS\nsnj+YD/TsUTe8elYgucP9lclkti9+zF27bqbJ57YnXc8kUjwl395P/ff/wAPPPC3PPLIw1y4MLzg\n/gShErq7e7jiik2EQqEFt9W0AnFuaHKWOGSYjiU4N7yw+e8vvbSXVavW8O53v5fvf/+7AOze/W98\n5zvf4cSJ46xevZYlS5bg8XjYsuUa9u17eUH9CUI5hoeHADAMA6XeRDDYtuA2m1YgJqbjpV+fKv16\nOX74wx/wzne+m3XrNuDxeHjttQPcffdbef/738/ExESeegeDbUxMVGdMKAiF0LqPPXue4fjxo1Vt\nt65Xcy6ENr+n9OuB0q+XYnx8nJ/97KeMjFzgoYe+w8REhO9//ztceeVVdtttbUxOTmTPn5ycqEq4\nJwiF0LqPQ4cOEgy20dPTW9W2m1YgVq4I4j/uLjjM8HvdrFwenHfbu3f/K+94x73cd9+nAZienuaX\nfuldjIyM0NXVzoYNl3D69CnGx8cIBILs2/cyH/jAL8+7P0EoRq447Nx5a1WGFbk0rUC4TJMbN/fw\n/MH8TK7f6+bGzT24zPmPrh599Ad87nN/cLFNv5/bb38zjz76MGvW9PLmN7+dT37yM/zWb30Ky7L4\nhV94F11d3fPuTxAK4bQ4ABj1vO19NTaMSVoW33vmGLG4xR3XrmLl8uCCxKEcXV3tDA6GHWvfSRrV\n9ka1G+Zvu2VZ7N37PGNjYwsWh1IbxjRtBJHBZZq8787La22GIFQV0zTZvv1GotEogUDAuX4ca1kQ\nhKqjdR+nTp0EbJFwUhygBSIIQWgWMjmHtrY2Vq1ajcvl/HZ1EkEIQgNwMSEZZMeOWxdFHKBFBOLJ\nU8/y5Klna22GIMyLXHHYufM2R6oVxWh6gUhaScai4wxODXMmco6klay1SYJQMSdPnqiZOECTC8Ro\ndIynTu/hTOQcg5ND7Bt4ladO72E0OlaV9kut5gR7AtUnPvER3njjRFX6E1qPlStXs2bNupqIAzSx\nQCStJHv79xFN5P94o4koe/v3VSWSKLaaE+DQoYPcd99HOXPmzIL7EVqPiQl7qr7H4+G6666viThA\nEwvE+cmBWeKQIZqI0j85uKD2S63mBIjFYvzxH/8Z69atX1A/QuuhdR9PPrmbwcGBWpvSvGXOyfhU\nydcn4gtb7l1sNWdmZtyWLVsX1L7QmuQmJNvaahM15NK0AhH0lJ5A0uaZ/2Ktcqs5BWE+1LJaUYym\nFYjeYDc+t6/gMMPn9tET7Jp32+VWcwrCXKlHcYAmzkG4TBfbe7bic/vyjvvcPrb3bMVlzn+iyaOP\n/oB77nl79nnuas5MDkIQKiWRSHDmzKm6EwdwcDWnUsoDfB3YACSBjwIJ4GtACjgA3Ke1Lro5ZHVW\ncyZ55Ni/EbPi3LrqJnqCXQsSh3K04srCWtOodsNF26enp7GsZE3EoVarOd8OuLXWO5RSdwF/BHiA\nz2qtn1ZK/TVwL+Dols8u08V7Lv8FJ7sQhHlx5IjG5boU8OD3+2ttTkGcFIjDgFspZQJLgDhwE/BM\n+vXHgLspIRCdnUHc7sWZc15NGjkP0ai2N5rdBw4c4MSJw4yPD3HPPfdgGEWdeE1xUiAi2MOLQ8AK\n4B3AbVrrzLAhDCwt1cDIyMJKkbWgGcLdRqPR7M5NSN5yyy0MDdV2Q+NS4upkkvIzwI+01lcA12Dn\nI7w5r7cDow72Lwh1x8xqRb1vZuykQIwAmUUPF7DzDy8rpe5IH3sbIEsshZbh9deP1GUpsxRODjG+\nCHxVKfUsduTwe8Be4O+UUl6gD3jIwf4Foa7o7e3l3LkzNV1bMVccEwitdQR4X4GXbneqT0GoR+Lx\nGB6Pl1ConVtuub1uE5KFaNqJUoJQD2jdx9NPP5G9kVIjiQOIQAiCY2QSko2MCIQgOEC9rq2YKyIQ\nglBlmkUcQARCEKrK1NQkr79+uCnEAZp4ubcg1IJAwN6W3ufzNbw4gEQQglAVTpw4nt28uLNzWVOI\nA4hACMKC0bqPV155iX37Xqq1KVVHBEIQFkBuQvLqq7fU2pyqIwIhCPOkmaoVxRCBEIR50AriAFLF\nEIR5sXTpUkKhEDfffEvTigOIQAjCnLAsC9M06e1dRXd3L6bZ3EF4c787QagiWvfx3HM/JZFIADS9\nOIAIhCBURCbnMDERIR6P1dqcRUMEQhDKMDMhGQjM/65sjYYIhCCUoFWqFcUQgRCEIoyOjrS0OIBU\nMQShKB0dnWzbdj3Lly9vSXEAiSAEYRbnzp3Fsuw7Qq5du65lxQFEIAQhD637+PnPf0Zf32u1NqUu\nEIEQhDS5CclLLrm01ubUBSIQgoBUK4ohAiG0PCIOxRGBEFoey7JEHIogZU6hZUmlUhiGwaZNm7n8\n8o14PN7yF7UYEkEILYnWffT1vZYVCRGHwohACC1HJudw5swp4vF41dp98tSzPHmquW5YLwIhtBQz\nE5Jer0QOpRCBaAKa0XM5QS2qFY3+txGBEFqC8+fPSilzHkgVQ2gJurt72bhRsWHDJSIOc0AiiCaj\n0UPaXKrxXkZGRgB7e7jNm6+SYcUcEYEQmhat+/jJT57k5MkTtTalYRGBWCSayas0ArkJyRUruiq6\nRv5GsxGBEJoOWVtRPUQgWohae8hC/VfbJhGH6iICITQNlmVx4cKwiEMVcbTMqZT6XeBdgBf4CvAM\n8DUgBRwA7tNaW07aILQOpmlyww03E4tFW2preidxLIJQSt0B7AB2ArcDa4H7gc9qrW8FDOBep/oX\n6genhzYHDhzg7NkzALhcrqw4LOaQKmklGYuOMzg1zFh0HCtl+70jI8c4MnLM0b6dfJ9ORhD3AK8C\nDwNLgP8KfBQ7igB4DLg7/bogzAut+zh58nVM00tv78qa3A5vNDrG3v59nImcyx5zm26292xddFuq\njZMCsQJYD7wDuAR4BDC11qn062FgaakGOjuDuN0uB010hq6u9lnHQhd8RV9bKLltl+qnEhuctK9Q\nH/O191/0E5w9epq2UQ9tbW3s2rWLUChU8fUz2wIItc39b5S0kjx3+HlcPvB6cr+rKQ5FDuH1uDAM\nw9HP3MnvlpMCMQwc0lrHAK2UmsYeZmRoB0ZLNTAyMumgec7Q1dXO4GB41vHIRBSg4GsLIWklOXNh\nkJgVZ5/3MOORKUzDnJcNxWxfKJl+M+T2UcqmUq8dfe04/cfPcOXKN7Fr1y6mplJMTYUrvr5S+8px\nJnKO4bExAGLxJKlUimgyRjKV5NjgaQbCF3C73Ow7fpjeYDcuM9/hVeMzX+h3q5SwOBmP7QHeqpQy\nlFKrgDbgiXRuAuBtgMxKWQCj0TGeOr2HM5FzDE4OsW/gVV4fPc5UYtrRfmtdLj1xwhYHr9/Hzp23\nzYocqkUl73MyPpV9HLcSXJgeIRwLE45FOHThCINTQ4xHx9k38CpPnd7DaHTMEVudwjGB0Fr/EHgZ\n+DnwKHAf8NvA55VSP8OubDzkVP/1RG4C60zkHEkrWZU29/bvI5q46P2slEUkNsGhC4c5NX6mKv3U\nilKf2apVq1na1cmlW1XVSplWyprX3yjoCaSvTzGeTk6mgGgymn08nYhipVJEE1H29u9rqL+Lo2VO\nrfXvFDh8u5N91hszE1j7Bl7F5/axvWcrHb6SKZiSnJ8cyBOHuJXgxPgpwjE7zPzp2edZEVy+4H5q\nQbHP7MolipUdvXi9XtZfdXnV+ptKTHMqfIaElcjrL5qMEXD7S17bG+zG5/YxODmcrVwkLXuoYRom\nBgYWKSbiE7R7Q0QTUfonB1kV6q2a/U4iE6VwLmQu5OWBqniS3NA2lfZemS84QMyK5/XjRBQzFyr1\n0MU+szP9Z3jo+e8zODxYVbuSVjJPHDJEE1FOhc9kf/TFcJkutvdsxcDIHkthYRgG7Z42DMM+HrMu\nbm03EW+c3JrsB+EgM718Lgv1JJnQFiCajM36IntNT7af10ePcypypupRTKUU89CFyoCFPrOh4SGG\nhgbxeDyEUxG6qGzxVSWcnxyYJQ4ZElaCcCxSto0O31J2rr6B05GzWKkkCSvJdGIat3nx55X5ewC0\neRpnEpdEEA6S6+ULsRBPkgltAaxUvjc2DZM2T1v6NYuXB151JIqphFIeem//vlnCNvMzi4yGs+Kw\nbu16kq4U1aTc3yjX85diVVsvIW8bQU+Qdm8or1phYmT/Hj63j55g9QTOaUQgHCTXyxdiIZ4kE9r6\n3D5MI+fLaJgs8S3BTIe2kfgk9sz22WSiGCc5GzlPJDbBRHyS6USUVOqiLVPxac5F+vOGHbmfWWQ0\nzMToeFYcPB5P1b1vub9Rrucvhct0sbZ9NW7TjWEYLPEtwW26MTHwu32YhpGNmmaWOuuZkkMMpZRF\n/rcrDliADxjXWnc6aFvDk/HyhYYZ1fAkHb6l3LnmFsaiYU6Mv8HKtp7sFzSDYUCbu3im38nx8Gh0\njJ+eez6bOAWYMMzssOps5DxTiSmCnmB22LGtaws+t4/J2CTTkUlcbndWHJzwvr3Bbtymu+Aww226\nafdWXkINuP1c3nEJ4ViEmBVn58obeOjwIySxuLZ7Cz3BroYSBygTQWitTa21C/hb4FeAgNY6CLyP\nFilRLoRcL59LNT2Jy3TR6V/KpmVXsCK4PE8cfG4fW7uuzkYThXBqPJxJNqZmBC9WyuJM5BxnIudI\nWIm86CeaiPLS4H6u7d5C0Buks3c5nT3Ls+KQ+cwqSbjOJSm71LeEaDKaF+H43D7Wtq/GNOYWZJuG\nyVLfEroCy1nZ1mMrNJAqEsUtFKeTz5UmKW/UWn8i80Rr/T2l1OeqakmTkvHyGa/ilCcJuP3cvnrH\nrH4ATkfOOhbFFCOTbAx5gpiGmZdrmIxPYQBBT3DWF/DcwDnW+9YU/cwKlUBPRU+yKbQpm3CttLSc\nOW946gJe00s0GSNmxbi841Iu77iEZ878v3m//6nENE+d3sNYdLykDQvBqRJ6LpXK44RS6teUUm1K\nqXal1G9iT6UWKsBlurJeZVWo17Ews1A/1Y5iKvVYmeTfxZzIxa+albLAMFjV1psX8QwNDzEw2M/+\nQ6/keeLMeylWAp1KTOeVcyspLc88zzDsXIHP5eNU5MyCPLOVsjgVPuNoYtjJEnoulUYQHwIeAL6M\nnYP4MfDLVbFAcJxqRTFz8Vi5yT+P6WaZvzNdjk2yzL+MNncAf45o5ZYyr73y2jzhyFBJ2ThFqqLS\ncqm2xqLj/PD4bs5OnC/7PgsRjkWKlk6rNVHKyRJ6LhUJhNb6DeCdSqllWusLC+5VWHQy0QUwry9O\nOY9155pb8gRnZoI246EB1i9ZQ4oU8aRdQoyMholOTeDxeNi4YSPrl60raEM1ysaZc4q1lUqlOBs5\nT0f6syr3PgtRrjRajcSwkyX0XCoaYiiltiqlDgH7lFKrlFKvK6W2VcWCGlPrGYaNQiUeK5diQxu3\n6eb63mu5oXcbPrePibFItpS5ccNGbl57Q9EfYCVl40pLy8XOiyZjJKwEngLlzUrLwuVKo9VIDDtZ\nQs+l0iHGl4H3AN/WWp9VSn0C+GvghqpYUSMWI8mTQZ+0V7a/eW2ZE+uUXI81NGavFr0kx8kW8lgz\nhzZe00O7N5T9bO9ccwtDo0OcHYd3bH0b65etK+mdKy0bV3JOsbaSqSQ+l5c2T5Ch6dnBciWeud0b\nyptFWcyGheB0CT1DpUnKoNa6L/NEa/049lyIhmWxkjzNwnw9Vm7idGlOsjKRSOAyXXR3dnPNNdu4\ndMUlZUP3YlFJwO3PJlwrTcoWPc/lZVVbb9HyZiWe2TRM1ravdry87XQJHSqPIC4opa4hPWlKKfVB\noKFzEYuV5GkWqumxtO7j9OmT7NhxG0DBhGQxCiVcr1p/KReGJ0ueUygpWyjCaUuXZRf6PjNl59eG\nDpFIJR0pby9GCb1SgfgE8HXgSqXUKHAE+GDVrKgBi5XkWQzevPZWx/vIeKy9/fvyjs/VY/WfOEv4\nwgDBYJBUmZWSpWzJTbgW6rvSpGzueRm2dW+Z1/vM/B0yK4NdpouruzaXtWEhLDT5XI5Khxh+rfUt\nwDJgndb6+vTjhmWxkjzNRMZjLTVXEDI7uLZ7C3euuaXifE3/ibP0Hz9T9/etyLzP1aGVdAVXzPl9\nNhPl1mLsBFzA3yul/iP2VvUopdzYScorHLfQIaoRMu9+4RQAd19fP5lHp6MJl+nCb9o/7Ll4rIw4\nZLaJq0QcFiMyKtbfQjzzYtvtJOWGGHdh7wC1EviDnOMJ4G+cMmoxqFbILJQnNhVl8OQ5vH5fVbeJ\nE5ynpEBorX8fQCn1y8A/aa0TSikP4NVaTyyCfY6yWOskAC7zVu8eCQuNXBY78vEGfGy4eiNevw9v\noDbFr1pHe7Xuf75UmoOIYm9AC7AOOKSUaoq7Yi3WOolW5I03ThCP27MKQ51LaiYOwvypVCA+C7wF\nQGt9FLgO+LxTRgmNj9Z97Nv3Iq+88lKtTREWQKVlTq/Wuj/zRGs9oJSqvHgtOMJffm8/AJ9675aK\nzs8kz3afPzWv/na/cArLSjE5YZBMGJweiLByRRDXjNvdad3HoUMHCQaDbN58VcmcQ6GE3sxwvFB4\nXmhmaqOG8fVMpQKxRyn1T8CD6efvA37mjElCvTIVTfDG+TAjMftrs1cP4D/u5sbNPXS228OHXHGo\n51KmUBmVCsR9wKeAj2FvO/cT4CtOGSVcpJpeMdPWfEhaFm+cDxNP5k9umo4leP6gHVwOnT9OO0NZ\ncdjz2gXgQknbnfT6lZYbnSxLLjRqm0sfTlBuHkSv1vo80AP8c/pfhl7gpGOWNRFJy2I0EiUWt4qG\n5fXOuaHJWeKQYTqWIJaw8PqCtLnb2LHj1nTk0NCz8QXKRxB/j3137mew12EYM/6/1FHrmoCRcJTn\nD/az74i9TDgUcM8KyxcbfXIEyPfau184RSjkY8ebugteMzF9cY+DTlf+xKFUKkUsbtG9vJc7t23D\n5XIV7Geu0UKpiKeaZWOhOOXmQbwj/f8li2NO41BJVJC0LJ4/2M90LH93oUxYftf1axomkmjzF97j\nYHh4iOnpaVZ02aKREQehOSg3xPhqqde11h+prjmNQSYqODVg33UpN1mXeyv1c0OTs8Qhw3Qswbnh\nSdZ0OXNn6mqzckUQj8vMG2YMp7eJC/g9BH1S1GpGyg0xnkn//w6gHfgW9jTr9wONdR/zKlEuKrhs\n/cU1bBPTcY6eHS/a1sRUZXdtmtl/JnKZiibwe+fmsS0rxVQ0QdJKFYx8Sg0DOtq9nB6YYCQ8jRWL\nkAglCPg9/Lu7r+fAG+VvUVcr5psDqlbyr5FzUOWGGF8HSO9ifbPW2ko//2fgOefNqz/KRQWnByIE\n3bY3LRaWZ2gLVHbXpgwzI5fRSBTTNBgJRyvKZ2TKlKMRe4FabuRTSb+Do9N4PSbx6QjxaIRVa9t4\n9z07aA+F6lYgSkV7i5EDqnX/C6VSGVtK/vLuHqAxYuMKePPaWyv2FrnJukKEJ2PZx5mwvBB+r5uV\nyytfUl4scrGsFM8f7Cdpld5boVyZ0rIK39glaVl880ea107YFYlEPEoyOobbhI6eywkG63dZfLlo\nr9xn1uj9V4NKBeKPgP1Kqe8qpb4P7MWeft1ylIsK2oPe7GOXabK+t32WSPi9tgeZS5g5M3KxhwpJ\nJqYTHD87xsn+cImry5cpx3OErdR1Hq8ff9syAku6SaZMzg3Pf2OdTOg9MDLF6YEIScvKOzYaiZJI\nlD7n9ECEZLLAMcuqKAfkJLXuvxpUuu39N5VSPwZ2YJc3P661HnDUsjpl5Yog/uPugn94v9fNmu4Q\nFy5cXOhq/3BTdIR8JK0U2zd1s3L53MeguZHLVDTBwMgk0Zi9b+bRc+N8+/EjfPDuK9jQu6Ts9YWI\nxQuLR+a62PQkXl8ADAOP/+LsyPnkUaBw6J06ZL+WORZPWEQmY7QHvbjdZsFz9uoBXjs1yuREbFYY\nv6LDX9KG+dpeKeU+c6f7rwaVbnvvBX4NuBd4Avh4+ljL4TJNbtzcg9+br63ZqKDAkMIwDAI+N6GA\nhzVdoYrFIdcrRibjpFIpkpbFwMjkrCFBImnxg2ePE0sU9lh+n4vpWIKJqThT0STReJLxiRipVIpU\nKkU0nuDc0ASjkWhe221+DxPhEUaHzzMRHpnVblvAU9B7xxIJRsJRBkenePHQQJ5dhULvVCrF0bNj\nHD07RiqVwrIshsemmJhOMDQ2hWWlZp2TuU6/MZJ3DGwP/frp/GOFbHeSauegakGlU63/DzAIbMOe\nan058A+06N21Ott93HX9GsYnY8Ti1ryjglLM9LBBv4tT/fbjmeJgGgZej8lULMGrRy9wneqe1dZr\nx0YYj8SYjCaIxZMYhh3duF0mbpcJBkzGksRiSTwuM5tECw+fZDo8jMvtxR9sz2vX73UT8Lp4/IXT\ned57an+CM4MTDI7a+34+/uIpnt1/lraAh4DPw77DQ0xG43k3FQ5PxkmkhzKxmMX4VIzpqC0g8QQM\njEzSEfJmz0laKfxet31dwiKRtLLHMng9JrGEhc8zu9Iz1xzQfCgXbTrdfzWo9Bt9ndb694C41noS\n+07f1zpnVv3jMk06Qj66OwNFo4KkZZciw5MxRsJRwpOxrIctRSEPaxoGPZ0BBsem8u4TbRgGS9o8\n2Z2hR8LRgm1F44m861IpGJ+MMzoRzeYYUqkU07EEo5EoP/r5SV577TUOH+5jRVuSru7VuNwXPZ5p\nGmzf1M0Lhwbz7ExaFvuPDjMwOpX13omkxdnhSV49Nsy54Qn2vT7EsbPjWQEAezhhGwFjE1GSsyKk\nlD1mTx9OJlP51+Ucy/3MLlu1tHi053CpsWy02QClzkojiFR6SJH5C6zIeSwUIBMBDI9PMx1NYKXs\nO8HvefUcHSFfyTJXseSW3+dm5bIgkckYhmmSSrfpzhnWzGwz01Z4Mo5hQMDnssN1wOM28LhMDMMe\nL49OxIjH7bxG37HzHD06zlWr27hi86W4Pf5sxJRMWvi9LqajiVl2DlyYyhEcgBTjE3GsVArLShGZ\nitHTGSCRtDg9ZOdqDMPA47bfQzSeBGO25zJNW8CiiSQ+jwuXyxbEzHVA9lguvcuDXHvFCkejvVIs\nRrTpJJVa+SXsG/b2KqW+hF3F+KJjVjU4Ga89GY1nxQHsH8zpwQhT0XjJMlep5Fbv8iAetxuv28Q0\njbx7SgS8bq6+LH+z8UxbGU9rGAamaeAyDdymiZWyf3jjE/GLQ5eUfYdqCy/m0stxe/yYppGNmAI+\nN4ZhFLRzKnbxhkOp9Hu2cvIA8USK9qAHt8skkbRsQYDsMcOwbXPlvDf7mInX48LAFsTMsKE96MHt\nNvOOZciE8ZVEe05S6/4XQqWWPgZ8HLvceQx4p9a65DTsVibjtSOTcWZOL0gkLcJT8ZJlrlLJLZdp\notZ1zCqdBrxu7r31Erzu/KAw01aup82QEQq7gpFjqAFtbSHWrFlFMmXmlUAzMzEjU/Fs4jTfDldu\nM7PCTI/b/uGv6QrhdpnZYYGRHg6s72nHMOxzfF4Xpmn/bxh2tLS+dwmXrVqaJx5qfWfeMWisML6e\nqXSI8azW+k3Awbk0rpTqBl7E3h07AXwN+ztzALgvMzOz2ZjptQGC/osfdTxdUixW5ipbSu0KsWp5\nGy/qAeJJi7u2r+Xqy5bNEofctlKpFKZh5HnzNr+HFCkmowlcpklseoqklcLr8+N2mXSE7DJhpgSa\nmYnpcpm4XPDGQJhT/RF6ll1MtnUvC3ByIELSSmHMUAjTNAgFvOn34eKyVUu4EI7i87izoXfSsvji\nd/YTjSdxuQy8bjvRmEymuKR3Ce+5/RJcpsn3njmWDdm3qB4Gh8J5xxopjK9nKv0EX1FKfVjZrMv8\nK3VBevfrvwEyt7C6H/is1vpWbOfSFJveFqKU1wbweOzjxcpc5ZJbpmngdpt0tPvo6ghwneouKA65\nbQV8HjpCvmzlwDRgTVcbl6/uYEPPEuKxaaLRaRLxGIYBa7pCGIaRLYGevzDJoZMjxBIXhxCZxGn/\nhclsJOEyTbZctpzujkA6ErDP87hMujuDmGbOkMjnYfWKUF7o7XW72bS+k46QD7/XjWma+L1uOkI+\n7r5hHV63e3bI7jIbOoyvZyqNIG7EvpN3bhao3H4Qf459c53fTT+/jouLvx4D7gYertjSBiLjta1U\nCtMgb5jhdpm0Bzxly1zVTG7lthWNJTlxbhyXy+CWa1axcnmQQ319HOgbw+P10b5kGQG/105CxpL0\nX7BnUmbmTZiGnYtwp8XP73OztifESIFI4E+/9TLxpMVbrltDeCrOofQ+knBR7F44NHu+XcDn5oq1\nHdn37vWYLAl6G2LtQrNRbrn3KuABYALYA/x3rfVoqWvS1/0qMKi1/pFSKiMQhtY681MJY6/vKEln\nZxC3u373F/jg2zcXPN7bs5S33+rl6RdP4/W6mUonKn0eFxvXdbJ8aYDbrl3N8qWlb/8HsCY9M/La\nzSuzx0IdzaRzAAANwElEQVQh+4fi89l/vtwl5pW0NZVODF67eSUHDhzg7NljrF+exFiyGisdVAYC\nXs6NjLFhdQcn+8eJJ63sJLDwdJzuziBtbRd/sCnDZOWKtjw7e1bYMy7fe9cmkkmLbzzWRzSW5O6b\nN7Cm2/b8fafHZr2HzPtbsiR/JmShczLHurraZx3LpdRri0Gp/mtlUyWUiyD+ETuH8LfYS7zvByrZ\nA+Ij2KXRtwBbgW8AubN32oGyQjMyUv9z1WfS1dXO4GCYkXCUiYkopFL43CYBnwuf28XGVUtQ6zqw\nYgkGB0uvnwCIpFde5p6bORZNzyOopJ1C173ySh8vvLCXYLCNDZcpOjqXcnYgTCxusaLdy0SHn1TS\nnjxlWSmSmUlKSQhPpN9fmmTCIhKJ5tky0z6vaeD1uwm6jex09FLvbyaFzhkcDGc/80JtFTq/FhTr\nP2N7LSklUOUEYrXW+h4ApdQTwL4y5wOgtb4t81gp9TR2BeTPlFJ3aK2fBt4GPFVJW43IxclJyWwe\nYmna257sj6DWdSyo/cxeDZkt3eaKWtcJQG/vKi699HIuu+xy9rx2IVvKBAgFPHkzHX0eu6KQKYXm\nTkrye934W3LifeU06lb85Qa02fqW1jqe+3we/DbweaXUzwAv8NAC2qpr6n0V3/Sk7bFM0+Tqq68p\nuDX9zFKrYRgsbfNlk4yZSUm5iVOh+ag0SZlhzrMntdZ35Dy9fa7XNyKLtYovEwnMhYGzxxg4e4y1\nl14NFPdqhUqtm9Z1YKVSHDk1xqoVbXVXTizlpRvVg9eacgJxpVLqWM7z1ennBpDSWsuu1gWo11V8\nA2ePMXDmKF5fgEBb4WXhGTLl0cw9LzIEfR42re8k4HM3zH6awvwpJxBXLIoVTUY1V/FVyytq3Uc7\nQ/RsXEU0cKm9t0MZipVan3jxTMX9lqPQe6hkm3yJCBaHcntSvrFYhjQTxbxvZhn1YofkF2+H18bO\nnbem73hVGZkJSIBEDC1IfQwem5CM9+0I+WgPelnbHeKKtR2LPtnHsiwGBwey4iD3yhTmwlyTlMIc\ncJkmgfRkpvfdeXlNbDBNk5tu2kk8HiMQsIc2xUL2rq52HvzXOS23ybt+JjPvOj7XYYEMI2qPRBBN\nyuHDh+jvPw+A2+3OioMgzAWJIBxmPqXIhZLJObS3t9PV1Y05h5xHJV5bPHvrIBFEk5GbkLzppp1z\nEgdBmIlEEE3EzGpFtRKSEjG0LuJemoQTJ445Ig5CayMRRJOwcuUqzp8/x5YtW0UchKohAuEwTofn\n09PT+P1+fD4/N92009G+hNZDhhgNjNZ9PPnk44yOzm/ZtyCUQwSiQckkJD0eD16vbMYgOIMIRAPi\nVLVCEGYiAtFgiDgIi4kIRAMRi8U4ceK4iIOwaEgVo4Hwer3s3HkbpmmIOAiLgkQQDcDRo0eYmLB3\ngQ6FQiIOwqIhAlHnaN3HgQP7efnlvbU2RWhBRCDqmNyE5LZt22ttjtCCiEDUKVKtEOoBEYg65MgR\nLeIg1AVSxahDVqzoYunSpdxww80iDkJNEYGoIxKJBG63m87OZdx++y4MQ+5WJdQWGWLUCVr38eyz\nTxON2jd5FXEQ6gERiDogk5BMJBIkk8lamyMIWUQgaszsaoXsPi3UDyIQNURKmUK9IwJRIyKRMIcP\nHxJxEOoaqWLUiFConRtv3CFrK4S6RiKIRebUqTdIJOy7fnd394g4CHWNCMQionUfL720l/3799Xa\nFEGoCBGIReJiQjLIpk1vqrU5glARIhCLQK447Nx5mwwrhIZBBMJhRByERkaqGA7j9/ullCk0LCIQ\nDpFKpTAMg/XrL2H16rW43fJRC42HI99apZQH+CqwAfABfwgcBL4GpIADwH1aa8uJ/muN1n2Ew2G2\nbduOaZoiDkLD4lQO4kPAsNb6VuCtwAPA/cBn08cM4F6H+q4pBw4c4NChg4yMDBOLxWptjiAsCKcE\n4rvA59KPDSABXAc8kz72GPAWh/quGVr3sX///mxC0u/319okQVgQjsS+WusIgFKqHXgI+Czw51rr\nVPqUMLC0XDudnUHcbpcTJladAwcOcPLk67S1tbFr1y5CoVCtTZoXXV3ttTZhXjSq3VDftjs2OFZK\nrQUeBr6itf62UuoLOS+3A6Pl2hgZmXTKvKoyNDTIc8/tJRgMsmvXLqamUkxNhWtt1pzp6mpncFDs\nXkzqwfZSAuXIEEMp1QPsBv6b1vqr6cMvK6XuSD9+G/CsE33XghUrurjqqi3s3Hlbw0YOglAIpyKI\n3wM6gc8ppTK5iE8DX1ZKeYE+7KFHQzM4OMCKFV0YhsFll22stTmCUHWcykF8GlsQZnK7E/3VgswM\nyc2br2LjRlVrcwTBEWSq9TzInT69evWaWpsjCI4hAjFHZG2F0EqIQMwBEQeh1RCBqJBUKsXU1KSI\ng9BSyCKBCjEMg2uu2UYsFsXnkxmSQmsgEUQZtO7j2LHXAVskRByEVkIEogSZnMPRo69nN5oVhFZC\nBKII+QnJW2XJttCSiEAUQKoVgmAjAjGD06dPijgIQhqJm2ewcuVqNmwYYuNGJeIgtDwSQaQJh8cB\ncLlcXHPNNhEHQUAEArBzDk899WPOnz9ba1MEoa5oeYHIJCQDgQBLlpTd5EoQWoqWFgipVghCaVpW\nIEQcBKE8LSkQyWSS8+fPiTgIQhlasszpcrnYseMW4vEEwWCw1uYIQt3SUhHEkSOaoaFBADwer4iD\nIJShZQRC6z4OHjzA/v37sKymvOOfIFSdlhCI3ITkTTftwDRb4m0LwoJp+l+KVCsEYf40tUAcO3ZU\nxEEQFkBTVzFWrlzJ+fNn2bpV1lYIwnxoSoGIRqfx+fwEAkF27Li11uYIQsPSdEOMzMKrzOpMQRDm\nT1MJRCYh6XK5cLlctTZHEBqephEIqVYIQvVpCoEQcRAEZ2h4gYhGpzl+/KiIgyA4QMNXMXw+Pzt2\n2NvSizgIQnVpeIEAZCcoQXCIhh9iCILgHCIQgiAURQRCEISiiEAIglAUEQhBEIoiAiEIQlFEIARB\nKIoIhCAIRTFSqVStbRAEoU6RCEIQhKKIQAiCUBQRCEEQiiICIQhCUUQgBEEoigiEIAhFEYEQBKEo\nTbFhTK1QSnmArwIbAB/wh8BB4GtACjgA3Ke1rsu7BSuluoEXgbuABI1j9+8C7wK8wFeAZ2gA29Pf\nl69jf1+SwEep889dIoiF8SFgWGt9K/BW4AHgfuCz6WMGcG8N7StK+sv6N8BU+lCj2H0HsAPYCdwO\nrKVBbAfeDri11juAPwD+iDq3XQRiYXwX+Fz6sYHtDa7D9mgAjwFvqYFdlfDnwF8DZ9PPG8Xue4BX\ngYeBR4Ef0ji2HwbcSikTWALEqXPbRSAWgNY6orUOK6XagYeAzwKG1jozfz0M1N2GmUqpXwUGtdY/\nyjlc93anWQFsB34J+DjwIGA2iO0R7OHFIeDvgC9T55+7CMQCUUqtBZ4Cvqm1/jaQO35sB0ZrYlhp\nPgLcpZR6GtgKfAPoznm9Xu0GGAZ+pLWOaa01ME3+j6qebf8Mtu1XANdg5yO8Oa/Xne0iEAtAKdUD\n7Ab+m9b6q+nDL6fHyQBvA56thW2l0FrfprW+XWt9B7AP+DDwWL3bnWYP8FallKGUWgW0AU80iO0j\nwFj68QXAQ51/X2Q15wJQSv0F8H7skDHDp7FDRy/QB3xUa52sgXkVkY4iPo4d+fwdDWC3UuoLwJ3Y\nDu73gOM0gO1KqRB21Wsltq1/Aeyljm0XgRAEoSgyxBAEoSgiEIIgFEUEQhCEoohACIJQFBEIQRCK\nIou1WpR0ye1PsacuTwDjwO9rrZ9Ilz5/X2v9dO0szEcpldJaG7W2o9WQCKIFUUoZ2OsYYsBmrfU1\nwH8CvpkzaUcQZB5EK5IWga8Cl+WsA0Ap9ZvAewEX9iKuTdiL0D6jtX5aKbUL+AL20uQR4ANa6yGl\n1IeB/4ztcF7EXrI8rZQaTD/vxZ7M9KDW+qF0X3uB38COXP4KWA5MAp/SWr+slNoAfAsIAc8BH5MI\nYvGRCKI1uR7YmysOaX6Sfg0gorXeBvwKdmThw16M9nGt9XbsCGSbUupK7H0NdmittwIDwH9Jt7EC\n+JP08W8A/x5AKbURCGitX8Jej/A76b5+A/i/6WsfAL6Wvvan1X37QqWIQLQmKQrnn3IXDv0DgNZ6\nPzCIHU08AjyslHoA6NNa78ae8rwReE4ptQ97P4NNOe08n/7/X4Cb0itfPwA8mM6DXA/8Y/rabwMh\npdRy4A7gO+lrH8ReGi0sMiIQrcnzwPb0pjG53Ay8kH6cyDluAHGt9Rexf7ivA19QSv0P7OHIP2ut\nt6a9/Q3AJzMXaq2n0v/HsPdueBfwPuwfvQuYzlybvv5G7IVMKS5+P1Pkr5IVFgkRiBZEa/0s8Brw\npYxIKKWuwx5C/K/0aR9MH9+OvbnJEaXU80C71vpLwBeBbcDTwHuUUt3p5OdfYecjCvFN4LeBC1rr\nN7TWY+l2P5Tu6y7sYQ7Aj7F37AL4Rewt/YRFRgSidflFIAocUEodxF5Z+KGc0mZIKfUy9q5T/0Fr\nHcdeOfk1pdSL2PmC/6m1fgX4PPAktuiYwJ8U6lBr/VPsvRu+lXP4g8CvK6X2A/8beH86N/JJ4L3p\n42/H3kxFWGSkiiEIQlEkghAEoSgiEIIgFEUEQhCEoohACIJQFBEIQRCKIgIhCEJRRCAEQSjK/wfh\nNTtthkngNQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0xf43f048>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# 実測値と予測値の散布図\n",
    "fig = plt.figure()\n",
    "ax = fig.add_subplot(1,1,1)\n",
    "for a, df_sub in df_pred.groupby(by='A'):\n",
    "    ax.errorbar(x=df_sub.M, y=df_sub.p50, \n",
    "                yerr=[df_sub.p50-df_sub.p20, df_sub.p80-df_sub.p50],\n",
    "                fmt='o', ms=8, alpha=0.5, marker='o', label='A:'+str(a))\n",
    "ax.legend()\n",
    "ax.set_aspect('equal')\n",
    "ax.set_xlabel('Observed')\n",
    "ax.set_ylabel('Predicted')\n",
    "ax.plot([5, 90], [5, 90], ls=\"--\", c=\".3\", alpha=.5)\n",
    "\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}