Chap7/Chap7.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 7. 一般化線形混合モデル(GLMM) -個体差のモデリング-"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* 「GLMではうまく説明できない」現象をうまく扱えるようにGLMを強化する。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Q. 一般化線形モデル(GLM)とはなんだったか？"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A.\n",
    "* 確率分布、リンク関数、線形予測子を組み合わせて、応答変数$y_i$と説明変数$x_i$を関連づけるシステム\n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* しかし、単純なGLMでは現実のデータ解析に応用できない。\n",
    "* 理由は、実験・調査で得られたカウントデータのばらつきは、ポアソン分布や二項分布**だけ**ではうまく説明できないから。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* 例えば、植物個体の種子数の調査。説明変数がどの個体でも同じ値だとしても、現実の生物では「説明変数以外は全部均質」といった条件は満たされないため、平均$\\lambda$の値が全個体共通とはならないだろう。\n",
    "* 「説明変数が同じなら平均も同じ」というGLMの仮定が成立しないので、集団全体の種子数の分布はポアソン分布で期待されるよりも大きくばらつく"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* この話の要点は「データにばらつきをもたらす『個体間の差異』なるものを、定量化できない」ということ。\n",
    "* 人間は自然の全てを測定できるわけではない。\n",
    "* しかし「何か原因不明な個体差がある」ことは統計モデルとして表現できる。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Q. 一般化線形混合モデル(generalized linear mixed model,GLMM)とは？"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A.\n",
    "* 「人間が測定できない・測定しなかった個体差」を組み込んだGLM\n",
    "* 複数の確率分布を部品とする統計モデル。\n",
    "    * 例えば、データのばらつきは、二項分布・ポアソン分布、個体のばらつきは正規分布など。\n",
    "* 多くの実際のデータ解析ではこれを用いるのが良い\n",
    "\n",
    "--- "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Q. 個体差、場所差とは？"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A.\n",
    "* 観測者がデータ化していないけれど、個体や調査地に由来する原因不明の差異"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 7.1. 例題：GLMでは説明できないカウントデータ"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* 架空植物の各個体から8個の種子を取ってきて、そのうち幾つが生存しているかを調べたとする。\n",
    "* データ解析の目的：生存種子数が葉数とともにどのように増大するのかを解明すること"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/takizawa/dev/python/books/midori/lib/python3.6/site-packages/statsmodels/compat/pandas.py:56: FutureWarning: The pandas.core.datetools module is deprecated and will be removed in a future version. Please use the pandas.tseries module instead.\n",
      "  from pandas.core import datetools\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import statsmodels.api as sm\n",
    "import statsmodels.formula.api as smf\n",
    "%matplotlib inline\n",
    "\n",
    "d = pd.read_csv(\"data7.csv\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style>\n",
       "    .dataframe thead tr:only-child th {\n",
       "        text-align: right;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: left;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>N</th>\n",
       "      <th>y</th>\n",
       "      <th>x</th>\n",
       "      <th>id</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>8</td>\n",
       "      <td>0</td>\n",
       "      <td>2</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>8</td>\n",
       "      <td>1</td>\n",
       "      <td>2</td>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>8</td>\n",
       "      <td>2</td>\n",
       "      <td>2</td>\n",
       "      <td>3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>8</td>\n",
       "      <td>4</td>\n",
       "      <td>2</td>\n",
       "      <td>4</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>8</td>\n",
       "      <td>1</td>\n",
       "      <td>2</td>\n",
       "      <td>5</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   N  y  x  id\n",
       "0  8  0  2   1\n",
       "1  8  1  2   2\n",
       "2  8  2  2   3\n",
       "3  8  4  2   4\n",
       "4  8  1  2   5"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "d.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style>\n",
       "    .dataframe thead tr:only-child th {\n",
       "        text-align: right;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: left;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>N</th>\n",
       "      <th>y</th>\n",
       "      <th>x</th>\n",
       "      <th>id</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>count</th>\n",
       "      <td>100.0</td>\n",
       "      <td>100.000000</td>\n",
       "      <td>100.000000</td>\n",
       "      <td>100.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>mean</th>\n",
       "      <td>8.0</td>\n",
       "      <td>3.810000</td>\n",
       "      <td>4.000000</td>\n",
       "      <td>50.500000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>std</th>\n",
       "      <td>0.0</td>\n",
       "      <td>3.070534</td>\n",
       "      <td>1.421338</td>\n",
       "      <td>29.011492</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>min</th>\n",
       "      <td>8.0</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>2.000000</td>\n",
       "      <td>1.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>25%</th>\n",
       "      <td>8.0</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>3.000000</td>\n",
       "      <td>25.750000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>50%</th>\n",
       "      <td>8.0</td>\n",
       "      <td>3.000000</td>\n",
       "      <td>4.000000</td>\n",
       "      <td>50.500000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>75%</th>\n",
       "      <td>8.0</td>\n",
       "      <td>7.000000</td>\n",
       "      <td>5.000000</td>\n",
       "      <td>75.250000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>max</th>\n",
       "      <td>8.0</td>\n",
       "      <td>8.000000</td>\n",
       "      <td>6.000000</td>\n",
       "      <td>100.000000</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "           N           y           x          id\n",
       "count  100.0  100.000000  100.000000  100.000000\n",
       "mean     8.0    3.810000    4.000000   50.500000\n",
       "std      0.0    3.070534    1.421338   29.011492\n",
       "min      8.0    0.000000    2.000000    1.000000\n",
       "25%      8.0    1.000000    3.000000   25.750000\n",
       "50%      8.0    3.000000    4.000000   50.500000\n",
       "75%      8.0    7.000000    5.000000   75.250000\n",
       "max      8.0    8.000000    6.000000  100.000000"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "d.describe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.PathCollection at 0x10a29b780>"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAW4AAAD8CAYAAABXe05zAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAFu9JREFUeJzt3W9sJHd9x/HP15tFbMIfQ7NCueTIBQkZFaKQdJWUgqIm\nKRgEhCgFKZGoVJ7cE0qhrYwwigqpQHngCkHVCunEfyUcfxJjQQQYJPKgPMjBOiYyycUSTUMue6Ex\nUgz5syKbvW8feNdn+3Y9M/aMZ34z75cU5Tw7M/f7fnf2c76Zn+9n7i4AQDgm8h4AACAZghsAAkNw\nA0BgCG4ACAzBDQCBIbgBIDAENwAEhuAGgMAQ3AAQmPOyOOmFF17oR44cyeLUAFBKS0tLv3f3Zpx9\nMwnuI0eOqN1uZ3FqACglM/tt3H25VQIAgSG4ASAwBDcABIbgBoDAENwAEJhYwW1m/2RmD5nZr83s\nuJm9NOuBAQBGi5wOaGYXS/pHSX/u7l0z+46kWyR9LeOxYRcLyx3NLa7q9HpXhyYbmpme0k1XXrzv\nfUOSZg9C6FEa9S4sd3T7Dx7S08/3NvedbNT16RvfWLh6k9pZ83VvaOreB5/Uenej1ledX9en3vtG\nSTqnB5I0YdIZly7e0a8iXhcWtXTZILjvl3SFpD9KWpD0H+7+k3HHtFotZx53dhaWO5qdX1G319/c\n1qjXdMfNl59zUSXZNyRp9iCEHqVR79/+xcX69i9Pqdc/9zNfnzDNfeCKwtSb1KiaR5kwyczUP7N7\n7g37dc9S58CuCzNbcvdWnH0jb5W4e0fSv0t6XNKTkv6wW2gje3OLq+dcoN1eX3OLq/vaNyRp9iCE\nHqVR7/ETo0NbknpnvFD1JjWq5lHOuCJDWzrbr6JeF5HBbWavkvQ+SZdJOiTpAjP74Ij9jppZ28za\na2tr6Y8Um06vd2NvT7JvSNLsQQg9SqPefsTfrotUb1JZjH1cv4rQpzgPJ/9G0v+6+5q79yTNS/qr\nnTu5+zF3b7l7q9mM9eP22KNDk43Y25PsG5I0exBCj9Kot2a2p98jBFmMfVy/itCnOMH9uKS/NLPz\nzcwk3SDpZLbDwm5mpqfUqNe2bWvUa5qZntrXviFJswch9CiNem+95rDqtdFhVJ+wQtWb1KiaR5kw\nqTax+x9g0tl+FfW6iJxV4u4nzOxuSQ9IelHSsqRjWQ8M4w0fjMR52p1k35Ck2YMQepRWva1LX13K\nWSWjak5jVknr0lcX8rqInFWyF8wqAYBkUp1VAgAoFoIbAAJDcANAYAhuAAgMwQ0AgSG4ASAwBDcA\nBIbgBoDAENwAEBiCGwACQ3ADQGAIbgAIDMENAIEhuAEgMAQ3AASG4AaAwESugGNmU5K+vWXT6yT9\nq7t/Ps2BLCx3Ul1pIu3zoTii3tudr1/3hqbue2StEtfCqN5I41fOGe7fWe/KTBquqzJcLaasfUrL\n1v7VzNR337aCTlYSrYBjZjVJHUnXuPtvx+2XdAWcheWOZudX1O31N7c16jXdcfPleyo+7fOhOKLe\n21Gv71TWa2FU7fWaSS71zpz9nA/rl7Rrr+o109z7ryhdn9Ky27W2l2ssyxVwbpD0P7uF9l7MLa6e\nU3y319fc4mohzofiiHpvR72+U1mvhVG19/q+LbSls/VH9arX91L2KS279S/rayzyVskOt0g6PuoF\nMzsq6agkvfa1r0100tPr3UTbD/p8KI6o9zbue1zGayFJTVXuU1qiepNl72J/x21mL5F0o6Tvjnrd\n3Y+5e8vdW81mM9EgDk02Em0/6POhOKLe27jvcRmvhSQ1HZpsxNq/jH1KS1Rvsuxdklsl75L0gLv/\nX9qDmJmeUqNe27atUa9tPljJ+3wojqj3dtTrO5X1WhhVe71mqk/Ytm3D+qN6Va9ZKfuUlt36l/U1\nluRWya0ac5tkv4Y38NOaBZL2+VAcUe/tqNerMqtkXG9GbdtaP7NK9mZrvws5q8TMLpD0uKTXufsf\novZPOqsEAKouyaySWN9xu/tzkv5sX6MCAKSCn5wEgMAQ3AAQGIIbAAJDcANAYAhuAAgMwQ0AgSG4\nASAwBDcABIbgBoDAENwAEBiCGwACQ3ADQGAIbgAIDMENAIEhuAEgMLGC28wmzexuM3vEzE6a2Vuy\nHhgAYLS4S5d9QdKP3f39g0WDz89wTJlaWO5UZkmzKtUqJau3ar0ZhR6E24PI4DazV0q6VtLfS5K7\nvyDphWyHlY2F5Y5m51fU7fUlSZ31rmbnVyQpiDcriSrVKiWrt2q9GYUehN2DOLdKLpO0JumrZrZs\nZl8arEEZnLnF1c03aajb62tucTWnEWWnSrVKyeqtWm9GoQdh9yBOcJ8n6SpJX3T3KyU9J+kTO3cy\ns6Nm1jaz9traWsrDTMfp9W6i7SGrUq1Ssnqr1ptR6EHYPYgT3E9IesLdTwy+vlsbQb6Nux9z95a7\nt5rNZppjTM2hyUai7SGrUq1Ssnqr1ptR6EHYPYgMbnf/naRTZjY12HSDpIczHVVGZqan1KjXtm1r\n1GuamZ4ac0S4qlSrlKzeqvVmFHoQdg/izir5iKS7BjNKHpX0oeyGlJ3hA4cQnyInVaVapWT1Vq03\no9CDsHtg7p76SVutlrfb7dTPCwBlZWZL7t6Ksy8/OQkAgSG4ASAwBDcABIbgBoDAENwAEBiCGwAC\nQ3ADQGAIbgAIDMENAIEhuAEgMAQ3AASG4AaAwBDcABAYghsAAkNwA0BgCG4ACEysFXDM7DFJz0jq\nS3ox7j/2nbaF5U6Qq1VkIUkvqtq3qtadBD1Kpij9irt0mSRd5+6/z2wkERaWO5qdX1G315ckdda7\nmp1fkaTKXWhJelHVvlW17iToUTJF6lcwt0rmFlc3GzbU7fU1t7ia04jyk6QXVe1bVetOgh4lU6R+\nxQ1ul/QTM1sys6OjdjCzo2bWNrP22tpaeiMcOL3eTbS9zJL0oqp9q2rdSdCjZIrUr7jB/TZ3v0rS\nuyR92Myu3bmDux9z95a7t5rNZqqDlKRDk41E28ssSS+q2req1p0EPUqmSP2KFdzu3hn8/ylJ35N0\ndZaDGmVmekqNem3btka9ppnpqYMeSu6S9KKqfatq3UnQo2SK1K/Ih5NmdoGkCXd/ZvDrd0j6t8xH\ntsPw5n8RnujmLUkvqtq3qtadBD1Kpkj9MnfffQez12nju2xpI+i/6e6f3e2YVqvl7XY7nRECQAWY\n2VLcqdaR33G7+6OSrtj3qAAAqQhmOiAAYAPBDQCBIbgBIDAENwAEhuAGgMAQ3AAQGIIbAAJDcANA\nYAhuAAgMwQ0AgSG4ASAwBDcABIbgBoDAENwAEBiCGwACEzu4zaxmZstmdm+WAwIA7C5yIYUtPirp\npKRXZDSWTC0sdwqx5FDRlb1PZa9vp93qrVovyiRWcJvZJZLeLemzkv450xFlYGG5o9n5FXV7fUlS\nZ72r2fkVSeJC3aLsfSp7fTvtVq+kSvWibOLeKvm8pI9LOpPhWDIzt7i6eYEOdXt9zS2u5jSiYip7\nn8pe30671Vu1XpRNZHCb2XskPeXuSxH7HTWztpm119bWUhtgGk6vdxNtr6qy96ns9e20W71V60XZ\nxPmO+62SbjSzxyR9S9L1Znbnzp3c/Zi7t9y91Ww2Ux7m/hyabCTaXlVl71PZ69tpt3qr1ouyiQxu\nd59190vc/YikWyT9zN0/mPnIUjQzPaVGvbZtW6Ne08z0VE4jKqay96ns9e20W71V60XZJJlVEqzh\nwxaeoO+u7H0qe307xam3Kr0oG3P31E/aarW83W6nfl4AKCszW3L3Vpx9+clJAAgMwQ0AgSG4ASAw\nBDcABIbgBoDAENwAEBiCGwACQ3ADQGAIbgAIDMENAIEhuAEgMAQ3AASG4AaAwBDcABAYghsAAkNw\nA0Bg4iwW/FIz+4WZPWhmD5nZ7QcxMADAaHGWLvuTpOvd/Vkzq0v6uZn9yN3vz3hskqSF5c6ellfa\n63FlUOXape31T55fl7v0h26vcr1I4zooy7V028KKjp84pb67ama69ZrD+sxNl+c9rD2LDG7fWNvs\n2cGX9cF/6a93NsLCckez8yvq9vqSpM56V7PzK5K068Wz1+PKoMq1S+fW//Tzvc3XqtSLNK6DslxL\nty2s6M77H9/8uu+++XWo4R3rHreZ1czsV5KekvRTdz+R7bA2zC2ubl40Q91eX3OLq5kcVwZVrl0a\nXf9WVelFGtdBWa6l4ydOJdoegljB7e59d3+zpEskXW1mb9q5j5kdNbO2mbXX1tZSGdzp9W6i7fs9\nrgyqXLsUr84q9CKN66As11J/zILo47aHINGsEndfl3SfpHeOeO2Yu7fcvdVsNlMZ3KHJRqLt+z2u\nDKpcuxSvzir0Io3roCzXUs0s0fYQxJlV0jSzycGvG5LeLumRrAcmSTPTU2rUa9u2Neo1zUxPZXJc\nGVS5dml0/VtVpRdpXAdluZZuveZwou0hiDOr5CJJXzezmjaC/jvufm+2w9owfACS9Kn2Xo8rgyrX\nLp1bf1VnlaRxHZTlWho+gCzTrBLzDO7ztFotb7fbqZ8XAMrKzJbcvRVnX35yEgACQ3ADQGAIbgAI\nDMENAIEhuAEgMAQ3AASG4AaAwBDcABAYghsAAkNwA0BgCG4ACAzBDQCBIbgBIDAENwAEhuAGgMAQ\n3AAQmMgVcMzssKRvSHqNJJd0zN2/kPXAhhaWO8GvwHFQbltYOWeVj9alrw62f3t97xeWO7r9Bw/p\n6ed7m9smG3V9+sY3BlN7muL2kc9atKL0KHIFHDO7SNJF7v6Amb1c0pKkm9z94XHHpLUCzsJyR7Pz\nK+r2+pvbGvWa7rj5ci6oHW5bWNGd9z9+zvbahKl/5ux7HEr/9vreLyx3NHP3g+r1z72u6xOmuQ9c\nUfja0xS3j3zWomXdo1RXwHH3J939gcGvn5F0UtKBvJNzi6vbmiRJ3V5fc4urB/HbB+X4iVMjt28N\nbSmc/u31vZ9bXB0Z2pLUO+NB1J6muH3ksxatSD1KdI/bzI5IulLSiRGvHTWztpm119bWUhnc6fVu\nou1V1k+wdmgI/dvre7/f18smbh/5rEUrUo9iB7eZvUzSPZI+5u5/3Pm6ux9z95a7t5rNZiqDOzTZ\nSLS9ympmsfcNoX97fe/3+3rZxO0jn7VoRepRrOA2s7o2Qvsud5/PdkhnzUxPqVGvbdvWqNc0Mz11\nUEMIxq3XHB65vTaxPdBD6d9e3/uZ6SnVa6P/EKtPWBC1pyluH/msRStSj+LMKjFJX5Z00t0/l/2Q\nzhre8C/CU9yi+8xNl0tSaWaV7PW9H77OrJINcfvIZy1akXoUZ1bJ2yT9t6QVSWcGmz/p7j8cd0xa\ns0oAoCqSzCqJ/I7b3X8uKf4NVABApvjJSQAIDMENAIEhuAEgMAQ3AASG4AaAwBDcABAYghsAAkNw\nA0BgCG4ACAzBDQCBIbgBIDAENwAEhuAGgMAQ3AAQGIIbAAITGdxm9hUze8rMfn0QAwIA7C5yIQVJ\nX5P0n5K+ke1QtltY7mwuEfTKRl29/hk990JfUrxlqLYeX4VlmJLWW6b+lKmWtNCTvdnZt+ve0NR9\nj6wVro+RS5dJkpkdkXSvu78pzkn3u3TZwnJHs/Mr6vb6Y/epT5jmPnDFyCaOOr5Rr+mOmy8vRNPT\nlrTeMvWnTLWkhZ7sTZzcybKPSZYuK+Q97rnF1V2bJ0m9M665xdXYx3d7/bH7hy5pvWXqT5lqSQs9\n2Zs4uVOUPqYW3GZ21MzaZtZeW1vb17lOr3f3tV/S7aFLqw8h9qdMtaSFnuzNfnPnIKUW3O5+zN1b\n7t5qNpv7Otehyca+9ku6PXRp9SHE/pSplrTQk73Zb+4cpELeKpmZnlKjXtt1n/qEaWZ6KvbxjXpt\n7P6hS1pvmfpTplrSQk/2Jk7uFKWPkbNKzOy4pL+WdKGZPSHpU+7+5SwHNbzxv9dZJTuPL9LT4Cwk\nrbdM/SlTLWmhJ3szqm9BzypJar+zSgCgaoKfVQIAGI/gBoDAENwAEBiCGwACQ3ADQGAIbgAIDMEN\nAIEhuAEgMAQ3AASG4AaAwBDcABAYghsAAkNwA0BgCG4ACAzBDQCBIbgBIDCxgtvM3mlmq2b2GzP7\nRNaDAgCMF2fpspqk/5L0dklPSPqlmX3f3R/OenDjLCx3NLe4qs6Y1ZZfdX5dn3rv9qXNhscUbQmi\nvdpaz+T5df2p19fzvTPn7DdqmbfbFlZ0/MQp9d1VM9Ot1xzWZ266/CCHn7mF5Y4+/f2HtN7tbW4b\ndV1UVdk+D0lF1V/0/kQGt6SrJf3G3R+VJDP7lqT3ScoluBeWO5qdX1G31x+7z9PP9zRz94OSNtaR\n23lMZ72r2fmVzddDs7Oep5/vjd13vdvTzHfP9uK2hRXdef/jm6/33Te/Lkt4Lyx3NPPdB9U7s31Z\nvp3XRVWV7fOQVFT9IfQnzq2SiyWd2vL1E4NtuZhbXN01tId6fdfc4urYY7q9/ubroYnbg6HembO9\nOH7i1Mh9xm0P0dzi6jmhPbT1uqiqsn0ekoqqP4T+pPZw0syOmlnbzNpra2tpnfYcp8fcHtlt33HH\nJDlXkexl3MNj+mMWhx63PURR/Qn1fU9L2T4PSUXVH0J/4gR3R9LhLV9fMti2jbsfc/eWu7eazWZa\n4zvHoclG4n3HHZPkXEWyl3EPj6mZjXx93PYQRfUn1Pc9LWX7PCQVVX8I/YkT3L+U9Hozu8zMXiLp\nFknfz3ZY481MT6lRr0XuV6+ZZqanxh7TqNc2Xw9N3B4M1SfO9uLWaw6P3Gfc9hDNTE+pPjH6D6Kt\n10VVle3zkFRU/SH0J/LhpLu/aGb/IGlRUk3SV9z9ocxHNsbw4UCSWSVbjynqU+IkdtaTZFbJ8AFk\nmWeVDGtlVsloZfs8JBVVfwj9Mc/g3mar1fJ2u536eQGgrMxsyd1bcfblJycBIDAENwAEhuAGgMAQ\n3AAQGIIbAAKTyawSM1uT9NvUT7y7CyX9/oB/zyKjH9vRj+3ox3ZF6Mel7h7rpxczCe48mFk77lSa\nKqAf29GP7ejHdqH1g1slABAYghsAAlOm4D6W9wAKhn5sRz+2ox/bBdWP0tzjBoCqKNN33ABQCcEH\nt5kdNrP7zOxhM3vIzD6a95jyZGYvNbNfmNmDg37cnveYisDMama2bGb35j2WvJnZY2a2Yma/MrPK\n/2twZjZpZneb2SNmdtLM3pL3mKLEWXOy6F6U9C/u/oCZvVzSkpn9NM/FjHP2J0nXu/uzZlaX9HMz\n+5G735/3wHL2UUknJb0i74EUxHXunve85aL4gqQfu/v7B2sOnJ/3gKIE/x23uz/p7g8Mfv2MNj6c\nxfmHcw+Yb3h28GV98F+lH2SY2SWS3i3pS3mPBcViZq+UdK2kL0uSu7/g7uv5jipa8MG9lZkdkXSl\npBP5jiRfg9sCv5L0lKSfunul+yHp85I+LunclSaqySX9xMyWzOxo3oPJ2WWS1iR9dXAr7UtmdkHe\ng4pSmuA2s5dJukfSx9z9j3mPJ0/u3nf3N2tjfdCrzexNeY8pL2b2HklPuftS3mMpkLe5+1WS3iXp\nw2Z2bd4DytF5kq6S9EV3v1LSc5I+ke+QopUiuAf3cu+RdJe7z+c9nqIY/JXvPknvzHssOXqrpBvN\n7DFJ35J0vZndme+Q8uXuncH/n5L0PUlX5zuiXD0h6Yktfyu9WxtBXmjBB7eZmTbuT51098/lPZ68\nmVnTzCYHv25IerukR/IdVX7cfdbdL3H3I9pY6Ppn7v7BnIeVGzO7YPAQX4NbAu+Q9Ot8R5Ufd/+d\npFNmNlwJ+AZJhZ/YUIZZJW+V9HeSVgb3dSXpk+7+wxzHlKeLJH3dzGra+IP5O+5e+Slw2PQaSd/b\n+H5H50n6prv/ON8h5e4jku4azCh5VNKHch5PJH5yEgACE/ytEgCoGoIbAAJDcANAYAhuAAgMwQ0A\ngSG4ASAwBDcABIbgBoDA/D9IXNyv3/CeeQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10a0d0978>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# そのまま散布図にすると点が重なってしまい数の多さがわからないため、\n",
    "# 正規分布のズレを加えている。\n",
    "d[\"x_bure\"] = np.random.normal(loc=0,scale=0.1,size=100)\n",
    "plt.scatter(d.x+d.x_bure,d.y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* 第6章と同じように、GLMを使って、データから種子の生存確率を推定してみる。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### GLMのあてはめ"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* 「$N$個のうち$y$個が生存していた」といった構造のカウントデータは二項分布を用いる。\n",
    "* 二項分布を使ったGLMの一つであるロジスティック回帰の統計モデルを用いる。\n",
    "* ロジスティック回帰ではリンク関数はロジットリンク関数。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import statsmodels.api as sm\n",
    "import statsmodels.formula.api as smf\n",
    "model = smf.glm('y + I(N-y)~ x',data=d,family=sm.families.Binomial())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table class=\"simpletable\">\n",
       "<caption>Generalized Linear Model Regression Results</caption>\n",
       "<tr>\n",
       "  <th>Dep. Variable:</th>  <td>['y', 'I(N - y)']</td> <th>  No. Observations:  </th>  <td>   100</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Model:</th>                 <td>GLM</td>        <th>  Df Residuals:      </th>  <td>    98</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Model Family:</th>       <td>Binomial</td>      <th>  Df Model:          </th>  <td>     1</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Link Function:</th>        <td>logit</td>       <th>  Scale:             </th>    <td>1.0</td>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Method:</th>               <td>IRLS</td>        <th>  Log-Likelihood:    </th> <td> -322.80</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Date:</th>           <td>Thu, 06 Jul 2017</td>  <th>  Deviance:          </th> <td>  513.84</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Time:</th>               <td>01:11:08</td>      <th>  Pearson chi2:      </th>  <td>  53.5</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>No. Iterations:</th>         <td>4</td>         <th>                     </th>     <td> </td>   \n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "      <td></td>         <th>coef</th>     <th>std err</th>      <th>z</th>      <th>P>|z|</th>  <th>[0.025</th>    <th>0.975]</th>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Intercept</th> <td>   -2.1487</td> <td>    0.237</td> <td>   -9.057</td> <td> 0.000</td> <td>   -2.614</td> <td>   -1.684</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>x</th>         <td>    0.5104</td> <td>    0.056</td> <td>    9.179</td> <td> 0.000</td> <td>    0.401</td> <td>    0.619</td>\n",
       "</tr>\n",
       "</table>"
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.summary.Summary'>\n",
       "\"\"\"\n",
       "                 Generalized Linear Model Regression Results                  \n",
       "==============================================================================\n",
       "Dep. Variable:      ['y', 'I(N - y)']   No. Observations:                  100\n",
       "Model:                            GLM   Df Residuals:                       98\n",
       "Model Family:                Binomial   Df Model:                            1\n",
       "Link Function:                  logit   Scale:                             1.0\n",
       "Method:                          IRLS   Log-Likelihood:                -322.80\n",
       "Date:                Thu, 06 Jul 2017   Deviance:                       513.84\n",
       "Time:                        01:11:08   Pearson chi2:                     53.5\n",
       "No. Iterations:                     4                                         \n",
       "==============================================================================\n",
       "                 coef    std err          z      P>|z|      [0.025      0.975]\n",
       "------------------------------------------------------------------------------\n",
       "Intercept     -2.1487      0.237     -9.057      0.000      -2.614      -1.684\n",
       "x              0.5104      0.056      9.179      0.000       0.401       0.619\n",
       "==============================================================================\n",
       "\"\"\""
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "result = model.fit()\n",
    "result.summary()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x10a409550>]"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAW4AAAD8CAYAAABXe05zAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3X2UFNWZP/DvM22rjRpHcdwwCEHdX8agLho7vhIUfBkM\niqiYiKIbd89hVzfRvOxsJMeNya5ZN5nsmuz5JTlhk2zCi6DC0CIKIxuNRrMQB0czKkxiFIFGl1EZ\nXhtoep79407PG93TVT3VU3Wrvp9zOM5UV7f3PlX9pai6VVdUFUREZI8qvxtARETuMLiJiCzD4CYi\nsgyDm4jIMgxuIiLLMLiJiCzD4CYisgyDm4jIMgxuIiLLHFGJDz3ppJN03LhxlfhoIqJQWr9+/fuq\nWuNk3YoE97hx49DS0lKJjyYiCiURecfpujxVQkRkGQY3EZFlGNxERJZhcBMRWYbBTURkGUfBLSJf\nFpHXReQ1EVksIkdXumFERFRYyeGAIjIawN0AxqtqRkQeBXAzgF9UuG00iFRrGo3N7djWmUFtdQIN\n9XWYce7oIa9rEy9rYEONvOhvqjWNbz3xOnbsy/asW52I45vTzwxcf90a2OfJZ9Rg5avvojNj+nrC\niDjuv/ZMADisBgBQJUCXAqMH1CuI+4WUmrqsO7jXApgAYBeAFID/UNWni70nmUwqx3FXTqo1jblN\nbchkcz3LEvEYHrzh7MN2Kjfr2sTLGthQIy/6e+N5o/HIS1uQzR3+nY9XCRpvmhCY/rpVqM+FVAkg\nIsh1DZ57+XotW58etv1CRNaratLJuiVPlahqGsD3AGwG8C6AnYOFNlVeY3P7YTtoJptDY3P7kNa1\niZc1sKFGXvR38brCoQ0A2S4NVH/dKtTnQroUJUMb6K1XUPeLksEtIicAuA7AqQBqARwjIrMLrDdH\nRFpEpKWjo8P7llKPbZ0Zx8vdrGsTL2tgQ4286G+uxL+ug9RftyrR9mL1CkKdnFycvALA26raoapZ\nAE0ALh64kqrOU9WkqiZrahzdbk9lqq1OOF7uZl2beFkDG2rkRX9jImX9P2xQibYXq1cQ6uQkuDcD\nuFBERoiIALgcwIbKNosG01Bfh0Q81m9ZIh5DQ33dkNa1iZc1sKFGXvR31gVjEI8VDqN4lQSqv24V\n6nMhVQLEqgb/CwzorVdQ94uSo0pUdZ2ILAXwMoBDAFoBzKt0w6i4/IURJ1e73axrEy9rYEONvOpv\n8mMnhnJUSaE+ezGqJPmxEwO5X5QcVVIOjiohInLH01ElREQULAxuIiLLMLiJiCzD4CYisgyDm4jI\nMgxuIiLLMLiJiCzD4CYisgyDm4jIMgxuIiLLMLiJiCzD4CYisgyDm4jIMgxuIiLLMLiJiCzD4CYi\nskzJGXBEpA7AI30WnQbgG6r6fS8bkmpNezrThNefR8FRatsOfH3yGTV4dmNHJPaFQrUBis+ck18/\n3ZmBCJCfVyU/W0xY6+SVvvWLiSCn2m8GnUpxNQOOiMQApAFcoKrvFFvP7Qw4qdY05ja1IZPN9SxL\nxGN48Iazy+q8159HwVFq2xZ6faCw7guF+h6PCaBAtqv3e57vP4BBaxWPCRpnTghdnbwy2L5Wzj5W\nyRlwLgfwp8FCuxyNze2HdT6TzaGxuT0Qn0fBUWrbFnp9oLDuC4X6ns1pv9AGevtfqlbZnIayTl4Z\nrH6V3sdKnioZ4GYAiwu9ICJzAMwBgLFjx7r60G2dGVfLh/vzKDhKbVun2ziM+4KbPkW5Tl4pVZtK\n1s7xEbeIHAlgOoDHCr2uqvNUNamqyZqaGleNqK1OuFo+3J9HwVFq2zrdxmHcF9z0qbY64Wj9MNbJ\nK6VqU8nauTlVcjWAl1X1f71uREN9HRLxWL9liXis58KK359HwVFq2xZ6faCw7guF+h6PCeJV0m9Z\nvv+lahWPSSjr5JXB6lfpfczNqZJZKHKaZKjyJ/C9GgXi9edRcJTatoVej8qokmK1KbSsb/85qqQ8\nfesdyFElInIMgM0ATlPVnaXWdzuqhIgo6tyMKnF0xK2qewGMHFKriIjIE7xzkojIMgxuIiLLMLiJ\niCzD4CYisgyDm4jIMgxuIiLLMLiJiCzD4CYisgyDm4jIMgxuIiLLMLiJiCzD4CYisgyDm4jIMgxu\nIiLLMLiJiCzjKLhFpFpElorIRhHZICIXVbphRERUmNOpy34AYLWqzuyeNHhEBdtUUanWdGSmNItS\nXwF3/Y1abQphDeytQcngFpHjAUwC8HkAUNWDAA5WtlmVkWpNY25TGzLZHAAg3ZnB3KY2ALBiY7kR\npb4C7vobtdoUwhrYXQMnp0pOBdAB4L9EpFVEfto9B6V1GpvbezZSXiabQ2Nzu08tqpwo9RVw19+o\n1aYQ1sDuGjgJ7iMAfBLAj1X1XAB7Adw7cCURmSMiLSLS0tHR4XEzvbGtM+Nquc2i1FfAXX+jVptC\nWAO7a+AkuLcC2Kqq67p/XwoT5P2o6jxVTapqsqamxss2eqa2OuFquc2i1FfAXX+jVptCWAO7a1Ay\nuFX1PQBbRKSue9HlAN6oaKsqpKG+Dol4rN+yRDyGhvq6Iu+wV5T6Crjrb9RqUwhrYHcNnI4q+SKA\nRd0jSt4CcEflmlQ5+QsONl5FditKfQXc9TdqtSmENbC7BqKqnn9oMpnUlpYWzz+XiCisRGS9qiad\nrMs7J4mILMPgJiKyDIObiMgyDG4iCrcPPwQWLABmzgRWr/a7NZ5wOqqEiMge77wDPP44kEoBzz8P\n5HJAbS0wY4bfLfMEg5uI7KcKtLUBy5ebsH7lFbN8/Hjga18zgX3eeUBVOE4yMLiJyE5dXcDatSas\nm5qAt94CRICLLgK++13guuuAj3/c71ZWBIObiOyRzQK//rUJ6lQKeO89IB4HpkwxR9bTpwMf/ajf\nraw4BjcRBdv+/cCaNcCyZcCKFcCOHcCIEcDVVwM33AB85jNAdbXfrRxWDG4iCp69e4FVq0xYr1wJ\n7NkDHH+8OaK+8UbgqquARPAfBlUpDG4iCobdu4EnnwSWLgWeegrIZICaGmDWLHNkPWUKcOSRfrcy\nEBjcROSfnTvNEfVjj5kx1gcOmHPUf/VXZtz1xInAEYypgVgRIhpeu3aZc9WPPgo0NwMHDwKjRwN/\n8zcmrC++GIjFSn9OhDG4iajydu0CnniiN6wPHABOOQW46y7gppuACy8MzRjr4cDgJqLK2LOnN6xX\nrTJhPXo0cOedDOshYnATkXf27TMXFh95xFxozGTMreZ/+7cmrC+6iGHtAUfBLSKbAOwGkANwyOnD\nvr2Wak1bOVtFJbipRVTrFtV+u+FJjQ4cMBcWH3nEnLveuxc4+WRzgfGznzUXGEMS1kHZp9wccU9W\n1fcr1pISUq1pzG1qQyabAwCkOzOY29QGAJH7MrqpRVTrFtV+uzGkGh06BDzzDLBkibmLcedOYORI\n4NZbgc99Drj00tBdYAzSPmXNX4ONze09BcvLZHNobG73qUX+cVOLqNYtqv12w3WNurrMk/buugsY\nNQqorzc3yFx/vTnifvdd4Cc/MeOtQxbaQLD2KadH3ArgaRFRAD9R1XkDVxCROQDmAMDYsWO9a2G3\nbZ0ZV8vDzE0tolq3qPbbDUc1UgXWrwcWLzanQtJpc8fi9OnAzTcDU6cCRx89TC32V5D2KafBPVFV\n0yJyMoA1IrJRVZ/vu0J3mM8DzGTBHrcTtdUJpAsUqLY6ere9uqlFVOsW1X67MWiNNm40Yb14MfDH\nP5oHOU2dCjQ2AtdeCxx7rA8t9leQ9ilHp0pUNd393+0AlgM4v5KNKqShvg6JeP9/fiXiMTTU1w13\nU3znphZRrVtU++3GwBp9dNf7+LuW5XjyF3cDn/gE8M//DIwZA/znf5qn8K1YYW4/j2BoA8Hap0oe\ncYvIMQCqVHV3989XAfinirdsgPzJ/yBc0fWbm1pEtW5R7bcbM84djfjOHXjtP36OS1vW4Pytr6NK\nFfjUp4CHHjIjQmpr/W5mYARpnxLVwc9qiMhpMEfZgAn6h1X124O9J5lMaktLizctJCJv7dtnbox5\n+GFzY0w2C9TVAbfcYv78+Z/73cJIEpH1TodalzziVtW3AEwYcquIyD+5nBm+t3ChGb63Z485mr77\nbhPW555rZo8hK/DOSaKwUgVeftmE9ZIl5jz1Rz5ixlnfeiswaVIoh+1FAYObKGzefhtYtMj82bjR\nPMN62jQT1tOmRWb4XpgxuInC4IMPzDOtFy4EXnzRLJs0CfjKV8yjUk84wd/2kacY3ES22r/fTEKw\ncKF5sFM2C4wfDzz4oDlvXYEb4SgYGNxENunqAl54AViwwBxh79xpZoz54heB2bOBc87hRcYIYHAT\n2WDjRhPWixYB77wDHHOMmYdx9mzg8st5kTFiGNxEQdXRYUaDzJ8PtLSYR6NeeSXw7W8DM2aY8KZI\nYnATBcn+/ebW8gULzBP3Dh0ypz/+7d/M7eajRvndQgoABjeR37q6zEiQ+fN7z1vX1poRIbfdBpx1\nlt8tpIBhcBP55Y9/NEfWCxYAmzaZUx833mjCevJknremohjcRMPpww/N5Lnz5wP/8z9mBMgVV5gn\n8V1/Pc9bkyMMbqJKO3jQnK+eP9883OngQeDMM4HvfMfczTiaTywkdxjcRJWQnzlm/nwzGcH77wM1\nNcCddwJ/+Zccb01DwuAm8tLWrWas9fz5wBtvAEcdZab5uv12M0djPO53CykEGNxEQ7V3L7B8OfDL\nXwK/+pU52r7kEjNx7k038Tkh5DnHwS0iMQAtANKqek3lmkRkga4u4LnnzJH10qXm+dbjxgH/+I9m\nVAgnI6AKcnPEfQ+ADQA+UqG2VFSqNR2IKYeCLux1GnL//vAHE9YLFgCbNwPHHWeeb3377cDEiebu\nxgAZrL9h39Zh5ii4ReQUANMAfBvAVyraogpItaYxt6kNmWwOAJDuzGBuUxsAcEftI+x1Krt/H34I\nPPKICey1a004X3UV8K//Clx3HTBixHA037XB+gsg1Ns67JweHnwfwD8A6KpgWyqmsbm9ZwfNy2Rz\naGxu96lFwRT2OrnqXzZrhu7NnGluM7/rLnM6pLHRXIBctcrcgh7Q0AYG72/Yt3XYOZnl/RoA21V1\nvYhcNsh6cwDMAYCxAXsO8LbOjKvlURX2OpXsnyrQ2mqOrB9+2DzkqabGhPbtt1s3hK+c7RmWbR12\nTo64LwEwXUQ2AVgCYIqILBy4kqrOU9WkqiZramo8bubQ1FYnXC2PqrDXqVg/zq7aB3zve8Bf/AVw\n3nnAj38MXHqpOeJOp4GHHrJyMt3BtmfYt3XYlQxuVZ2rqqeo6jgANwN4RlVnV7xlHmqor0Mi3v+5\nD4l4DA31dT61KJjCXqe+/Ts6ux/T33gO85d+E49/52agoQE49lgT2u+9Zx72dM01Vo+7Hmx7hn1b\nh10kxnHnL7bwCvrgwl6nGRNGYeT6tej8yc9w2e9/jeMOZrDvo6Mhc+eaUyEf/7jfTfSUk+0Z1m0d\ndqKqnn9oMpnUlpYWzz+XqCxvvtk7hG/TJnNkfdNNJqwnTQrcED6KJhFZr6pJJ+tG4oibImjHjt6n\n8P32t+b89JVXAg88wNljyHoMbgqPbLb3KXwrVpin8I0fz6fwUegwuMlu+afwLVhgnsKXH8J3553m\nVIiFo0GISmFwk522bOl9Ct+GDcCRR/Y+hW/qVKtHgxCVwuAme+zeDTQ1mbB+9lk+hY8ii8FNwXbo\nkHlU6vz55tGpmQxw2mnA/fcDs2cDp5/udwuJhh2Dm4Lp1Vd7bz1/7z2gutqcBrntNuDii3nemiKN\nwU3BsXWrucC4YAHQ1mbOU3/mMyawp00zs8kQEYObfLZ7N7BsGbBwIfDMM+a89YUXAj/8oXnO9ciR\nfreQKHAY3DT8sllgzRpzZP344+a89emnA9/4hjlvzdljiAbF4KbhoQr87nfmyHrJEjPr+YknAnfc\nYcL6wgt53prIIQY3Vdabb5rx1gsXmp+POsrMGjN7tpn1/Mgj/W4hkXUY3OS97dvNVF+LFgHr1pkj\n6cmTga9/HbjhBuD44/1uIZHVGNzkjT17gFTKhPWaNUAuB0yYYJ4TcsstwCmn+N1CotBgcFP5slng\n6adNWD/+OLBvHzB2rJmU4NZbgbPO8ruFRKHE4CZ3urqAF180N8Y89hjwwQfmVvPbbjNhfcklfL41\nUYU5mSz4aADPAziqe/2lqnp/pRtGAaIK/P73JqwXLzYPeEokzHOtb7kFuOoqXmQkGkZOjrgPAJii\nqntEJA7gBRFZpaprK9w2AECqNV3W9Erlvi8MPOv7m2+aoF682DyBLxYzI0EefNCMDDn2WO8b74G+\n/a8eEYcqsDOT5X5QRt/D8j26L9WGxeu2IKeKmAhmXTAGD8w42+9mla1kcKuZ22xP96/x7j/ez3dW\nQKo1jblNbchkcwCAdGcGc5vaAGDQnafc94XBkPueTpsRIYsXA/np5z79aeBHPwJmzjTPug6wgf3f\nsS/b8xr3A3d9D8v36L5UGxau3dzze06153dbw9vRyUgRiYnIKwC2A1ijqusq2yyjsbm9Z6fJy2Rz\naGxur8j7wqCsvr//vnk06mWXAWPGAF/9qjmX3dgIbN4MPP+8mZgg4KENFO5/X9wPnPc9LN+jxeu2\nuFpuA0cXJ1U1B+AcEakGsFxEzlLV1/quIyJzAMwBgLFjx3rSuG2dGVfLh/q+MHDc9507zfC9JUt6\nh+/V1ZnHpd58s/nZQk62MfeD4fuMIMgVmRC92HIbuLr8r6qdAJ4FMLXAa/NUNamqyRqPjsxqqxOu\nlg/1fWEwaN/37DFBff31wMknA5//PLBxoxm+98or5jz2/fdbG9qAs20c+f1gGD8jCGJFHqVQbLkN\nSga3iNR0H2lDRBIArgSwsdINA4CG+jok4rF+yxLxGBrqBw+Wct8XBgP7fnR2P65787d49JmHTFjP\nmmXuZrzrLmDtWuCtt8zFxgkTQvGskELbvq+o7geA+76H5Xs064IxrpbbwMmpklEAfikiMZigf1RV\nV1a2WUb+Aojbq9rlvi8MZpw7GlUH9uO3P1yEi9Y/gyv/tA4jDu43oX3HHeZRqRMnhnas9cBtH9VR\nJV58B8LyPcpfgAzTqBLRCpznSSaT2pIfkUDDY/9+oLkZePRR4IknzHOuR44EbrzRhPWll5rhfEQU\nSCKyXlWTTtblnZM2KxTWJ54IfPazZvLcKVM42zlRCDG4bbNvH7B6tZk1hmFNFEkMbhvs3Qs89RSw\ndCnw5JPm95EjGdZEEcXgDqqdO4GVK4GmJmDVKjO918knm4c5zZxpzlkfwc1HFEX85gdJR4d5PGpT\nE/Df/20emzpqFPDXf23CeuJEXmAkIga377ZsMWG9bJm5tbyrCzj1VOCee8xsMRdcENqhe0RUHga3\nHzZsAJYvN3/ywyY/8Yneqb3OOScUN8MQUWUwuIdDV5cJ6HxYt3c/pOdTnwL+5V/MLehnnOFvG4nI\nGgzuSjl4EHj2WfMgpxUrgG3bzPnpSy8FvvAFMwkB52EkojIwuL20c6cZAZJKmeF7u3cDxxxjJh+Y\nMQOYNs2MuSYiGgIG91C9/bY5on7iCeC554BDh8ywvc99zoT15ZcDRx/tdyuJKEQY3G51dQEvvWTC\nesUK4LXux5KPHw/8/d8D115rRoJw2B4RVQiD24ldu8xkAytXmlMg27ebYJ40CXjoIRPWp5/udyuJ\nKCIY3MX86U8mqFeuNKdAslnghBOAq68256qvvtr8TkQ0zBjceQcOAL/5jTmifvJJ4A9/MMvHjwe+\n/GXgmmuAiy7ibeZE5Ltop1A6bYL6qafMLeZ79gBHHQVMnmyG7E2bBpx2mt+tJCLqp2Rwi8gYAPMB\n/BkABTBPVX9Q6YblpVrT3s3AcfAg8OKLZsje6tVAW5tZPnYsMHu2CerJk80QPgvdl2o7bJaP5MdO\ntHYGk3K3fao1jW898Tp27Mv2LKtOxPHN6Wda03cvOa2jp9+1kApKjUrOgCMiowCMUtWXReQ4AOsB\nzFDVN4q9x6sZcFKtacxtakMmm+tZlojH8OANZzsv1qZNJqRXrwZ+9StzVB2PA5/+NDB1qjlXfeaZ\n1t9ifl+qDQvXbj5seaxKkOvq3cau6+eTcrd9qjWNhqWvIps7fL+OVwkab5oQ+L57yWkdPfmuhVyl\na+RmBpySTy9S1XdV9eXun3cD2ABgWLZkY3N7vyIBQCabQ2Nze/E37dljLijefbeZrfzUU4E77zSz\nmM+ebR7o9MEHJsQbGoCzzrI+tAEzn14hfUMbcFC/gChr23e/r1BoA0C2S63ou5ec1rHcekdJkGrk\n6hy3iIwDcC6AdQVemwNgDgCMHTvWg6YB2zozzpa/8Ya5Aaa5GXjhBTMCJJEALrvMzGZeX29CPAQB\nXUzOxdyhxeoaJI63vcevh43TOpZb7ygJUo0cPy9URI4FsAzAl1R118DXVXWeqiZVNVlTU+NJ42qr\nE86W/+hHwL33miPpL33JjLn+8ENz0fGee8wDnEIc2gAQc9G/YnUNEsfb3uPXw8ZpHcutd5QEqUaO\ngltE4jChvUhVmyrbpF4N9XVIxPvfgZiIx9BQX9d/xa99zTzE6dVXge9+F7jiisjdZj7rgjEFl8eq\n+gd6wfoFkONtX+B98Vjhv8TiVWJF373ktI7l1jtKglQjJ6NKBMDPAGxQ1X+vfJN65U/4l7yKO6Zw\naEXJAzPOBoDQjCpxvO2LvI+jSgyndSy33lESpBo5GVUyEcBvALQB6Ope/HVVfarYe7waVUJEFBVu\nRpWUPOJW1RcAhPsEMRGRRTiZIRGRZRjcRESWYXATEVmGwU1EZBkGNxGRZRjcRESWYXATEVmGwU1E\nZBkGNxGRZRjcRESWYXATEVmGwU1EZBkGNxGRZRjcRESWYXATEVmmZHCLyM9FZLuIvDYcDSIiosE5\nmeX9FwD+P4D5lW1Kf6nWdM8UQccn4sjmurD3YA6As2mo+r4/CtMwue1vmOoTpr54hTUpz8C6TT6j\nBs9u7AhcHUtOXQYAIjIOwEpVPcvJhw516rJUaxpzm9qQyeaKrhOvEjTeNKFgEQu9PxGP4cEbzg5E\n0b3mtr9hqk+Y+uIV1qQ8TnKnknV0M3VZIM9xNza3D1o8AMh2KRqb2x2/P5PNFV3fdm77G6b6hKkv\nXmFNyuMkd4JSR8+CW0TmiEiLiLR0dHQM6bO2dWaGtJ7b5bbzqg421idMffEKa1KeoebOcPIsuFV1\nnqomVTVZU1MzpM+qrU4MaT23y23nVR1srE+Y+uIV1qQ8Q82d4RTIUyUN9XVIxGODrhOvEjTU1zl+\nfyIeK7q+7dz2N0z1CVNfvMKalMdJ7gSljiVHlYjIYgCXAThJRLYCuF9Vf1bJRuVP/Jc7qmTg+4N0\nNbgS3PY3TPUJU1+8wpqUp1DdrB5V4tZQR5UQEUWN9aNKiIioOAY3EZFlGNxERJZhcBMRWYbBTURk\nGQY3EZFlGNxERJZhcBMRWYbBTURkGQY3EZFlGNxERJZhcBMRWYbBTURkGQY3EZFlGNxERJZhcBMR\nWcZRcIvIVBFpF5E3ReTeSjeKiIiKczJ1WQzADwFcCWArgJdEZIWqvlHpxhWTak2jsbkd6SKzLZ8w\nIo77r+0/tVn+PUGbgqhcfftTPSKOA9kc9mW7Dluv0DRv96XasHjdFuRUERPBrAvG4IEZZw9n8ysu\n1ZrGN1e8js5MtmdZof0iqsL2fXCrVP+DXp+SwQ3gfABvqupbACAiSwBcB8CX4E61pjG3qQ2ZbK7o\nOjv2ZdGw9FUAZh65ge9Jd2Ywt6mt53XbDOzPjn3Zout2ZrJoeKy3Fvel2rBw7eae13OqPb+HJbxT\nrWk0PPYqsl39p+UbuF9EVdi+D26V6r8N9XFyqmQ0gC19ft/avcwXjc3tg4Z2XjanaGxuL/qeTDbX\n87ptnNYgL9vVW4vF67YUXKfYchs1NrcfFtp5ffeLqArb98GtUv23oT6eXZwUkTki0iIiLR0dHV59\n7GG2FTk9Mti6xd7j5rOCpJx259+TKzI5dLHlNipVH1u3u1fC9n1wq1T/baiPk+BOAxjT5/dTupf1\no6rzVDWpqsmamhqv2neY2uqE63WLvcfNZwVJOe3OvycmUvD1YsttVKo+tm53r4Tt++BWqf7bUB8n\nwf0SgP8nIqeKyJEAbgaworLNKq6hvg6JeKzkevGYoKG+ruh7EvFYz+u2cVqDvHhVby1mXTCm4DrF\nltuoob4O8arCfxH13S+iKmzfB7dK9d+G+pS8OKmqh0TkCwCaAcQA/FxVX694y4rIXxxwM6qk73uC\nepXYjYH9cTOqJH8BMsyjSvJ95aiSwsL2fXCrVP9tqI9oBc5tJpNJbWlp8fxziYjCSkTWq2rSybq8\nc5KIyDIMbiIiyzC4iYgsw+AmIrIMg5uIyDIVGVUiIh0A3vH8gwd3EoD3h/n/GWSsR3+sR3+sR39B\nqMfHVNXR3YsVCW4/iEiL06E0UcB69Md69Md69GdbPXiqhIjIMgxuIiLLhCm45/ndgIBhPfpjPfpj\nPfqzqh6hOcdNRBQVYTriJiKKBOuDW0TGiMizIvKGiLwuIvf43SY/icjRIvI7EXm1ux7f8rtNQSAi\nMRFpFZGVfrfFbyKySUTaROQVEYn80+BEpFpElorIRhHZICIX+d2mUpzMORl0hwB8VVVfFpHjAKwX\nkTV+TmbsswMApqjqHhGJA3hBRFap6lq/G+azewBsAPARvxsSEJNV1e9xy0HxAwCrVXVm95wDI/xu\nUCnWH3Gr6ruq+nL3z7thvpzBeXDuMFNjT/ev8e4/kb6QISKnAJgG4Kd+t4WCRUSOBzAJwM8AQFUP\nqmqnv60qzfrg7ktExgE4F8A6f1vir+7TAq8A2A5gjapGuh4Avg/gHwAcPtNENCmAp0VkvYjM8bsx\nPjsVQAeA/+o+lfZTETnG70aVEprgFpFjASwD8CVV3eV3e/ykqjlVPQdmftDzReQsv9vkFxG5BsB2\nVV3vd1sCZKKqfhLA1QD+TkQm+d0gHx0B4JMAfqyq5wLYC+Bef5tUWiiCu/tc7jIAi1S1ye/2BEX3\nP/meBTAq6HHCAAABA0lEQVTV77b46BIA00VkE4AlAKaIyEJ/m+QvVU13/3c7gOUAzve3Rb7aCmBr\nn3+VLoUJ8kCzPrhFRGDOT21Q1X/3uz1+E5EaEanu/jkB4EoAG/1tlX9Uda6qnqKq42Amun5GVWf7\n3CzfiMgx3Rfx0X1K4CoAr/nbKv+o6nsAtohIfibgywEEfmBDGEaVXALgNgBt3ed1AeDrqvqUj23y\n0ygAvxSRGMxfzI+qauSHwFGPPwOw3Bzv4AgAD6vqan+b5LsvAljUPaLkLQB3+NyeknjnJBGRZaw/\nVUJEFDUMbiIiyzC4iYgsw+AmIrIMg5uIyDIMbiIiyzC4iYgsw+AmIrLM/wEB5GmljH4wSwAAAABJ\nRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10a39f668>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "df_test = pd.DataFrame()\n",
    "df_test[\"x\"]=np.arange(2,6,0.1)\n",
    "plt.scatter(d.x+d.x_bure,d.y)\n",
    "plt.plot(df_test.x,result.predict(df_test)*8,c=\"red\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* 真のモデルと異なる結果となった。\n",
    "* 改めて、データを確認してみよう。$x_i=4$である20個体の生存種子数は下記。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "y\n",
       "0    3\n",
       "1    1\n",
       "2    4\n",
       "3    2\n",
       "4    1\n",
       "5    1\n",
       "6    2\n",
       "7    3\n",
       "8    3\n",
       "Name: N, dtype: int64"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "d[d.x == 4].groupby([d.y]).count()[\"N\"]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 図7.3.(B)を作成する（できてない）"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from collections import Counter\n",
    "df_test_2 = pd.DataFrame()\n",
    "df_test_2[\"x\"]= np.array([4])\n",
    "p = result.predict(df_test_2)\n",
    "# binomi = np.random.binomial(n=8,p=p,size=10000)\n",
    "# ys,counts = zip(*Counter(binomi).items())\n",
    "# indexes=np.arange(len(ys))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "ename": "TypeError",
     "evalue": "'module' object is not callable",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------\u001b[0m",
      "\u001b[0;31mTypeError\u001b[0m               Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-11-a00d92c339e2>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[1;32m      2\u001b[0m \u001b[0mp_list\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mlist\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      3\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0mi\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mrange\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m9\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 4\u001b[0;31m     \u001b[0mp_list\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mappend\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mscm\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m8\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mi\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;31mTypeError\u001b[0m: 'module' object is not callable"
     ]
    }
   ],
   "source": [
    "import scipy.misc as scm\n",
    "p_list = list()\n",
    "for i in range(9):\n",
    "    p_list.append(scm(8,i,1)*)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'indexes' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------\u001b[0m",
      "\u001b[0;31mNameError\u001b[0m               Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-12-2e00ebb4a14e>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[1;32m      1\u001b[0m \u001b[0mplt\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mplot\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0md\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0md\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mx\u001b[0m \u001b[0;34m==\u001b[0m \u001b[0;36m4\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mgroupby\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0md\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0my\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcount\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m\"N\"\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mindex\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0md\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0md\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mx\u001b[0m \u001b[0;34m==\u001b[0m \u001b[0;36m4\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mgroupby\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0md\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0my\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcount\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\"o\"\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 2\u001b[0;31m \u001b[0mplt\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mplot\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mindexes\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0marray\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mcounts\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m/\u001b[0m\u001b[0;36m500\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      3\u001b[0m \u001b[0mplt\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mylim\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m6\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mNameError\u001b[0m: name 'indexes' is not defined"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAEZxJREFUeJzt3X+sX3d93/Hny/b1NRjmLDdXa5YfGKkR0oomTK8CjKVC\nyZiARmHSqBqk4nZi8mjXLay1qtZ/UMEfkSpZrKoiUaUJI14plAVaZSh0jUSqGmkErn8ESMymjNEm\nWbbcXhcHW/j613t/3JPKvti533vv9/qc+8nzIX2V8z3fj7+fl3yc1z3+fM/XJ1WFJKktm/oOIEka\nP8tdkhpkuUtSgyx3SWqQ5S5JDbLcJalBlrskNchyl6QGWe6S1KAtfU183XXX1c6dO/uaXpI2pEOH\nDv1NVU0vN663ct+5cyezs7N9TS9JG1KSvxplnMsyktQgy12SGmS5S1KDLHdJapDlLkkNGvlqmSSb\ngVng+aq6c8lrk8AB4KeBeeDnq+r7Y8ypVTqwdx+n53dxduu1TJw5zrapI+zef2/fsSSts5Wcud8D\nHLvCax8G/raqfhL4j8DvrDWY1u7A3n2cOnEbZyenIOHs5BSnTtzGgb37+o4maZ2NVO5JbgR+Fnjg\nCkPeDzzUbT8M3JEka4+ntTg9v4sLmycv2Xdh8ySn53f1lEjS1TLqmfvvAr8BXLjC6zcAzwJU1Tng\nBDC1dFCSPUlmk8zOzc2tIq5W4uzWa1e0X1I7li33JHcCL1bVobVOVlX3V9VMVc1MTy/77Vmt0cSZ\n4yvaL6kdo5y5vxO4K8n3gc8Dtyf5wyVjngduAkiyBdjB4ger6tG2qSNsOr9wyb5N5xfYNnWkp0SS\nrpZly72qfquqbqyqncDdwFer6heWDHsE+MVu+wPdmBprUq3Y7v33sn3HQSYW5qGKiYV5tu846NUy\n0qvAqv/hsCSfAGar6hHgQeA/J3kGOM7iDwENwI8X+c/1kkPS1bWicq+qvwD+otv+2EX7T2NrSNJg\n+A1VSWqQ5S5JDbLcJalBlrskNchyl6QGWe6S1CDLXZIaZLlLUoMsd0lqkOUuSQ2y3CWpQZa7JDXI\ncpekBlnuktQgy12SGmS5S1KDLHdJapDlLkkNstwlqUGWuyQ1yHKXpAZZ7pLUIMtdkhpkuUtSg5Yt\n9yTbknwjyZNJnkry8cuM+aUkc0mOdo9/vT5xJUmj2DLCmAXg9qo6mWQC+FqSr1TV15eM++Oq+tXx\nR5QkrdSy5V5VBZzsnk50j1rPUJKktRlpzT3J5iRHgReBx6rqicsM+5dJvpXk4SQ3XeF99iSZTTI7\nNze3htiSpFcyUrlX1fmqegtwI3BrkjcvGfJfgZ1V9Y+Bx4CHrvA+91fVTFXNTE9PryW3JOkVrOhq\nmar6AfA48J4l++eraqF7+gDw0+OJJ0lajVGulplOck23/Rrg3cB3l4y5/qKndwHHxhlSkrQyo1wt\ncz3wUJLNLP4w+EJVfTnJJ4DZqnoE+PdJ7gLOAceBX1qvwJKk5WXxYpirb2ZmpmZnZ3uZW5I2qiSH\nqmpmuXF+Q1WSGmS5S1KDLHdJapDlLkkNstwlqUGWuyQ1yHKXpAZZ7pLUIMtdkhpkuUtSgyx3SWqQ\n5S5JDbLcJalBlrskNchyl6QGWe6S1CDLXZIaZLlLUoMsd0lqkOUuSQ2y3CWpQZa7JDXIcpekBi1b\n7km2JflGkieTPJXk45cZM5nkj5M8k+SJJDvXI6wkaTRbRhizANxeVSeTTABfS/KVqvr6RWM+DPxt\nVf1kkruB3wF+ftxhD+zdx+n5XZzdei0TZ46zbeoIu/ffO+5pJG1wQ+2Kq5lr2TP3WnSyezrRPWrJ\nsPcDD3XbDwN3JMnYUrL4m3LqxG2cnZyChLOTU5w6cRsH9u4b5zSSNrihdsXVzjXSmnuSzUmOAi8C\nj1XVE0uG3AA8C1BV54ATwNQ4g56e38WFzZOX7LuweZLT87vGOY2kDW6oXXG1c41U7lV1vqreAtwI\n3JrkzauZLMmeJLNJZufm5lb0a89uvXZF+yW9Og21K652rhVdLVNVPwAeB96z5KXngZsAkmwBdgDz\nl/n191fVTFXNTE9PryjoxJnjK9ov6dVpqF1xtXONcrXMdJJruu3XAO8Gvrtk2CPAL3bbHwC+WlVL\n1+XXZNvUETadX7hk36bzC2ybOjLOaSRtcEPtiquda5Qz9+uBx5N8C/gmi2vuX07yiSR3dWMeBKaS\nPAP8GvCb4w66e/+9bN9xkImFeahiYmGe7TsODuITcEnDMdSuuNq5MuYT7JHNzMzU7OxsL3NL0kaV\n5FBVzSw3zm+oSlKDLHdJapDlLkkNstwlqUGWuyQ1yHKXpAZZ7pLUIMtdkhpkuUtSgyx3SWqQ5S5J\nDbLcJalBlrskNchyl6QGWe6S1CDLXZIaZLlLUoMsd0lqkOUuSQ2y3CWpQZa7JDXIcpekBlnuktQg\ny12SGrRsuSe5KcnjSZ5O8lSSey4z5l1JTiQ52j0+tj5xJUmj2DLCmHPAr1fV4SSvBw4leayqnl4y\n7mBV3Tn+iJKklVr2zL2qXqiqw932D4FjwA3rHUyStHorWnNPshPYBTxxmZffkeTJJF9J8lNjyCZJ\nWqVRlmUASPI64IvAR6vqpSUvHwbeUFUnk7wP+FPglsu8xx5gD8DNN9+86tCSpFc20pl7kgkWi/2z\nVfWlpa9X1UtVdbLbfhSYSHLdZcbdX1UzVTUzPT29xuiSpCsZ5WqZAA8Cx6rqk1cY8xPdOJLc2r3v\n/DiDSpJGN8qyzDuBDwHfTnK027cPuBmgqn4f+ADwy0nOAT8C7q6qWoe8kqQRLFvuVfU1IMuMuQ+4\nb1yhJElr4zdUJalBlrskNchyl6QGWe6S1CDLXZIaZLlLUoMsd0lqkOUuSQ2y3CWpQZa7JDXIcpek\nBlnuktQgy12SGmS5S1KDLHdJapDlLkkNstwlqUGWuyQ1yHKXpAZZ7pLUIMtdkhpkuUtSgyx3SWqQ\n5S5JDVq23JPclOTxJE8neSrJPZcZkyS/l+SZJN9K8tb1iStJGsWWEcacA369qg4neT1wKMljVfX0\nRWPeC9zSPd4GfKr7r3RZB/bu4/T8Ls5uvZaJM8fZNnWE3fvv7TuWVsjjOFzLnrlX1QtVdbjb/iFw\nDLhhybD3Awdq0deBa5JcP/a0asKBvfs4deI2zk5OQcLZySlOnbiNA3v39R1NK+BxHLYVrbkn2Qns\nAp5Y8tINwLMXPX+OH/8BIAFwen4XFzZPXrLvwuZJTs/v6imRVsPjOGwjl3uS1wFfBD5aVS+tZrIk\ne5LMJpmdm5tbzVuoAWe3Xrui/Romj+OwjVTuSSZYLPbPVtWXLjPkeeCmi57f2O27RFXdX1UzVTUz\nPT29mrxqwMSZ4yvar2HyOA7bKFfLBHgQOFZVn7zCsEeA3d1VM28HTlTVC2PMqYZsmzrCpvMLl+zb\ndH6BbVNHekqk1fA4DtsoZ+7vBD4E3J7kaPd4X5KPJPlIN+ZR4HvAM8AfAL+yPnHVgt3772X7joNM\nLMxDFRML82zfcdCrLDYYj+Owpap6mXhmZqZmZ2d7mVuSNqokh6pqZrlxfkNVkhpkuUtSgyx3SWqQ\n5S5JDbLcJalBlrskNchyl6QGWe6S1CDLXZIaZLlLUoMsd0lqkOUuSQ2y3CWpQZa7JDXIcpekBlnu\nktQgy12SGmS5S1KDLHdJapDlLkkNstwlqUGWuyQ1yHKXpAZZ7pLUoGXLPcmnk7yY5DtXeP1dSU4k\nOdo9Pjb+mJKkldgywpjPAPcBB15hzMGqunMsiSRJa7bsmXtV/SVw/CpkkSSNybjW3N+R5MkkX0ny\nU2N6T0nSKo2yLLOcw8AbqupkkvcBfwrccrmBSfYAewBuvvnmMUwtSbqcNZ+5V9VLVXWy234UmEhy\n3RXG3l9VM1U1Mz09vdapJUlXsOZyT/ITSdJt39q95/xa31eStHrLLssk+RzwLuC6JM8Bvw1MAFTV\n7wMfAH45yTngR8DdVVXrlliStKxly72qPrjM6/exeKmkJGkg/IaqJDXIcpekBlnuktQgy12SGmS5\nS1KDLHdJapDlLkkNstwlqUGWuyQ1yHKXpAZZ7pLUIMtdkhpkuUtSgyx3SWqQ5S5JDbLcJalBlrsk\nNchyl6QGWe6S1CDLXZIaZLlLUoMsd0lqkOUuSQ2y3CWpQcuWe5JPJ3kxyXeu8HqS/F6SZ5J8K8lb\nxx9TkrQSW0YY8xngPuDAFV5/L3BL93gb8Knuv68aB/bu4/T8Ls5uvZaJM8fZNnWE3fvv7TuWVmGo\nx3KouTRcy565V9VfAsdfYcj7gQO16OvANUmuH1fAoTuwdx+nTtzG2ckpSDg7OcWpE7dxYO++vqNp\nhYZ6LIeaS8M2jjX3G4BnL3r+XLfvVeH0/C4ubJ68ZN+FzZOcnt/VUyKt1lCP5VBzadiu6geqSfYk\nmU0yOzc3dzWnXjdnt167ov0arqEey6Hm0rCNo9yfB2666PmN3b4fU1X3V9VMVc1MT0+PYer+TZy5\n/IrVlfZruIZ6LIeaS8M2jnJ/BNjdXTXzduBEVb0whvfdELZNHWHT+YVL9m06v8C2qSM9JdJqDfVY\nDjWXhm2USyE/B/x34E1Jnkvy4SQfSfKRbsijwPeAZ4A/AH5l3dIO0O7997J9x0EmFuahiomFebbv\nOOiVDBvQUI/lUHNp2FJVvUw8MzNTs7OzvcwtSRtVkkNVNbPcOL+hKkkNstwlqUGWuyQ1yHKXpAZZ\n7pLUoN6ulkkyB/zVKn/5dcDfjDHOuAw1Fww3m7lWxlwr02KuN1TVst8C7a3c1yLJ7CiXAl1tQ80F\nw81mrpUx18q8mnO5LCNJDbLcJalBG7Xc7+87wBUMNRcMN5u5VsZcK/OqzbUh19wlSa9so565S5Je\nwYYr9yTvSfI/uhty/2bfeWD5m4j3JclNSR5P8nSSp5Lc03cmgCTbknwjyZNdro/3neliSTYnOZLk\ny31neVmS7yf5dpKjSQbzL+4luSbJw0m+m+RYkncMINObut+nlx8vJflo37kAkvyH7s/8d5J8Lsm2\ndZtrIy3LJNkM/E/g3Szezu+bwAer6umec/0McJLFe8m+uc8sF+vuZXt9VR1O8nrgEPAvBvD7FWB7\nVZ1MMgF8Dbinuwdv75L8GjAD/L2qurPvPLBY7sBMVQ3qmu0kDwEHq+qBJFuB11bVD/rO9bKuM54H\n3lZVq/1ezbiy3MDin/V/VFU/SvIF4NGq+sx6zLfRztxvBZ6pqu9V1Rng8yzeoLtXI9xEvBdV9UJV\nHe62fwgcYwD3t+1upn6yezrRPQZxlpHkRuBngQf6zjJ0SXYAPwM8CFBVZ4ZU7J07gP/Vd7FfZAvw\nmiRbgNcC/2e9Jtpo5f6qvhn3WiTZCewCnug3yaJu6eMo8CLwWFUNIhfwu8BvABf6DrJEAX+e5FCS\nPX2H6bwRmAP+U7eM9UCS7X2HWuJu4HN9hwCoqueB/cBfAy+weNe6P1+v+TZauWsVkrwO+CLw0ap6\nqe88AFV1vqrewuI9d29N0vtyVpI7gRer6lDfWS7jn1bVW4H3Av+2Wwrs2xbgrcCnqmoXcAoYxOdg\nAN0y0V3Af+k7C0CSv8/iSsMbgX8IbE/yC+s130Yr95Fvxq1F3Zr2F4HPVtWX+s6zVPfX+MeB9/Sd\nBXgncFe3vv154PYkf9hvpEXdWR9V9SLwJywuUfbtOeC5i/7W9TCLZT8U7wUOV9X/6ztI558B/7uq\n5qrqLPAl4J+s12Qbrdy/CdyS5I3dT+W7WbxBty6j++DyQeBYVX2y7zwvSzKd5Jpu+zUsfkD+3X5T\nQVX9VlXdWFU7Wfyz9dWqWrczq1El2d59IE637PHPgd6vzKqq/ws8m+RN3a47gF4/rF/igwxkSabz\n18Dbk7y2+3/zDhY/B1sXW9brjddDVZ1L8qvAfwM2A5+uqqd6jvXyTcTfBVyX5Dngt6vqwX5TAYtn\noh8Cvt2tbwPsq6pHe8wEcD3wUHclwybgC1U1mMsOB+gfAH+y2AdsAf6oqv6s30h/598Bn+1Otr4H\n/Kue8wB/90Pw3cC/6TvLy6rqiSQPA4eBc8AR1vGbqhvqUkhJ0mg22rKMJGkElrskNchyl6QGWe6S\n1CDLXZIaZLlLUoMsd0lqkOUuSQ36/z89vlUToUvCAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10a413e10>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(d[d.x == 4].groupby([d.y]).count()[\"N\"].index,d[d.x == 4].groupby([d.y]).count(),\"o\")\n",
    "plt.plot(indexes,np.array(counts)/500,)\n",
    "plt.ylim(0,6)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* 見て見ると二項分布に従ってはいないようだ。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 7.2 過分散と個体差"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* 「N個のうちy個が生存した」というカウントデータを扱っているのに、そのばらつきが二項分布で説明できないように見える。\n",
    "* 二項分布で期待されるよりも大きなばらつきを過分散（または過大分散,overdispresion)と言う。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 7.2.1. 過分散：ばらつきが大きすぎる"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* 葉数$x_i=4$の個体の生存種子数のデータが、どのように「二項分布に合っていない」のかを調べる。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "d_4 = pd.DataFrame()\n",
    "d_4 = d[d.x==4]\n",
    "# d_4[\"seizon\"] = d[d.x == 4].groupby([d.y]).count()[\"N\"].index\n",
    "# d_4[\"samplesuu\"] = d[d.x == 4].groupby([d.y]).count()[\"N\"].values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x10a6cb240>]"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl0VPd99/H3V/suIWkGhFgESBoh8E4INjbWeEkwBqdL\n2ibPk6RJ08d167Rxk57UiRNck7pZ+mQ5WdrUjdMkT90krZO0AWPHdgzGu40x2BpJCLGjQUhCMKN9\nGf2ePzTCCpbQSJqZO/fO93XOnDPSXM39HCG+uvre3yLGGJRSSjlLitUBlFJKRZ8Wd6WUciAt7kop\n5UBa3JVSyoG0uCullANpcVdKKQfS4q6UUg6kxV0ppRxIi7tSSjlQmlUnLi0tNRUVFVadXimlbOn1\n11/vNMa4pjvOsuJeUVHB3r17rTq9UkrZkogcj+Q4bcsopZQDaXFXSikH0uKulFIOpMVdKaUcSIu7\nUko5UMTFXURSReQNEdkxyWuZIvIzEWkRkVdEpCKaIZVSSs3MTK7cPwk0TvHax4FzxphK4BvAV+Ya\nTCml1OxFVNxFZBFwO/D9KQ55H/Cj8PNHgZtFROYeT83VcGiU/3jlBEMjo1ZHUUrFUaRX7t8EPgNM\nVSHKgZMAxpgRIACUXHyQiNwpIntFZG9HR8cs4qqZ2n2wg8/98i22H/BbHUUpFUfTFncR2Qy0G2Ne\nn+vJjDEPGWPWGGPWuFzTzp5VUVDfGgBgd7P+MlUqmURy5b4euENEjgE/BW4SkX+/6JhWYDGAiKQB\nhcDZKOZUs+TzBwHY09zBSEhbM0oli2mLuzHms8aYRcaYCuADwDPGmA9ddNivgD8OP39/+BgT1aRq\nVhr8AQqz0wn0D7P/5Hmr4yil4mTW49xFZJuI3BH+8GGgRERagE8B90YjnJqbrt4h/IEBPrxuKakp\nwq6D7VZHUkrFyYyKuzFmtzFmc/j5VmPMr8LPB4wxf2CMqTTGrDXGHIlFWDUzPv9Yv/3aFSVcs3Qe\nuw9q312pZKEzVB1svN++amEBXo8bnz/ImeCAxamUUvGgxd3BfP4g5UXZFOVk4K0ZG530rF69K5UU\ntLg7mK81wKqFBQB45uezoCBL++5KJQkt7g7VOzjC0bO9rFpYCICI4K1x8dyhToZ1SKRSjqfF3aEa\nTwcxhgtX7gB1Hjc9gyPsPXbOwmRKqXjQ4u5Q4zNTV5cXXvjc+spS0lOF3dqaUcrxtLg7lM8fpCQ3\ng/kFmRc+l5eZxtplxdp3VyoJaHF3KJ8/SO3CAi5enNPrcdN8podT5/osSqaUigct7g40OBLiUHv3\nhZupE9V53AA6oUkph9Pi7kCHzvQwHDKsLi94x2srXLksLs7W4q6Uw2lxd6DxZQcmu3IXEbweNy+0\ndDI4Eop3NKVUnGhxdyCfP0heZhpLi3Mmfd3rcdM/HOLVo11xTqaUihct7g7k8wepLSsgJWXynQ7X\nLS8hIy2FXU3amlHKqbS4O0xo1NAQHikzleyMVK5dXqLj3ZVyMC3uDnO0s5f+4dBvzUydjNfj4khn\nL8c6e+OUTCkVT1rcHeZSN1MnentIpF69K+VEWtwdxucPkpGaQtX8vEseV1Gay/LSXHbpkEilHEmL\nu8P4/AE8C/JJT53+n7bO4+alI2fpH9IhkUo5jRZ3BzHG4PMHp+23j/PWuBgaGeWlI50xTqaUijct\n7g7iDwxwvm844uK+dlkx2empOiRSKQfS4u4g48v8riq/9M3UcZlpqayvLGXXwXaMMbGMppSKMy3u\nDuLzB0kRWLkgsit3GGvNnDrXz+EOHRKplJNocXeQBn+A5a48sjNSI/4aHRKplDNpcXcQnz/I6gj7\n7ePKi7Kpnp+nG3go5TDTFncRyRKRV0XkgIj4ROSBSY75qIh0iMj+8ONPYxNXTeVszyCnAwPTTl6a\njNfj5tWjXfQMjsQgmVLKCpFcuQ8CNxljrgCuBDaKyLpJjvuZMebK8OP7UU2ppuXzBwEiHikzUZ3H\nzXDI8EKLDolUyimmLe5mTE/4w/TwQ4dWJJjx4n6pBcOmsqZiHnmZadp3V8pBIuq5i0iqiOwH2oGn\njDGvTHLY74vImyLyqIgsnuJ97hSRvSKyt6NDx1ZHU70/wKJ52RTlZMz4a9NTU7ihqpRdTR06JFIp\nh4iouBtjQsaYK4FFwFoRWX3RIduBCmPM5cBTwI+meJ+HjDFrjDFrXC7XXHKrizTMYGbqZLweN23B\nAZrauqOYSilllRmNljHGnAd2ARsv+vxZY8xg+MPvA9dEJ56KRM/gCEc7e2d1M3XcjZ6xX7Y6akYp\nZ4hktIxLRIrCz7OBW4Gmi44pm/DhHUBjNEOqS2s8PfubqePmF2SxamGBbpytlENEcuVeBuwSkTeB\n1xjrue8QkW0ickf4mL8KD5M8APwV8NHYxFWTGV92YHWEyw5Mxetx8/rxcwT6h6MRSyllobTpDjDG\nvAlcNcnnt054/lngs9GNpiLl8wcpzcvAnZ85p/ep87j4zq4Wnj/Uye2Xl03/BUqphKUzVB3A5w9S\nu7AQkck3xI7UlYuLKMxO1767Ug6gxd3mBkdCHDrTPeNlByaTlprChmoXuw92MDqqQyKVsjMt7jbX\n3NbDyKiZ00iZibweF509gxcmRSml7EmLu829vSH23K/cATZUuxDRIZFK2Z0Wd5vz+YPkZaaxpDgn\nKu9XmpfJ5YuKtLgrZXNa3G2u3h+gdmEBKSlzu5k6kdfjYv/J83T1DkXtPZVS8aXF3cZCo4am091R\na8mM83rcGAN7mnVCk1J2pcXdxo529tA/HIrazdRxl5UXUpKboa0ZpWxMi7uNjY9oWV0e3Sv3lBTh\nRo+LPc0dhHRIpFK2pMXdxupbA2SkpbDClRf1967zuDnXN8yBU+ej/t5KqdjT4m5jPn+QmgX5pKdG\n/59xQ1UpKQK7m7Q1o5QdaXG3KWMMvjmu4X4pRTkZXL1kHrt0lUilbEmLu021nu8n0D8c9ZupE3lr\n3LzVGqC9eyBm51BKxYYWd5uqb537Gu7TqQtv4PGsXr0rZTta3G2qwR8gRaBmQeyKe21ZAe78TN3A\nQykb0uJuUz5/kBWuPLIzUmN2DhHB63Gz51AHw6HRmJ1HKRV9Wtxtqt4fmPPOS5Hw1rjoHhhh3/Fz\nMT+XUip6tLjbUGfPIGeCgzHtt49bX1lKWoroqBmlbEaLuw2Nz0ytjUNxz89K510VxezWpQiUshUt\n7jb09hrusW/LwNiomaa2bk4H+uNyPqXU3GlxtyFfa5DFxdkUZqfH5XzeGjeAjppRyka0uNuQzx9g\nVVl8rtoBqtx5lBdls0uXIlDKNrS420z3wDDHzvbF5WbqOBGhzuPihZZOBkdCcTuvUmr2pi3uIpIl\nIq+KyAER8YnIA5MckykiPxORFhF5RUQqYhFWQePpboC4DIOcyOtx0zsUYu8xHRKplB1EcuU+CNxk\njLkCuBLYKCLrLjrm48A5Y0wl8A3gK9GNqcbVt0Z3Q+xIXVdZQkZqirZmlLKJaYu7GdMT/jA9/Lh4\nB4f3AT8KP38UuFlEorepp7rA5w9SmpeJuyArrufNyUjj3cuLdXcmpWwiLZKDRCQVeB2oBL5rjHnl\nokPKgZMAxpgREQkAJUBnFLMC8MaJc3z/uaP84x9cTk5GRPEdxecPxP2qfZzX42bbjgZOnO1jSUmO\nJRmUioQxhof2HKGlvWf6gy3grXGz6bKymJ4joupojAkBV4pIEfBLEVltjKmf6clE5E7gToAlS5bM\n9MsB6BsK8dhbp/ndq8q5pXb+rN7DrgaGQxxq7+HmlW5Lzu+tGSvuu5vb+ci1FZZkUCoS2988zZce\nb8Kdn0laSuI1ESpKc2N+jhld+hpjzovILmAjMLG4twKLgVMikgYUAmcn+fqHgIcA1qxZM6vNOddU\nzCM3I5VdB9uTrrg3n+kmNGriNnnpYstKc6koyWFXkxZ3lbj6h0J8aWcjq8sL+J+7ryc1AYt7PEQy\nWsYVvmJHRLKBW4Gmiw77FfDH4efvB54xxsRkZ+XMtFTWV5ay+2AHMTpFwhpfdsCqtgyM7a360pGz\nDAzrkEiVmL737GFOBwbYunlV0hZ2iGy0TBmwS0TeBF4DnjLG7BCRbSJyR/iYh4ESEWkBPgXcG5u4\nY7w1blrP93MoQftpseLzB8jPSmNJsXX97jqPi4HhUV4+8o4/zJSyXOv5fr737GE2X17G2mXFVsex\n1LRtGWPMm8BVk3x+64TnA8AfRDfa1MZ3CNrV1E71/Px4ndZy9a1BassKsHIg0rrlJWSlp7D7YAd1\nHmt6/0pN5cuPjzUVPrtppcVJrGfLGaplhdnULMhPqmF5oVFDU1vQsn77uKz0VK5bUcozTe1J1xZT\nie3Vo11sP+DnrhtXUF6UbXUcy9myuMNYa2bvsXN0DwxbHSUujnT0MDA8amm/fZzX4+JEVx9HO3ut\njqIUMHbx88B2H2WFWdx14wqr4yQE+xZ3j5uRUcMLLVEfSp+Q6sPL/MZ72YHJjLdjdAMPlSgeff0k\nPn+Qe2+rienWk3Zi2+J+9ZIi8rPS2NWUHAXG1xokMy2FFa7Yj4+dzuLiHCrdebqBh0oI3QPD/OOv\nD7Jm6TzuuGKh1XEShm2Le1pqChuqXew6mBy9X58/SM2CfNJSE+OfzOtx8cqRLnoHR6yOopLcd55p\n4WzvEPdvWWXpYINEkxiVYpbqql20dw/ScDpodZSYMsbg8weotfhm6kRej5uh0CgvHtYhkco6Rzt7\n+cELR/mDaxZx2aLE+f+RCGxd3G8MD4l0+g5Bp871ExwYYXW59TdTx62pKL4wU1gpqzz4WAOZaan8\nzXs9VkdJOLYu7u78LC4rL3T8MrTx3jM1EhlpKWMzhXVIpLLIs80dPN3YziduqsSdH99VUu3A1sUd\nxnq/+06c43zfkNVRYsbnD5KaItQsSKwJW94aN/7AQNLNFFbWGw6N8sUdDSwtyeFj6yusjpOQbF/c\n62rcjBrYc8i5QyJ9/iCVrjyy0hNriNfEmcJKxdMjLx+npb2Hz99eS2ZaYv2/SBS2L+5XLCpiXk46\nux1cYOpbrVvD/VKScaawsl5X7xBff6qZG6pKucWi5a/twPbFPTVFuLHaxe7mDkZHndf77egepL17\nkNoELO7w9kzhYJLMFFbW+8ZTzfQOhfjC5lod+ngJti/uMFZgunqHeDO8v6iTJOLN1IkuzBR2cFtM\nJY6mtiCPvHKcD717SVItGjgbjijuG6pciODIGZPja7gn6pX7hZnCDvzeq8RijGHb9gYKstP561ur\nrY6T8BxR3OflZnDV4iJHrnXi8wdYUpxDYXa61VEm9fZM4eTbPEXF15MNZ3jx8Fk+dWs1RTkZVsdJ\neI4o7jDWHnjz1Hk6ewatjhJVPn8wIW+mTuT1uOnoHrzwV4ZS0TYwHOLBxxqpnp/H/1o7u/2Xk41j\ninudx40xsKfZOVfvwYFhjp/tS/jifmP1+Exhbc2o2PjBC0c50dXH/VtWJcz6SonOMd+lVQsLKM3L\ndFRrpmF8z9QEWOb3Ulz5mWMzhR30vVeJ40xwgO8808KttfNZX1lqdRzbcExxT0kR6jwu9jR3MBIa\ntTpOVCTChtiR8npcvOHwmcLKGl994iAjIcN9unXejDimuMNY7zfQP8z+k+etjhIVPn8AV36mLdbN\nSIaZwir+Dpw8z8/3neJPrl9GRan1exnYiaOK+/VVpaSmiGOG5TX4g6y2wVU7JMdMYRVfxhj+brsP\nV34mn7ip0uo4tuOo4l6Ync41S+c5YnemgeEQh9p7Enby0sWcPlNYxd//7PfzxonzfOa9HvIy06yO\nYzuOKu4w1pppOB2kLTBgdZQ5OdjWTWjU2KLfPs7JM4VVfPUOjvClxxu5fFEhv3/1Iqvj2JLzinvN\n2LC8Z5vt3R54+2aqPa7c4e2ZwrpKpJqr7z17mDPBQe7fUktKiq4fMxvTFncRWSwiu0SkQUR8IvLJ\nSY6pE5GAiOwPP7bGJu70PPPzKSvMsv3uTPX+APlZaSwuzrY6SsTGZwrreHc1Fye7+nhozxHed+VC\nrllabHUc24rkyn0E+LQxphZYB9wtIrWTHPecMebK8GNbVFPOgIhQ53Hz3KFOhm08JHJ8ZqrdVr3z\netwcOBWgo9tZM4VV/Hz58SZSRLj3thqro9jatMXdGHPaGLMv/LwbaATKYx1sLrweFz2DI+w9ds7q\nKLMyEhql6XTQVi2Zcd6asfW1nTRTWMXPy0fO8thbp/nzuhWUFdrnr9ZENKOeu4hUAFcBr0zy8rUi\nckBEHheRVVHINmvXVZaSniq2bQ8c6exlcGTUVjdTx9WWjc8Utuf3XlknNGp4YHsD5UXZ3LlhudVx\nbC/i4i4iecDPgXuMMRevELUPWGqMuQL4NvDfU7zHnSKyV0T2dnTE7souLzONtcuKbVtg6sOjTVYn\n+LIDk3HiTGEVHz977SSNp4N8dlNNwm0paUcRFXcRSWessD9ijPnFxa8bY4LGmJ7w851Auoi8YxEI\nY8xDxpg1xpg1LpdrjtEvzetx03ymh1Pn+mJ6nljw+YNkpqWw3KYz8rweN8GBEcfMFFaxF+gf5mtP\nHmRtRTG3X1ZmdRxHiGS0jAAPA43GmK9PccyC8HGIyNrw+56NZtCZqvOM9X7tOGrG5w9QU1Zg29Xv\nnDZTWMXet39ziK6+IbZu0a3zoiWS6rEe+DBw04ShjptE5C4RuSt8zPuBehE5AHwL+ICxeOeGFa5c\nFhdn267vbozBZ6NlBybjpJnCKvYOd/TwwxeP8YF3LbZlKzJRTTun1xjzPHDJX6XGmO8A34lWqGgQ\nEbweN/+19xQDwyHb9PBOdvXTPTBiy5EyE3k9br7yRBNtgQEWFCb+wmfKOn+/o4Hs9FQ+/R6P1VEc\nxZ5/90fI63HTPxzi1aNdVkeJ2NsbYtv3yh2cM1NYxdaug+3sOtjBX91cRWleptVxHMXRxX3d8hIy\n01Js1Xf3+YOkpgieBfbe2X18prC2ZtRUhkOjfHFHA8tKc/nj6yqsjuM4ji7u2RmpXLuixFZ993p/\ngCp3nm3aSFMZnyn8fEsnQyM6JFK9049fOs6Rjl6+sHklGWmOLkWWcPx31Otxc6Szl2OdvVZHiYjP\nH6TW5i2ZcRdmCh+3T1tMxcfZnkG++XQzN1a78IZHtqnoSoriDvbYvLm9e4CO7kHb30wd9/ZMYW3N\nqN/2taea6R8K8YXNK3XoY4w4vrgvKclheWmuLTZvttOeqZG4MFNYlwBWEzT4g/z01RN8+NqlVLrt\nfW8pkTm+uMPYhKaXjpylfyhkdZRL8oWXHXBKWwbG/nI61G7PmcIq+owxbNvhozA7nXturrY6jqMl\nRXH31rgYGhnlpSOJvXmzzx9kaUkOBVnpVkeJGjvPFFbR90R9Gy8f6eLT7/FQmOOcn/NElBTFfe2y\nYrLTUxN+WN74Gu5OYteZwir6BoZDPLizkZoF+Xxw7RKr4zheUhT3zLRU1leWsutgOxavijClQP8w\nJ7r6HHMzddz4TOEXWs4yMJzYbTEVW99/7ginzvWzdUstqbp1XswlRXGHsdbMqXP9HO7osTrKpBoc\ndjN1IjvOFFbR1RYY4Lu7DrNx1QKuW/GOBWNVDCRNcR/v/SZqa+btZQecdeUOb88U1lUik9dXn2gi\nZAyf27TS6ihJI2mKe3lRNp75+exO0LVOGvxB3PmZuPKdt77G2zOFE/MXq4qtfSfO8Ys3Wvk/Nyxj\nSUmO1XGSRtIUd4C6GhevHu2iZ3DE6ijvUO8POHq5U6/HzdHOXo7aZKawio7R8NZ57vxM/qKu0uo4\nSSWpirvX42Y4ZHihJbGGRA4Mhzjc0evIfvu4Os/YKpE6aia5/PKNVg6cPM+9t9WQmzntCuMqipKq\nuF+zdB75mWkJV2Ca2roJjRpHF/elJbm2mSmsoqN3cISvPNHEFYuL+J0ry62Ok3SSqrinp6ZwfVUp\nu5o6EmpIpJNvpk5U53Hzsg1mCqvo+KfdLbR3D3L/llpSdOhj3CVVcYex1kxbcICmtm6ro1xQ3xqk\nMDudRfOyrY4SU3aZKazm7sTZPv71uaP83lXlXL1kntVxklLSFfcbw73fRBqW1+APUFtW4PjV8ewy\nU1jN3T/sbCQtRfjMxhqroyStpCvu8wuyWLWwgN0JUmBGQqM0tXU7ut8+zg4zhdXcvdjSyRO+Nu72\nVur+uRZKuuIOY62Z10+cI9A3bHUUDnf0Mjgy6uhhkBMl+kxhNTcjoVG27Whg0bxsPn79MqvjJLXk\nLO41LkKjhudarL96r291xobYkUr0mcJqbn762kma2rq5b9NK228VaXdJWdyvXDyPopz0hCgwPn+Q\nrPQUlrvyrI4SF+MzhRPpnoeKjkDfMF978iDrlhezcfUCq+MkvaQs7qkpwoYqF882dzA6am3v1+cP\nULOgIKlWyaurcfHasS66B6xvi6no+eZvmgn0D7N18yrHDw6wg2mLu4gsFpFdItIgIj4R+eQkx4iI\nfEtEWkTkTRG5OjZxo8db46KzZ/DC1nZWGB01NPiDrC5PjpbMuLrq8ZnCZ62OoqKkpb2bH790nA+u\nXeKoncTsLJIr9xHg08aYWmAdcLeI1F50zG1AVfhxJ/DPUU0ZAxuqXIhYOyTy5Lk+ugdHHD956WJr\nKuaRl4AzhdXsjG2d10hORiqfulW3zksU0xZ3Y8xpY8y+8PNuoBG4eC7x+4AfmzEvA0UiUhb1tFFU\nkpfJ5YuKLC3uTtsQO1LpqSncUKVDIp1i18F29jR3cM8t1ZTkOW9VU7uaUc9dRCqAq4BXLnqpHDg5\n4eNTvPMXQMLxelzsP3mert4hS87v8wdISxGq5yffDvBej5szQWvbYmruhkZG+eKORla4cvnItUut\njqMmiLi4i0ge8HPgHmPMrP5HisidIrJXRPZ2dFg/UsXrcWMM7Gm2Jkt9a5BKd15SDhm7pXY+eZlp\nfPPpQ1ZHUXPwoxePcbSzl89vriU9NSnHZySsiP41RCSdscL+iDHmF5Mc0gosnvDxovDnfosx5iFj\nzBpjzBqXyzWbvFF1WXkhJbkZlrVmxjbETq5++7ji3Aw+cVMlTzeeseyXq5qbju5BvvWbQ9xU48Yb\nnr+gEkcko2UEeBhoNMZ8fYrDfgV8JDxqZh0QMMacjmLOmEhJEW70jA2JDMV5SGR7cIDOnsGk67dP\n9LH1FSwtyeGLOxoYDo1aHUfN0NeePEj/cIj7btet8xJRJFfu64EPAzeJyP7wY5OI3CUid4WP2Qkc\nAVqAfwX+IjZxo8/rcXO+b5j9J8/H9bzjveZkWXZgMplpqdy3aSWH2nt45OXjVsdRM1DfGuBne0/y\n0esqWJEkE/DsZtqtUYwxzwOXnJFgxoY83B2tUPG0ocpFioztEHTN0vgtTTq+7MDKsuS7mTrRrbXz\nub6ylG88fYj3XVnOvNwMqyOpaRhj2La9geKcDP7y5iqr46gpJP0dkMKcdK5ZOi/ufXefP0hFSQ75\nWelxPW+iERG+sLmWnsERvvF0s9VxVAQee+s0rx7r4m/e66EwO7l/fhNZ0hd3GFvMqr41SHv3QNzO\n6TsdSNqbqRfzLMjnQ+9ewr+/fJymNh0amcj6h0J8aWcTtWUF/OGaxdN/gbKMFne4cKf/2Tjt7xno\nG+ZkVz+rkmzZgUu555Zq8rPS2ba9QSc2JbCH9hyh9Xw/92+pTar1kOxIiztjfe/5BZnsjlNx951O\njj1TZ2JebgafurWaFw+f5cmGM1bHUZPwn+/nn59t4fbLynj38hKr46hpaHFnrO9bV+1mz6GOuAzJ\na0jSZQem87/fvYTq+Xk8+FgjgyO6iXai+coTTRgD996mW+fZgRb3MG+Ni+6BEfYdPxfzc/n8QRYU\nZFGq63D8lrTUFLZuXsWJrj5+8Pwxq+OoCfYe6+J/9vv5sw3LWVycY3UcFQEt7mHrK0tJSxF2xaE1\nU98a0Kv2KVxfVcqttfP5zjOHaA/G7wa3mtroqOGB7Q0sKMjirroVVsdREdLiHpaflc67Kopjvgxt\n/1CIwx09Wtwv4b5NKxkOGb7664NWR1HAo/tO8VZrgM9uqiEnY9qpMSpBaHGfwFvjoqmtG//5/pid\no6ktyKiBWr2ZOqWK0lw+dn0Fj75+igNxnjmsflv3wDBffeIgVy8p4o4rFlodR82AFvcJxodExnLU\nTP2FZQf0yv1SPuGtpDQvkwe2+3RopIW+u+swnT2D3L9Ft86zGy3uE1S68ygvyo7pbNUGf4DC7HTK\ni7Jjdg4nyM9K5zMbPew7cZ5fHfBbHScpHevs5QfPH+X91yziisVFVsdRM6TFfQIRwVvj4oWWzpgN\nxRtb5rdAr4Ii8P6rF3FZeSFf2tlE39CI1XGSzoM7G0lPFT7zXo/VUdQsaHG/iNfjpm8oxN5j0R8S\nORwapamtW2+mRiglRfi7O2ppCw7wvd2HrY6TVJ471MFTDWf4xE1VuAuyrI6jZkGL+0WuXVFCRloK\nu5qi35ppae9haGQ0qZf5nalrlhZzxxUL+Zc9RzjZ1Wd1nKQwEhpl2/YGlhTn8CfXV1gdR82SFveL\n5GSksW55SUz67sm6IfZc3XtbDSLw5cebrI6SFB555QSH2nu47/aVZKYl3xaQTqHFfRJ11S4Od/Ry\n4mx0rxR9/gDZ6aksK9XNDWZiYVE2f35jJY+9dZpXjpy1Oo6jnesd4utPNbO+soT31M63Oo6aAy3u\nk/DWhIdENkf36t3nD7KyLF9X05uFOzcsp7womwe2N8R9S8Rk8s2nm+keGGbrZh36aHda3CexrDSX\nipKcqPbdR0cNDUm8IfZcZWek8tlNNTScDvKfe09aHceRDrZ18++vnOBD65biWZDcO4Q5gRb3KdR5\n3Lx4+CwDw9EZEnmiq4+ewRHtt8/B7ZeVsbaimP/764ME+oetjuMoxhi27fCRl5nGX99SbXUcFQVa\n3KfgrXEzODLKS1Hq8b59M1Wv3GdLRNi6pZauviG+/ZtDVsdxlKcazvBCy1n++pYq3cfWIbS4T+Hd\ny4rJSk9hd5RaM/X+AGkpQvUCvZk6F6vLC/mjNYv54YvHONzRY3UcRxgcCfHgzkaq3Hn873VLrY6j\nokSL+xQvv3C2AAAOy0lEQVSy0lNZv6KUXQc7orK2ic8fpGp+vg4ti4JPv8dDdnoqDz7WaHUUR/i3\nF45x/GwfW7fUkp6qJcEp9F/yEupq3Jzo6uNoZ++c3scYQ4Nf13CPFld+Jn91cxXPNLXHdB2gZNDe\nPcC3f3OIW1bO54Yql9VxVBRpcb+EuuqxH/a5buDR3j1IZ8+QFvco+uPrKlhWmssXdzTEZWtEp/rH\nJw4yFBrlvttXWh1FRdm0xV1EfiAi7SJSP8XrdSISEJH94cfW6Me0xuLiHKrceXPewKO+dWxDbF12\nIHoy0lL4/O0rOdLRy49fOm51HFs6cPI8//X6Kf5k/TKWleZaHUdFWSRX7j8ENk5zzHPGmCvDj21z\nj5U4vDVuXjnSRe/g7Fcl9PmDiMDKMr1yj6abatxsqHbxzaebOdszaHUcWzHG8MB2H6V5GXzipkqr\n46gYmLa4G2P2AF1xyJKQ6qpdDIVGefHw7IdE+vwBKkpyycvULcqiSUTYunklfUMhvv5Us9VxbOVX\nB/zsO3Gez7y3hvysdKvjqBiIVs/9WhE5ICKPi8iqKL1nQlhTUUxuRuqcbtyNr+Guoq/Snc9Hrl3K\nT149QUN4LoG6tL6hEb78eBOXlRfy/msWWR1HxUg0ivs+YKkx5grg28B/T3WgiNwpIntFZG9HR+y2\nsoumjLQUrq8qZXdT+6yGRJ7vG+LUuX6dvBRD99xcTWF2Ott26JZ8kfjes0c4HRhg65ZaUnSdI8ea\nc3E3xgSNMT3h5zuBdBEpneLYh4wxa4wxa1wu+wy78nrc+AMDNJ+Z+aSZBl3mN+YKc9L51Hs8vHyk\niyfq26yOk9BOnevjX549zJYrFvKuimKr46gYmnNxF5EFEl4+TkTWht/TUeuy1oU3zp5Na0bXcI+P\nD75rMTUL8nlwZ2PU1gNyoi893oTI2Br5ytkiGQr5E+AlwCMip0Tk4yJyl4jcFT7k/UC9iBwAvgV8\nwDjsb+MFhVmsLCuY1SqR9f4AZYVZlORlxiCZGpeWmsLWzbWcOtfPw88ftTpOQnrlyFkee/M0d924\nQjdoTwLTDt8wxnxwmte/A3wnaokSlNfj4l/2HCE4MEzBDEYX6M3U+LmuspSNqxbw3V0t/P7Vi1hQ\nqHt/jguNGh7Y3sDCwiz+bMMKq+OoONAZqhHy1rgJjRpeONQZ8df0D4U40tFDrd5MjZvPbVrJSMjw\n1Sd0S76J/mvvSRpOB7l300qyM3R9o2SgxT1CVy0uoiArbUZ998a2IKNG++3xtKQkhz+9YRm/eKOV\nfSfOWR0nIQQHhvnHXx9kzdJ5bLm8zOo4Kk60uEcoLTWFDdWuGa0S6dNlByzxF95K3PmZPLC9gVHd\nko9v/+YQXX1D3L9Ft85LJlrcZ6DO46aje/DCCJjp+PxBinLSWai937jKy0zjbzfWcODkef57f6vV\ncSx1uKOHf3vhGH94zWIuW6QXGclEi/sM3BheJTLShcTGb6bq1VL8/e5V5VyxuIgvP940p3WB7O7B\nxxrJSk/lb97rsTqKijMt7jPgys/k8kWFES0BPBwa5WBbN6v1ZqolUlKE+7fU0t49yD/tbrE6jiV2\nH2znmaZ2/vKmSlz5OhQ32Whxn6E6j5s3TpzjXO/QJY87dKaHodAotXoz1TJXL5nH715Vzr8+d5QT\nZ/usjhNXw6FRvrijgYqSHD66vsLqOMoCWtxnyOtxMWpgz6FLX737/GM3U3VNGWv97cYaUkX4h53J\ntSXf/3vpOIc7evn87bW6tWOS0uI+Q5cvKqI4N4Pd07RmfP4g2empugmCxRYUZnG3dwVP+Np48XDk\ncxTs7GzPIN94upkbqkq5eaXb6jjKIlrcZyg1Rbix2sWzzR2ELjHMzucPULuwgFRddc9yf3rDchbN\ny2bb9gZGkmBLvq8/1UzfUIitm2v1Zn4S0+I+C3UeF129Q7x56vykr4+OGhp02YGEkZWeyn2bVtLU\n1s1PXztpdZyYajwd5CevnuDD65ZSNT/f6jjKQlrcZ2FDlYsUYcrWzPGuPnqHQlrcE8jG1Qt497Ji\nvvbkQQJ9w1bHiQljDNu2N1CQnc49t1RZHUdZTIv7LMzLzeCqJfOmHO+uN1MTj4iwdUstgf5hvvkb\nZ27J92tfGy8dOcunb62mKCfD6jjKYlrcZ8nrcXHgVICO7nduzFzfGiQ9VajWP4sTyqqFhXxg7RL+\n30vHaWnvtjpOVA0Mh/j7xxrxzM/ng2uXWB1HJQAt7rM0voHHnuZ3tmZ8/gBV7nwy0vTbm2g+fWs1\n2RmpbNvR6Kgt+R5+/iinzvWzdUstaan6c6e0uM9abVkBrvzMd6wSaYzeTE1kJXmZ3HNLNXuaO+a0\n6XkiaQsM8N1dLbx31XzWV066w6VKQlrcZyklRairdrGnueO3htedCQ5ytndIV4JMYB+5dinLXbl8\ncUcjQyP2Hxr51SeaGAkZ7ttUa3UUlUC0uM+Bt8ZNcGCEN06+PSSyvnX8ZqpeuSeq9NQUvrC5lqOd\nvfzoxWNWx5mTN06c4xdvtPLxG5axpCTH6jgqgWhxn4Prq0pJTZHf2lvV5w8iAivLtLgnMq/Hjdfj\n4lu/OURnzztvitvBaHjrPFd+Jnd7K62OoxKMFvc5KMhKZ83Seb+1SqTPH2BZSS65mdNuT6ss9vnN\ntfQPh/jakwetjjIr/72/lf0nz/O3G2vI0583dREt7nPkrXHTeDpIW2AACK/hrv12W1jhyuOj11Xw\n09dOXmin2UXv4AhffryJKxYV8ntXlVsdRyUgLe5z5A0Pidx9sJ1zvUO0nu/XfruN/OXNVczLyWDb\n9gZbDY38p90ttHcPsnXLKlJ0/SI1CS3uc1Q9P4+FhVnsPthBw+mx7fe0uNtHYXY6f/MeD68e6+Kx\nt05bHSciJ7v6+NfnjvI7Vy7kmqXzrI6jEpQW9zkSEepq3Dzf0sn+8KgZXXbAXv7oXYtZWVbAl3Y2\nMTAcsjrOtP5hZyOpIvztbTVWR1EJbNriLiI/EJF2Eamf4nURkW+JSIuIvCkiV0c/ZmLzetz0DI7w\nyMvHWViYRXGuruthJ6nhLflaz/fz0J4jVse5pBcPd/J4fRt/UbeCssJsq+OoBBbJlfsPgY2XeP02\noCr8uBP457nHspfrVpSQkZqCPzBArV6129K65SXcflkZ/7S7Bf/5fqvjTGokNMq27Q2UF2XzfzYs\ntzqOSnDTFndjzB6g6xKHvA/4sRnzMlAkImXRCmgHuZlprF1WDGi/3c7uva2GUQNfeaLJ6iiT+ulr\nJ2lq6+Zzm1aSla5b56lLi8bg2HJg4g4Ip8Kfs8fdqSip87h4vqVTlx2wscXFOfzZhuV8+5mWsclo\nVge6yKlz/axdVsymyxZYHUXZQFxnPojInYy1bliyxFnLkv7+1Ys4HRjgel24ydb+vG4F3QMjtHcP\nWB3lHS4rL+STt1Tp1nkqIhLJ2F4RqQB2GGNWT/LavwC7jTE/CX98EKgzxlzyyn3NmjVm7969s8ms\nlFJJS0ReN8asme64aAyF/BXwkfComXVAYLrCrpRSKrambcuIyE+AOqBURE4B9wPpAMaY7wE7gU1A\nC9AHfCxWYZVSSkVm2uJujPngNK8b4O6oJVJKKTVnOkNVKaUcSIu7Uko5kBZ3pZRyIC3uSinlQFrc\nlVLKgSKaxBSTE4t0AMdn+eWlQGcU40RLouaCxM2muWZGc82ME3MtNca4pjvIsuI+FyKyN5IZWvGW\nqLkgcbNprpnRXDOTzLm0LaOUUg6kxV0ppRzIrsX9IasDTCFRc0HiZtNcM6O5ZiZpc9my566UUurS\n7HrlrpRS6hJsV9xFZKOIHAxvyH2v1Xlg+k3ErSIii0Vkl4g0iIhPRD5pdSYAEckSkVdF5EA41wNW\nZ5pIRFJF5A0R2WF1lnEickxE3hKR/SKSMBshiEiRiDwqIk0i0igi1yZAJk/4+zT+CIrIPVbnAhCR\nvw7/zNeLyE9EJCtm57JTW0ZEUoFm4FbGtvN7DfigMabB4lwbgB7G9pJ9x4YmVgnvZVtmjNknIvnA\n68DvJMD3S4BcY0yPiKQDzwOfDO/BazkR+RSwBigwxmy2Og+MFXdgjTEmocZsi8iPgOeMMd8XkQwg\nxxhz3upc48I1oxV4tzFmtvNqopWlnLGf9VpjTL+I/Cew0xjzw1icz25X7muBFmPMEWPMEPBTxjbo\ntlQEm4hbwhhz2hizL/y8G2hkbH9bS4U3U+8Jf5gefiTEVYaILAJuB75vdZZEJyKFwAbgYQBjzFAi\nFfawm4HDVhf2CdKAbBFJA3IAf6xOZLfiPtVm3Goa4a0SrwJesTbJmHDrYz/QDjxljEmIXMA3gc8A\no1YHuYgBnhSR18N7ESeCZUAH8G/hNtb3RSTX6lAX+QDwE6tDABhjWoH/C5wATjO2a92TsTqf3Yq7\nmgURyQN+DtxjjAlanQfAGBMyxlwJLALWiojl7SwR2Qy0G2NetzrLJK43xlwN3AbcHW4FWi0NuBr4\nZ2PMVUAvkBD3wQDCbaI7gP+yOguAiMxjrNOwDFgI5IrIh2J1PrsV91Zg8YSPF4U/p6YQ7mn/HHjE\nGPMLq/NcLPxn/C5go9VZgPXAHeH+9k+Bm0Tk362NNCZ81Ycxph34JWMtSqudAk5N+KvrUcaKfaK4\nDdhnjDljdZCwW4CjxpgOY8ww8AvgulidzG7F/TWgSkSWhX8rf4CxDbrVJMI3Lh8GGo0xX7c6zzgR\ncYlIUfh5NmM3yJusTQXGmM8aYxYZYyoY+9l6xhgTsyurSIlIbviGOOG2x3sAy0dmGWPagJMi4gl/\n6mbA0pv1F/kgCdKSCTsBrBORnPD/zZsZuw8WE9PuoZpIjDEjIvIJ4NdAKvADY4zP4liTbiJujHnY\n2lTA2JXoh4G3wv1tgM8ZY3ZamAmgDPhReCRDCvCfxpiEGXaYgOYDvxyrB6QB/2GMecLaSBf8JfBI\n+GLrCPAxi/MAF34J3gr8mdVZxhljXhGRR4F9wAjwBjGcqWqroZBKKaUiY7e2jFJKqQhocVdKKQfS\n4q6UUg6kxV0ppRxIi7tSSjmQFnellHIgLe5KKeVAWtyVUsqB/j8G9014gxx7pAAAAABJRU5ErkJg\ngg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10a3b0128>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(d[d.x == 4].groupby([d.y]).count()[\"N\"].index,d[d.x == 4].groupby([d.y]).count()[\"N\"].values)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* 平均と分散を調べる。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style>\n",
       "    .dataframe thead tr:only-child th {\n",
       "        text-align: right;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: left;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>N</th>\n",
       "      <th>y</th>\n",
       "      <th>x</th>\n",
       "      <th>id</th>\n",
       "      <th>x_bure</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>40</th>\n",
       "      <td>8</td>\n",
       "      <td>6</td>\n",
       "      <td>4</td>\n",
       "      <td>41</td>\n",
       "      <td>0.088943</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>41</th>\n",
       "      <td>8</td>\n",
       "      <td>2</td>\n",
       "      <td>4</td>\n",
       "      <td>42</td>\n",
       "      <td>-0.017723</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>42</th>\n",
       "      <td>8</td>\n",
       "      <td>0</td>\n",
       "      <td>4</td>\n",
       "      <td>43</td>\n",
       "      <td>-0.032386</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>43</th>\n",
       "      <td>8</td>\n",
       "      <td>5</td>\n",
       "      <td>4</td>\n",
       "      <td>44</td>\n",
       "      <td>-0.024742</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>44</th>\n",
       "      <td>8</td>\n",
       "      <td>6</td>\n",
       "      <td>4</td>\n",
       "      <td>45</td>\n",
       "      <td>-0.072645</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "    N  y  x  id    x_bure\n",
       "40  8  6  4  41  0.088943\n",
       "41  8  2  4  42 -0.017723\n",
       "42  8  0  4  43 -0.032386\n",
       "43  8  5  4  44 -0.024742\n",
       "44  8  6  4  45 -0.072645"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "d_4.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "count    20.000000\n",
       "mean      4.050000\n",
       "std       2.892367\n",
       "min       0.000000\n",
       "25%       2.000000\n",
       "50%       3.500000\n",
       "75%       7.000000\n",
       "max       8.000000\n",
       "Name: y, dtype: float64"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "d_4.y.describe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "8.3657894736842131"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "d_4.y.describe()[\"std\"]**2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# heikin = (d_4.seizon * d_4.samplesuu).sum() / d_4.samplesuu.sum()\n",
    "# heikin\n",
    "# (((d_4.seizon - heikin) ** 2) * d_4.samplesuu).sum()/d_4.samplesuu.sum()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* 生存確率の平均は$4.05/8=0.5$ぐらい。二項分布に従うのであれば、その分散は、$8\\times0.5\\times(1-0.5)=2$ぐらいとなるはず。上記だと分散は$8.37$となっており比較的大きくなっている。\n",
    "     * 二項分布の分散は$N \\times p \\times (1-p)$で求められる。\n",
    "* このデータは二項分布と呼ぶには「ばらつきが大きすぎる」ので、二項分布を使っては説明できないだろうということです。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "### 7.2.2. 観測されていない個体差がもたらす過分散"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* 単純な二項分布モデルからのずれの原因の一つとして、観測されていない個体差があげられる。\n",
    "* 個体差があれば過分散は生じうる。\n",
    "* 過分散は現象の状態を記述しているのではなく、統計モデリングしている人間の錯誤を表している。\n",
    "    * 「この個体たちは皆均質」という過程を設け、「二項分布で説明できう」と期待したのに、失敗した、ということ。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 7.2.3. 観測されていない個体差"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Q. 個体差とは？"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A.\n",
    "* データとしては定量化も識別もされていないけれど、「各個体（観測の単位）の何かに起因している様に見える差」\n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* 個体差をもたらしうる要因は少なくとも2種類ある\n",
    "* 生物的な要因\n",
    "    * 個体の遺伝子、年齢や過去に経験した履歴の相違\n",
    "* 非生物的な要因\n",
    "    * 生育環境の微妙な違い\n",
    "    * 非生物的な局所環境の影響のことを、この本では場所差、あるいはブロック差と呼ぶ。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Q. 個体差、場所差は識別できるか？"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A.\n",
    "* 例えば、1個体ずつ植木鉢に植えている場合、その個体差と植木鉢の差は識別できない。\n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* 個体に影響を与えている要因すべてを定量・特定することは不可能。個体差や場所差を原因不明のまま、統計モデルを作る必要がある。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 7.3. 一般化線形混合モデル"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Q. 一般化線形混合モデル（GLMM）とは？"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A.\n",
    "* 個体差や場所差の効果をGLMに組み込んだ統計モデル\n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 7.3.1.個体差を表すパラメータの追加"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* 架空植物の種子の生存確率$q_i$を表す式に、個体$i$の個体差を表すパラメーター$r_i$を追加する。\n",
    " $$\n",
    " logit(q_i) = \\beta_1+\\beta_2 x_i + r_i\n",
    " $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* 「観測されていない個体差などない」と仮定しているGLMは、種子数を調査した全個体で$r_i=0$と設定していることになる。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 7.3.2. 個体差のばらつきを表す確率分布"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Q. GLMMの特徴は？"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A.\n",
    "* 個体差を表すパラメーター{$r_1,r_3,\\cdots,r_{100}$}が何か確率分布に従っていると仮定するところ。\n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* ここではとりあえず、個体差$r_i$は平均0で標準偏差$s$の正規分布に従うと仮定して見る。正規分布である根拠は何もなく、どの様な確率分布に従うかわからないが、便利なのでそうする。\n",
    "* 各個体の$r_i$は個体間で相互に独立した確率変数であると仮定すると、確率密度関数は下記。\n",
    "$$\n",
    "p(r_i | s)= \\frac {1} {\\sqrt {2\\pi s^2}} \\exp \\left(- \\frac {r_{i}^2} {2s^2}\\right)\n",
    "$$\n",
    "* 確率密度$p(r_i|s)$は$r_i$の「出現しやすさ」を表している。\n",
    "    * $r_i$が$0$に近い個体はわりと「ありがち」。絶対値が大きな個体は「あまりいない」\n",
    "    \n",
    "* 標準偏差$s$は「集団内の$r_i$のばらつき」を表している。\n",
    "    * $s$が大きければ大きいほど、個体差は大きくなり、過分散がひどくなる。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 7.3.3. 線形予測子の構成要素：固定効果とランダム効果"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "Q. 一般化線形 **混合** モデル（GLMM)という名前の「混合」とは何か？"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A.\n",
    "* 線形予測子に、固定効果とランダム効果の表す項を持っているもの。\n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* 固定効果とランダム効果といった分類はわかりにくい。パラメーターにはデータのかなり広い範囲を説明する大域的なものと、ごく一部だけを担当している局所的なものがあると整理して統計モデルを作るのが良いかも。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 7.4. 一般化線形混合モデルの最尤推定"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* GLMMに含まれている個体差$r_i$は最尤推定できない。\n",
    "* 理由は、{$r_1,r_2,\\cdots,r_{100}$}の100個のパラメーターと、{$\\beta_1,\\beta_2,s$}の3個のパラメーター合わせて103個を100個のデータから推測するのはナンセンス。\n",
    "* この様な場合の対応策の1つが、**個体ごとの尤度**$L_i$の式の中で、$r_i$を積分すること。（それ以外の対応策は8章で扱う）\n",
    "$$\n",
    "L_i = \\int_{-\\infty}^{\\infty} p(y_i|\\beta_1,\\beta_2,r_i) p(r_i|s)dr_i\n",
    "$$\n",
    "* この様にすると、尤度から$r_i$は消えてしまう。\n",
    "    * 色々な$r_i$の値における尤度を評価し、その期待値の算出をしている。\n",
    "* 尤度$L_i$を評価するときに、二項分布$p(y_i|\\beta_1,\\beta_2,r_i)$と正規分布をかけて$r_i$で積分する操作は、「2種類の分布を混ぜていること」に相当。\n",
    "* 無限個の二項分布を混ぜることで、平均よりも分散の大きい過分散な確率分布を作り出せる。\n",
    "    * サンプルサイズは有限個だが、無限個の二項分布を混ぜる。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "* 確率分布を混ぜて新しい確率分布を作る、というのは、統計モデル作りの基本的な技法の1つ\n",
    "* 全データの尤度は$L_i$の100個体分の積なので、$$L(\\beta_1,\\beta_2,s)=\\prod_i L_i$$\n",
    "* 個体差$r_i$がどこにも入っていない対数尤度$\\log L(\\beta_1,\\beta_2,s)$が最大になる様なパラメータ$\\beta_1,\\beta_2,s$の最尤推定値を探すことになる。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* （自分が理解してないので、推薦できるがわからないが、）線形混合モデルの可視化\n",
    "http://qiita.com/kenmatsu4/items/da59a6fccde0cee4ef49"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 7.4.1. Rを使ってGLMMのパラメーターを推定"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* PythonにはGLMMのライブラリがない模様。\n",
    "* Chap7_4_1を参照"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 7.5. 現実のデータ解析にはGLMMが必要"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Q. GLMMの様な考え方が必要になるかどうかの判断ポイントは？"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A. \n",
    "* 同じ個体・場所などからなんどもサンプリングしているか\n",
    "* 個体差や場所差が識別できてしまう様なデータのとり方をしているか"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 7.5.1. 反復・擬似反復と統計モデルの関係"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* 個体差・場所差をどの様に統計モデルに組み込むかの方法は、「データをどのように取ったのか」に依存している。\n",
    "* 個体差・場所差が「見えてしまう」データの取り方をしているのかどうかが重要。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 具体例\n",
    "* 植木鉢で育てられた植物から種子を取って生死を調べたとする。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "#### 各植木鉢に1つの個体、各個体で一つだけ種子の生死を調べた場合(図7.11(A))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* 過分散生じない。\n",
    "* 個体差のばらつきも推定できない。\n",
    "* 個体差と植木鉢の区別がつかない。\n",
    "* このようなデータの取り方をすれば、統計モデルに個体差・場所差を入れようがない（？）ので、GLMMを使えない。\n",
    "* GLMで推定して問題ない\n",
    "* （滝沢コメント）同じ$x_i$でも、種子の生存数のばらつきが大きい場合、それはその植物において、$x_i$が説明変数として妥当でない、ということになる（？）"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* このようなデータの取り方を「**反復**」、あるいは「独立した反復」と呼ぶ。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Q. 反復をとる目的は？"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A.\n",
    "* データから個体差・場所差の効果を除去し、GLMで推定するなど、統計モデルを簡単にすること。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 各植木鉢に1つの個体、各個体から複数の種子の生死を調べた場合(図7.11(B))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* 同じ個体から複数の種子を調べると、「この個体では種子がやたらと死んでいる」「こっちの個体は同じ$x_i$だけれど、生存確率が高いようだ」と言った事がデータから読み取れる。\n",
    "* 個体差を入れたGLMMを使えば、「個体差のばらつき」$s$が推定される。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* 「個体から複数のデータを取る」という実験操作は「**擬似反復**」。「反復」に見えるが、個体差などを打ち消す反復になっていない。\n",
    "* 個体差が推定可能。個体差の影響を考慮しなければ、推定結果に偏りが生じる。\n",
    "* GLMMが必要。\n",
    "* 植木鉢は反復になっていて、個体に差があるのか植木鉢に差があるのか区別ができないので、植木鉢差だけを分離して推定できない。（上記の「個体差」に「植木鉢差」も含まれている。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 各植木鉢に複数の個体、各個体から1つだけ種子の生死を調べた場合(図7.11(C))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* 個体は反復\n",
    "\n",
    "* 植木鉢には複数の個体が入っているため、「こっちの植木鉢では死亡確率が高い」といった事が推定可能。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 各植木鉢に複数の個体、各個体から複数の種子の生死を調べた場合(図7.11(D))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* 個体も植木鉢も擬似反復。\n",
    "* 種子の生存確率は\n",
    "$$\n",
    "logit(q_i) = \\beta_1+ \\beta_2 x_i + r_i + r_j\n",
    "$$\n",
    "* 個体差・場所差を考慮しつつ、$\\beta_1$,$\\beta_2$に加えて、個体差・場所差それぞれのばらつきの大きさを最尤推定するのは、数値計算の問題としては難しい。\n",
    "* このような時にはGLMMと最尤推定の組み合わせではなく、階層ベイズモデルとMCMCを使って統計モデルをデータに当てはめます(次章)\n",
    "* 擬似反復を正しく統計モデリングすれば、個体差・場所差を定量化できる。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 7.6. 色々な分布のGLMM"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* ポアソン分布\n",
    "    * 正規分布と混ぜ合わせる事が多い\n",
    "    * ガンマ分布と混ぜ合わされた、負の二項分布\n",
    "* 観測データのばらつきが正規分布やガンマ分布である時には過分散を定義することはできない。\n",
    "    * これらの確率分布は平均とは独立に分散を決められる\n",
    "        * 分散が大きければ、正規分布・ガンマ分布の中で、分散を大きくすれば良い（？）。\n",
    "* GLMMが必要になるかどうかの本質は、過分散の有無ではなく、「1個体から何度もデータを取る」ような擬似反復が含まれているかどうかで決まる。\n",
    "* データのばらつきが正規分布やガンマ分布である時にも、サンプリングが擬似反復になっている可能性がある場合には、GLMMのように、個体差・場所差を考慮したモデルを当てはめなければならない。\n",
    "* 分布は正規分布で恒等リンク関数とする統計モデルは線形混合モデルと呼ばれる事がある。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 7.7. この章のまとめと参考文献"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* 「観測されなかった個体差」を組み込んだGLMである、一般化線形混合モデル(GLMM)を扱った。\n",
    "* GLMがあてはまらない場合がある。\n",
    "* それは「説明変数が同じならどの個体も均質」と仮定していたGLMが誤り。個体差があり、実際には二項分布で期待されるより過分散。\n",
    "* GLMMは線形予測子に個体差のばらつき$r_i$を追加し、全個体の$r_i$がある確率分布に従うと仮定した統計モデル。\n",
    "* 積分によって$r_i$を消去した尤度を最大化する事で、GLMMの$\\beta$,$r_i$のばらつきといった、大域的なパラメータを最尤推定できる。\n",
    "* 擬似反復(一つの個体から複数データ、一つの場所に多数の個体）では、個体差・場所差を組み込んだGLMMが必要。\n",
    "* データのばらつきを表す確率分布に依らず、個体差・場所差に影響されるデータであれば、これらの効果をランダム効果として組み込んだ統計モデルで推定しなければならない。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* 「人間には観測できなかった・できない、しかし観測データにばらつきをもたらす」要因を組み込んだ統計モデリングは重要。\n",
    "* 個体差＋場所差の統計モデル（第10章）\n",
    "* 空間相関のある場所差を考慮した統計モデル(第11章）"
   ]
  }
 ],
 "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
}