Skip to content

Instantly share code, notes, and snippets.

@kiwamizamurai

kiwamizamurai/generalized_linear_model.ipynb

Created February 28, 2019 01:29
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save kiwamizamurai/8d49670d148a89f7e280358a14eee6af to your computer and use it in GitHub Desktop.
Save kiwamizamurai/8d49670d148a89f7e280358a14eee6af to your computer and use it in GitHub Desktop.
Download ZIP
Display the source blob
Display the rendered blob
Raw
generalized_linear_model.ipynb
{
"nbformat": 4,
"nbformat_minor": 0,
"metadata": {
"colab": {
"name": "generalized_linear_model.ipynb",
"version": "0.3.2",
"provenance": [],
"collapsed_sections": []
},
"kernelspec": {
"name": "python3",
"display_name": "Python 3"
}
},
"cells": [
{
"metadata": {
"id": "VOyzvEaxiUhw",
"colab_type": "text"
},
"cell_type": "markdown",
"source": [
"# GLM をやってみる\n",
"- http://hosho.ees.hokudai.ac.jp/~kubo/ce/EesLecture2018.html"
]
},
{
"metadata": {
"id": "rTUrraEAtJNy",
"colab_type": "text"
},
"cell_type": "markdown",
"source": [
"## 先に結果と感想\n",
"\n",
"- Python3で実装した\n",
"- 実はGLMはPythonではあんまり使われているところを知らない\n",
"- 調べてみたがRだとライブラリが充実していてソースも多かった\n",
"- とはいえ上の北大のLinkでRで結果と同様のものが得られたからOK"
]
},
{
"metadata": {
"id": "l7xl3IjoUZkD",
"colab_type": "code",
"colab": {}
},
"cell_type": "code",
"source": [
"import pandas as pd\n",
"df = pd.read_csv('http://hosho.ees.hokudai.ac.jp/~kubo/stat/2018/Fig/poisson/data3a.csv')"
],
"execution_count": 0,
"outputs": []
},
{
"metadata": {
"id": "2jugZMl4byQM",
"colab_type": "text"
},
"cell_type": "markdown",
"source": [
"- y: 種子数\n",
"- x: 体サイズ\n",
"- T = 施肥処理 \n",
"- C = 無処理"
]
},
{
"metadata": {
"id": "OXtyUESiU2Jl",
"colab_type": "code",
"outputId": "d250a353-767a-4761-e1bb-202709568696",
"colab": {
"base_uri": "https://localhost:8080/",
"height": 204
}
},
"cell_type": "code",
"source": [
"df.columns = ['seeds_num', 'body_size', 'status']\n",
"df.head()"
],
"execution_count": 135,
"outputs": [
{
"output_type": "execute_result",
"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>seeds_num</th>\n",
" <th>body_size</th>\n",
" <th>status</th>\n",
" </tr>\n",
" </thead>\n",
" <tbody>\n",
" <tr>\n",
" <th>0</th>\n",
" <td>6</td>\n",
" <td>8.31</td>\n",
" <td>C</td>\n",
" </tr>\n",
" <tr>\n",
" <th>1</th>\n",
" <td>6</td>\n",
" <td>9.44</td>\n",
" <td>C</td>\n",
" </tr>\n",
" <tr>\n",
" <th>2</th>\n",
" <td>6</td>\n",
" <td>9.50</td>\n",
" <td>C</td>\n",
" </tr>\n",
" <tr>\n",
" <th>3</th>\n",
" <td>12</td>\n",
" <td>9.07</td>\n",
" <td>C</td>\n",
" </tr>\n",
" <tr>\n",
" <th>4</th>\n",
" <td>10</td>\n",
" <td>10.16</td>\n",
" <td>C</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"</div>"
],
"text/plain": [
" seeds_num body_size status\n",
"0 6 8.31 C\n",
"1 6 9.44 C\n",
"2 6 9.50 C\n",
"3 12 9.07 C\n",
"4 10 10.16 C"
]
},
"metadata": {
"tags": []
},
"execution_count": 135
}
]
},
{
"metadata": {
"id": "OQ0keUKIVYN6",
"colab_type": "code",
"outputId": "b4404aa0-96ac-4a2a-ea08-c05d1830e2c0",
"colab": {
"base_uri": "https://localhost:8080/",
"height": 297
}
},
"cell_type": "code",
"source": [
"df.describe()"
],
"execution_count": 136,
"outputs": [
{
"output_type": "execute_result",
"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>seeds_num</th>\n",
" <th>body_size</th>\n",
" </tr>\n",
" </thead>\n",
" <tbody>\n",
" <tr>\n",
" <th>count</th>\n",
" <td>100.000000</td>\n",
" <td>100.000000</td>\n",
" </tr>\n",
" <tr>\n",
" <th>mean</th>\n",
" <td>7.830000</td>\n",
" <td>10.089100</td>\n",
" </tr>\n",
" <tr>\n",
" <th>std</th>\n",
" <td>2.624881</td>\n",
" <td>1.008049</td>\n",
" </tr>\n",
" <tr>\n",
" <th>min</th>\n",
" <td>2.000000</td>\n",
" <td>7.190000</td>\n",
" </tr>\n",
" <tr>\n",
" <th>25%</th>\n",
" <td>6.000000</td>\n",
" <td>9.427500</td>\n",
" </tr>\n",
" <tr>\n",
" <th>50%</th>\n",
" <td>8.000000</td>\n",
" <td>10.155000</td>\n",
" </tr>\n",
" <tr>\n",
" <th>75%</th>\n",
" <td>10.000000</td>\n",
" <td>10.685000</td>\n",
" </tr>\n",
" <tr>\n",
" <th>max</th>\n",
" <td>15.000000</td>\n",
" <td>12.400000</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"</div>"
],
"text/plain": [
" seeds_num body_size\n",
"count 100.000000 100.000000\n",
"mean 7.830000 10.089100\n",
"std 2.624881 1.008049\n",
"min 2.000000 7.190000\n",
"25% 6.000000 9.427500\n",
"50% 8.000000 10.155000\n",
"75% 10.000000 10.685000\n",
"max 15.000000 12.400000"
]
},
"metadata": {
"tags": []
},
"execution_count": 136
}
]
},
{
"metadata": {
"id": "oaYBWGtUVbwx",
"colab_type": "code",
"outputId": "42af9bc0-8933-491c-f22f-bfeb1c9f2d7c",
"colab": {
"base_uri": "https://localhost:8080/",
"height": 34
}
},
"cell_type": "code",
"source": [
"import collections\n",
"\n",
"collections.Counter(df[\"status\"])"
],
"execution_count": 137,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"Counter({'C': 50, 'T': 50})"
]
},
"metadata": {
"tags": []
},
"execution_count": 137
}
]
},
{
"metadata": {
"id": "FGWtIrT1U9NZ",
"colab_type": "code",
"outputId": "b72317dc-f883-4f17-8442-739d3968ea23",
"colab": {
"base_uri": "https://localhost:8080/",
"height": 187
}
},
"cell_type": "code",
"source": [
"import numpy as np\n",
"\n",
"df[\"label\"] = 0\n",
"\n",
"# Convert status into categorical\n",
"\n",
"def to_categorial(data):\n",
" for i in range(100):\n",
" if data[i] == 'T':\n",
" df[\"label\"][i] = 1\n",
" else:\n",
" df[\"label\"][i] = 0\n",
"\n",
" \n",
"if __name__ == '__main__':\n",
" to_categorial(df[\"status\"]) "
],
"execution_count": 138,
"outputs": [
{
"output_type": "stream",
"text": [
"/usr/local/lib/python3.6/dist-packages/ipykernel_launcher.py:10: SettingWithCopyWarning: \n",
"A value is trying to be set on a copy of a slice from a DataFrame\n",
"\n",
"See the caveats in the documentation: http://pandas.pydata.org/pandas-docs/stable/indexing.html#indexing-view-versus-copy\n",
" # Remove the CWD from sys.path while we load stuff.\n",
"/usr/local/lib/python3.6/dist-packages/ipykernel_launcher.py:8: SettingWithCopyWarning: \n",
"A value is trying to be set on a copy of a slice from a DataFrame\n",
"\n",
"See the caveats in the documentation: http://pandas.pydata.org/pandas-docs/stable/indexing.html#indexing-view-versus-copy\n",
" \n"
],
"name": "stderr"
}
]
},
{
"metadata": {
"id": "eRaSLhIgVzuV",
"colab_type": "code",
"colab": {}
},
"cell_type": "code",
"source": [
"%matplotlib inline\n",
"import matplotlib.pyplot as plt\n",
"plt.style.use('seaborn-darkgrid')\n",
"import seaborn as sns\n",
"from statsmodels.genmod.generalized_estimating_equations import GEE\n",
"from statsmodels.genmod.cov_struct import (Exchangeable, Independence,Autoregressive)\n",
"from statsmodels.genmod.families import Poisson\n",
"\n",
"plt.rcParams['figure.figsize'] = 14, 6"
],
"execution_count": 0,
"outputs": []
},
{
"metadata": {
"id": "Uw6J8jQuafEL",
"colab_type": "code",
"outputId": "604d008a-c962-45b2-827e-99b1e9129091",
"colab": {
"base_uri": "https://localhost:8080/",
"height": 153
}
},
"cell_type": "code",
"source": [
"df.groupby(['status']).mean().unstack()"
],
"execution_count": 140,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
" status\n",
"seeds_num C 7.7800\n",
" T 7.8800\n",
"body_size C 9.8076\n",
" T 10.3706\n",
"label C 0.0000\n",
" T 1.0000\n",
"dtype: float64"
]
},
"metadata": {
"tags": []
},
"execution_count": 140
}
]
},
{
"metadata": {
"id": "zPyFM3aKcSaJ",
"colab_type": "code",
"colab": {}
},
"cell_type": "code",
"source": [
"mask_T = df[\"label\"] == 1\n",
"mask_C = df[\"label\"] == 0"
],
"execution_count": 0,
"outputs": []
},
{
"metadata": {
"id": "FAlQIMNmalAr",
"colab_type": "code",
"outputId": "480127bc-2d8b-4666-a2b5-781ee49bc8ed",
"colab": {
"base_uri": "https://localhost:8080/",
"height": 410
}
},
"cell_type": "code",
"source": [
"plt.plot(df[df.columns[1]][mask_T], df[df.columns[0]][mask_T], 'o', color='r', label='T', alpha=0.5)\n",
"plt.plot(df[df.columns[1]][mask_C], df[df.columns[0]][mask_C], 'o', color='g', label='C', alpha=0.5)\n",
"\n",
"plt.xlabel(\"{}\".format(df.columns[1]), size=17)\n",
"plt.ylabel(\"{}\".format(df.columns[0]), size=17)\n",
"plt.legend()"
],
"execution_count": 142,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"<matplotlib.legend.Legend at 0x7fbbacb0a0b8>"
]
},
"metadata": {
"tags": []
},
"execution_count": 142
},
{
"output_type": "display_data",
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAA0AAAAF4CAYAAAB91ZmrAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3X10k/X9//FX0pTehIIlDVgBq9Mp\nRUTr15tSbzYRcHPTqRNwCsp3N0cPougURKdnTCYouIni3BQ35lC3HUAr9atzKs7hanE/dSq2DIGp\niKVNQ4Ha9DbJ7w/X2NLctkmuJNfzcQ7nNJ9cN+/rc72vT/rmuvKpxe/3+wUAAAAAJmA1OgAAAAAA\nSBYKIAAAAACmQQEEAAAAwDQogAAAAACYBgUQAAAAANOwGR3AQLhcLRGXKSzMV3OzJwnRIF2RIwiH\n/EA45AfCIT8QCTmSHE5nQdD2jL0DZLNlGR0CUhw5gnDID4RDfiAc8gORkCPGytgCCAAAAAAORQEE\nAAAAwDQogAAAAACYBgUQAAAAANOgAAIAAABgGhRAAAAAAEyDAggAAACAaVAAAQAAADANm9EBAAAA\nAMChVq26X//+d5327XOrvb1dRxwxWsOGDdfSpSsGtV0KIAAAAACDllVXq+yaalndTfI5itRVXiFv\n6fgBb+/662+SJD3/fJV27dqpefNujEucFEAAAMPF+0MTAJBcWXW1yqmqDLy2uhqVU1WpDinlxnO+\nAwQAMFTPh6bV1Sj5fIEPzay6WqNDAwBEKbumOqZ2I1EAAQAMlU4fmgCA4KzuppjajUQBBAAwVDp9\naAIAgvM5imJqNxIFEADAUOn0oQkACK6rvCKmdiNRAAEADJVOH5oAgOC8pePVceHF8jlHSlarfM6R\n6rjw4pSbAEFiFjgAgMG8pePVITELHACkOW/p+ISM3RdccGFct0cBBAAwXKI+NAEAOBSPwAEAAAAw\nDQogAAAAAKZBAQQAAADANCiAAAAAAJgGBRAAAAAA02AWOAAAAAApaffuT/Tgg7/Q/v3N8np9OvHE\nibruuhs1ZMiQAW+TAggAAADAoNW5a1VTXy13W5MceUUqL65QqWPgf+LA6/XqjjsW6sYbF6is7H/k\n9/u1cuUKrVmzWtdcc92At0sBBAAAAGBQ6ty1qtpZGXjt8jQGXg+0CPrnP7foyCOPUlnZ/0iSLBaL\n5s69QRbL4L7Fw3eAAAAAAAxKTX11TO3R+OSTj/TVrx7Xpy0nJ3dQj79JFEAAAAAABsnd1hRTe3Qs\n8vl8g1g/OAogAAAAAIPiyCuKqT0aJSVHqbb2gz5tnZ2d2rVrx4C3KVEAAQAAABik8uKKmNqjcdpp\nZ6ihoV6vv/53SZLP59Ovf71Kr7zy0oC3KTEJAgAAAIBB6pnoIJ6zwFmtVv3iFw9p+fK7tWbNamVn\nZ+u0087Q//7vjwYVKwUQAAAAgEErdYwfVMETTFFRkZYvvz+u2+QROAAAAACmkfQCaPv27ZoyZYqe\neOKJPu2bN2/W8ccfn+xwAAAwRFZdrXLXPKb8++5R7prHlFVXa3RIAGAKSX0EzuPxaMmSJZo0aVKf\n9o6ODj366KNyOp3JDAcAAENk1dUqp+rLPxhodTUqp6pSHZK8pfF9fAQA0FdS7wANGTJEq1ev1siR\nI/u0/+Y3v9EVV1wx6D9qBABAOsiuCf6HAUO1AwDiJ6l3gGw2m2y2vrv8z3/+o23btmn+/PlasWJF\nVNspLMyXzZYVcTmns2BAccI8yBGEQ34gnEHlR3uLZM8J2l5A3mUExg9EQo4Yx/BZ4JYtW6Y77rgj\npnWamz0Rl3E6C+RytQw0LJgAOYJwyA+EM9j8yM0tkNXV2K/d5xypdvIu7TF+IBJyJDlCFZmGzgLX\n0NCgXbt26ZZbbtGMGTPU2NioWbNmGRkSAAAJ11Ue/A8DhmoHAMSPoXeARo0apZdffjnwevLkyf1m\nhwMAINN4S8erQ19858fqbpLPUaSu8gomQACAJEhqAbR161bde++92rNnj2w2m1588UWtWrVKhx12\nWDLDAADAcN7S8RQ8AGCApBZAEyZM0Nq1a0O+v2nTpiRGAwAAAMBsDP0OEAAAAAAkEwUQAAAAANOg\nAAIAAABgGhRAAAAAAEyDAggAAACAaVAAAQAAADANCiAAAAAApkEBBAAAAMA0KIAAAAAAmAYFEAAA\nAADToAACAAAAYBoUQAAAAABMgwIIAAAAgGlQAAEAAAAwDZvRAQAABi+rrlbZNdWyupvkcxSpq7xC\n3tLxRocFABmBMTazUAABQJrLqqtVTlVl4LXV1aicqkp1SHxAA8AgMcZmHh6BA4A0l11THVM7ACB6\njLGZhwIIANKc1d0UUzsAIHqMsZmHAggA0pzPURRTOwAgeoyxmYcCCADSXFd5RUztAIDoMcZmHiZB\nAIA05y0drw6JGYoAIAEYYzMPBRAAZABv6Xg+jAEgQRhjMwuPwAEAAAAwDQogAAAAAKZBAQQAAADA\nNCiAAAAAAJgGBRAAAAAA06AAAgAAAGAaFEAAAAAATIMCCAAAAIBpUAABAAAAMA0KIAAAAACmQQEE\nAAAAwDQogAAAAACYBgUQAAAAANOgAAIAAABgGkkvgLZv364pU6boiSeekCTV19drzpw5mjVrlubM\nmSOXy5XskAAAAACYhC2ZO/N4PFqyZIkmTZoUaFu5cqVmzJihCy64QE8++aTWrFmjhQsXJjMsAIhJ\nnbtWNfXVcrc1yZFXpPLiCpU6xhsdFoAkyqqrVXZNtazuJvkcReoqr5C3NL3GAcYy+sCskloADRky\nRKtXr9bq1asDbT/96U+Vk5MjSSosLNQHH3yQzJAAICZ17lpV7awMvHZ5GgOv+dAEzCGrrlY5VV+O\nA1ZXo3KqKtUhpU0RxFhGH5hZUh+Bs9lsys3N7dOWn5+vrKwseb1ePfXUU7rwwguTGRIAxKSmvjqm\ndgCZJ7sm+PUeqj0VMZbRB2aW1DtAoXi9Xi1cuFDl5eV9Ho8LpbAwXzZbVsTlnM6CeISHDEaOIJxg\n+dFubZHdntO/3dJCPpkM59vE2lukIOOA2ltU8N+8SPX8YCwzvg/M0s+pKCUKoNtuu00lJSWaN29e\nVMs3N3siLuN0FsjlahlsaMhg5AjCCZUfub4CuTyN/ZfPH0k+mQjjh7nl5hbI6uo/DvicI9XuakmL\n/GAsM7YP0iFHMkGoItPwabA3btyo7Oxs3XDDDUaHAgARlRdXxNQOIPN0lQe/3kO1pyLGMvrAzJJ6\nB2jr1q269957tWfPHtlsNr344otyu93KycnR7NmzJUnHHHOMFi9enMywACBqPV+MZdYgwLy8pePV\nIaX1LHCMZfSBmVn8fr/f6CBiFc0tQ24tIhJyBOGQHwiH/EA45AciIUeSI2UfgQMAAACAZKEAAgAA\nAGAaFEAAAAAATIMCCAAAAIBpUAABAAAAMA0KIAAAAACmQQEEAAAAwDQogAAAAACYBgUQAAAAANOg\nAAIAAABgGhRAAAAAAEyDAggAAACAaVAAAQAAADANCiAAAAAApkEBBAAAAMA0bEYHAABInjp3rSp3\nbNDWpvckSROKJuriY7+rUsd4gyPrL6uuVtk11bK6m+RzFKmrvELe0tSLM1PVuWtVU18td1uTHHlF\nKi+uSMk8CaYnd+r21eofw5rVOHqECo8sTatjiDeuJ+BLFEAAYBJ17lqt2bpa/963LdBW81m19rXv\n0/cn/CilfjHMqqtVTlVl4LXV1aicqkp1SPzSlgR17lpV7fyy/12exsDrVMqTYHpy5wOrS5W2bVKr\nZNm+V275VOVplJT6xxBvXE9AXzwCBwAmUVNfrU9bdvdr39OyWzX11QZEFFp2TfB4QrUjvkLlQ6rl\nSTA9OVJt7Zvr1k+/eJ0OxxBvXE9AXxRAAGAS7rYmebo9/do93R6525oMiCg0qzt4PKHaEV+h8iHV\n8iSYnhxpsvTNdUvbF6/T4RjijesJ6IsCCABMwpFXpHxbfr/2fFu+HHlFBkQUms8RPJ5Q7YivUPmQ\nankSTE+OFPn75ro/74vX6XAM8cb1BPRFAQQAJlFeXKExBWP7tY8uGKvy4goDIgqtqzx4PKHaEV+h\n8iHV8iSYnhyp8PXNdd+YL16nwzHEG9cT0BeTIACASZQ6xut/J/woLWaB85aOV4fErFUG6cmHdJwF\nrid3SmuqpX2WwCxwDhPPAsf1BPRl8fv9fqODiJXL1RJxGaezIKrlYF7kCMIhPxAO+YFwyA9EQo4k\nh9NZELSdR+AAAAAAmAYFEAAAAADToAACAAAAYBoUQAAAAABMgwIIAAAAgGlQAAEAAAAwDQogAAAA\nAKZBAQQAAADANCiAAAAAAJgGBRAAAAAA07DFsnBHR4c+/vhjHTx4UH6/v9/7p512WtwCAwAAAIB4\ni7oA+utf/6o77rhDLS0t/d7z+/2yWCyqq6uLa3AAAAAAEE9RF0D33XefSktLdeWVV2r48OGyWCyJ\njAsAAAAA4i7qAsjlcunRRx/VUUcdlcBwAAAAACBxoi6ATjjhBO3du3fQBdD27ds1d+5czZkzR7Nm\nzVJ9fb0WLlwor9crp9OpFStWaMiQIYPaBwBzqHPXqqa+Wu62JjnyilReXKFSx/io1qvcsUFbm96T\nJE0omqiLj/1uVOvGS1ZdrbJrqmV1N8nnKFJXeYW8pcnbfzoZ6HlOBdHEnuh8TET/Rdom+Y1wyA8Y\nLWvx4sWLo1nwpJNOChQnWVlZam9vV0tLS59/BQUFYbfh8Xi0YMECnXjiiSoqKtLEiRO1dOlSffvb\n39aiRYtUV1enTz75RCeeeGKE7XRGjNduz4lqOZgXOZLe6ty1qtpZKU9Xq/zyy9PVqu3N2zQi1yFn\nvjPsemu2rta7je+o09upTm+nPm3ZrU9aPtHYgiMD6yYyP7LqapVTVSmLp1Xy+2XxtMq2fZt8Ixzy\nO0PHbkYDPc+JFk1+RBN7tPk4UInov0jbJL/5fAmH/PgCOZIcdntO0Paop8FuaGjQ7t27tWDBAn37\n29/Weeed1+9fJEOGDNHq1as1cuTIQNuWLVsC65577rl64403og0JgInV1FfH1N77/U9bdvdr39Oy\nO+K68ZJdE3w/odrNbKDnORVEE3ui8zER/Rdpm+Q3wiE/kAqifgRuyZIlKigo0NVXXz3gSRBsNpts\ntr67bGtrCzzy5nA45HK5Im6nsDBfNltWxOWczvB3pAByJH21W1uC/s9Ou6Ul7Hltt7ao29KpIUP6\njkXd6lS7te+6CcuP9hYp2P9KtbeogJzsY6DnORki7T+a2GPJx4FIRP9F3Cb5LYnPl5DIjwByxDhR\nF0CfffaZ1q1bp+OOOy5hwQT720LBNDd7Ii7jdBbI5eo/ZTfQgxxJb7m+Ark8jf3anfkjw57XXF+B\nbP4hau1q7dNuz7Z/sc3/rpvI/MjNLZDV1T92n3Ok2snJPgZ6nhMtmvyIJvZo83GgEtF/kbZJfvP5\nEg758QVyJDlCFZlRPwJ39NFHq6urK24B9cjPz1d7e7ukLx6z6/14HACEUl5cEVN77/fHFIzt1z66\nYGzEdeOlqzz4fkK1m9lAz3MqiCb2ROdjIvov0jbJb4RDfiAVRD0Jwle+8hU99NBDKi0tlcPhGNRO\n33zzTeXl5WnixInasWOH2traNG7cOK1Zs0annHKKTjjhhLDrMwkC4oEcSW/OfKdG5DrU3LFP7d1t\nKsp36rwjp0ac3cqZ79SYgiP1edfn2t/RrCFZQ3TKqFP1vXGz+qybyPzwO53yjXDI2rxPlvY2+Yqc\n6jxvKrMgBTHQ85xo0eRHNLFHm48DlYj+i7RN8pvPl3DIjy+QI8kRahIEiz/K586mTZumgwcP6sCB\nA8rOzlZubm7fDVks2rJlS9htbN26Vffee6/27Nkjm82mUaNG6b777tOiRYvU0dGhI444QsuWLVN2\ndnbY7URzy5Bbi4iEHEE45AfCIT8QDvmBSMiR5Aj1CFzU3wE65ZRTBjTxQW8TJkzQ2rVr+7WvWbNm\nUNsFAAAAgGhEXQDdc889iYwDAAAAABIu6gLo888/j7jM0KFDBxUMAAAAACRS1AXQqaeeGvERuLq6\nukEHBAAAAACJEnUBdN111/UrgFpbW/XWW2/p4MGDuvrqq+MeHAAAAADEU9QF0PXXXx/yveXLl8vl\ncsUlIAAAAABIlKj/EGo406dP17p16+KxKQAAAABImLgUQPv3749qkgQAAAAAMFLUj8D94Q9/6Nfm\n9/u1b98+VVVV6dhjj41rYAAAAAAQb1EXQEuXLg353tFHH63FixfHIx4AAAAASJioC6BXXnmlX5vF\nYtGwYcP4+z8AAAAA0kLUBdDo0aMTGQcAAAAAJFzUBZAkvfTSS3rnnXe0f/9++f3+Pu9ZLJawj8kB\nAAAAgNGiLoB++ctf6tFHH5XFYpHdbu/3R1EpgAAAAACkuqgLoGeffVY/+MEPNH/+fA0ZMiSRMQEA\nAABAQkRdAB08eFDf+973KH4AxFWdu1Y19dVytzXJeaBLZ38sTWjOls9RpK7yCnlLxxsdIv4rq65W\n2TXVsrqbUv78DCTW3rnoyCtSeXGFJjQqpu0E20apI3IfxRpvtPsJtpykoOsGi2HryODLxsNA+wow\nSjqNgQgva3GU81e/9tprOuGEEzRmzJgEhxSZx9MZcRm7PSeq5WBe5Ijx6ty1qtpZKU9XqyxNjeqo\ne0f/7twthy9Po1ol2/Zt8o1wyO90Jj028qOvrLpa5VRVyuJplfx+WTythp6fcAYSa+9c9MsvT1er\ndny4WYdXv61Rreq3nbySMf3yI9g2tjdv04hch5z5ofso1nij3U+w5f7x2Wb9q/FtSeqzrvOz/Rrz\n4t/7xPDvHZv1bOfbas1WTMcTjYH2Vbpg/Mg88R4DyZHksNtzgrZbo93AbbfdpgcffFBbt26NW1AA\nzK2mvjrws/XT3YGfq61f/pxdUy0YL9R5SMXzM5BYe+diD+unu/vkYqTtBNtGuPZI2xvsfoItt6dl\ntz5t6X9Mb/7r6X5t1dbdfa7LSPuPxUD7CjBKOo2BiCzqR+BuvfVWHTx4UNOnT1d2drZyc3P7vG+x\nWLRly5a4Bwggc7nbmgI/W9o8gZ+bLF/+bHU3CcYLdR5S8fwMJNbeudjD0uZRk6X/sqG2E2wb4doj\nbW+w+wm2nKfb069NktytjZK+2qetyeKRpS36/cdioH0FGCWdxkBEFnUBdNJJJ/Wb+Q0ABsORVySX\np1GS5M/L/+LRAklF/vzAMj5HkSGxoS+fo0hWV2PQ9lQzkFh752IPf16+ij4Pvv1ot9HTHs94o91P\nsOXybfkKxmEfKR3yNE6RP18NecH3P1gD7SvAKOk0BiKyqAuge+65J+qNPvfcc5o8ebLy84MPtAAg\nSeXFFaraWSlJ8o0Zq6zt2yRJFb6xgWW6yisMiQ19dZVXKKeqMmh7qhlIrL1zsYdvzFhV1AbffrTb\n6GmPZ7zR7ifYcqMLxirYf2WefvKl0qvv92mr8I3VhiBf+410PNEYaF8BRkmnMRCRRT0JQiwuv/xy\nfec739Hw4cPjvWlJTIKA+CBHjOfMd2pErkPNHfvUNsSqouGjNa1trCa0D5OvyKnO86YaNsMO+dGX\n3+mUb4RD1uZ9srS3GX5+whlIrL1zsb27TUX5Tk0ed7HGjf6foNsJlh/BtnHekVMjzmwWa7zR7ifY\nchcec7HKRv5Pv3XHHXtWvxgOm3yxDju+/7LxmKltoH2VLhg/Mk+8x0ByJDlCTYJg8fv9/njvrKys\nTBs3btTYsWMjLzwALldLxGWczoKoloN5kSMIh/xAOOQHwiE/EAk5khxOZ0HQ9qhngQMAAACAdEcB\nBAAAAMA0KIAAAAAAmAYFEAAAAADToAACAAAAYBoUQAAAAABMgwIIAAAAgGkMugDat29fv7Zrr702\nYX8EFQAAAAAGKuoCqKWlRXPnztWOHTskSbt27dKUKVN05pln6uKLL1ZDQ0Ng2WuuuUbDhg2Lf7QA\nAAAAMAhRF0ArVqzQRx99pKFDh0qSli5dqu7ubi1cuFAWi0WrVq1KWJAAAAAAEA+2aBd87bXXtGLF\nCh1++OE6cOCA3njjDa1YsUIXXHCBJk6cqAULFiQyTgAAAAAYtKjvAO3bt08lJSWSpDfffFMWi0Vf\n+9rXJEljxoyRy+VKTIQAAAAAECdRF0CFhYXau3evJGnTpk2aOHGi7Ha7JMnlcgV+BgAAAIBUFfUj\ncKeffrruuOMOnXrqqdq4caPuvPNOSZLH49Gjjz6qE088MWFBAgAAAEA8RH0HaMGCBRo+fLieeeYZ\nXXDBBZo+fbok6eWXX9bmzZs1f/78hAUJAAAAAPEQ9R2gUaNG6YknnujXfs4552jTpk0qLCwcUACt\nra269dZbdeDAAXV1dem6667T2WefPaBtAUg/de5a1dRXy93WJEdekcqLK1TqGG90WHGTVVer7Jpq\nWd1N8jmK1FVeIW9p5hwfkq93Tm0t7NLmEsk1PDsjr59UE2686n1edNQYZZ1wSkzXerLGwlQYk4yK\nIZH7zfTPskhSIa9iEXUBFMphhx02qPWfeeYZHX300br55pvV0NCgq6++Wn/5y18GGxaANFDnrlXV\nzsrAa5enMfA6Ez44supqlVP15fFZXY3KqapUh5TSHwxIXb1z6gOrS8+2bJO2SjpunFxFvoy6flJN\nuPFqQqP6XOtqaFDOruiv9WSNhakwJhkVQyL3m+mfZZGkQl7FKmwBNG7cOFkslqg3VldXF3MAhYWF\n+ve//y1JOnjw4IDvJAFIPzX11SHbM+FDI7sm+PFl11Sn7IcCUlvvnKq27g78bP10t3xFTkmZc/2k\nmnDjVdk/g68T7bWerLEwFcYko2JI5H4z/bMsklTIq1iFLYBmz57dpwB64YUXZLfbVVZWJrvdroMH\nD+rtt99Wd3e3LrvssgEF8K1vfUtPP/20pk6dqoMHD+qRRx6JuE5hYb5stqyIyzmdBQOKCeZBjhir\n3doiuz2nf7ulJSXOzaBjaG+Rghyf2ltUkALHh8ExJEd75dRBdWpIz8d4d6dy/tueKtdPpgk3XhW0\nq9+1brfnRH2tJ20sTIUxyagYErjfgZ6/jLlOUyGvYhS2APrJT34S+PmRRx7ROeeco7vvvrvPMn6/\nX4sWLZLNNrCn6Z599lkdccQR+u1vf6tt27bp9ttv19NPPx12neZmT8TtOp0FcrlaBhQTzIEcMV6u\nr0AuT2O/dmf+SMPPTTzyIze3QFZX/+PzOUeqndxLa0aNH71zaphtiBqtrZIkf75d3a0dX8SWAtdP\nJgo3XrXkqs+1brfnqLW1I+prPVljYSqMSUbFkMj9DuT8ZdLvIKmQV6GEKjKjngXuj3/8o+bMmdOv\n3WKx6Pvf/77++Mc/Diiwt99+W2eddZakLx65a2xslNfrHdC2AKSX8uKKmNrTTVd58OMI1Q5E0jt3\nKnxjAz/7xnz5c6ZcP6km3Hg12Gs9WWNhKoxJRsWQyP1m+mdZJKmQV7GK+rbNvn371NnZGfS97u5u\n7du3b0ABlJSU6N1339X555+vPXv2yG63Kysr8uNtANJfz7PRmTpzjrd0vDqktJoZB6mtd06d4LbK\nXzAiMAucM8Oun1QTbrzyOtTnWteoUeqIYRa4ZI2FqTAmGRVDIveb6Z9lkaRCXsXK4vf7/dEs+N3v\nfldWq1VLlizRuHHjAu21tbVasmSJOjo6Ij66Fkxra6tuv/12ud1udXd3a/78+Zo0aVLYdaK5ZZhJ\ntxaRGOQIwiE/EA75gXDID0RCjiRHqEfgor4DdOedd+qaa67RJZdcoqysLOXm5qqjo0Pd3d3Ky8uL\navKCYOx2ux544IEBrQsAAAAAsYi6ADr55JP10ksv6cUXX9SOHTv0+eefKy8vT8cee6ymTZumESNG\nJDJOAAAAABi0mKZuGzZsmKZPn56oWAAAAAAgoaKeBU6SXC6XVq5cqR/+8Ie66KKLVF9fr66uLj33\n3HOJig8AAAAA4ibqO0AffvihZs2apY6ODo0bN067du1Sd3e3du/erVtvvVVWq1UXXHBBImMFAAAA\ngEGJ+g7Q8uXLVVpaqk2bNulPf/qTsrOzJUlf+cpXdOONN+p3v/tdwoIEAAAAgHiIugB65513dMst\ntwSd7GDatGn68MMP4xoYAAAAAMRb1AWQxWKR3W4P+l5nZ6csFkvcggIAAACARIi6ADr22GO1Zs2a\noO9t3LhRxx9/fNyCAgAAAIBEiHoShKuvvlo33nijamtrdeaZZ8rr9erJJ5/Url279Prrr+vBBx9M\nZJwAAAAAMGhR3wH6xje+oQceeEAdHR165JFH1NnZqccff1wNDQ26//77NWXKlETGCQAAAACDFtMf\nQj3//PN1/vnnq7W1VZ9//rkKCgqUn5+fqNgAAAAAIK5i+kOokuTxeLRjxw69++67gTafzxfXoAAA\nAAAgEaIugPx+v+6//35NmjRJM2fO1Pz58+V2u9XQ0KCLLrpITU1NiYwTAAAAAAYt6gLo0Ucf1e9/\n/3tdeeWVeuyxx5STkyNJys3NVX5+vh544IGEBQkAAAAA8RD1d4A2bNign/zkJ5oxY4YkBf7uz/Dh\nw3XLLbfoxz/+cWIiBAAkVJ27VjX11XK3NcmRV6Ty4gqVOsYbHVbCZNXVKrumWlZ3k9TV9UVjdrZ8\njiJ1lVfIW2rMsffElbWtVtbmZvkKR8g7rjTmmHofn9HHNBjJOo5M6S8A0Yu6ANq7d68mTZoU9L3R\no0dr//79cQsKAJAcde5aVe2sDLx2eRoDrzOxCMqqq1VO1RfHZ3G5ZNu+TZLkPX6c5PMpp6pSHVLS\nfwHuiat3TNaGvbL4fbK6GqOOqffxSZLV1WjYMQ1Gso4jU/oLQGyifgRu1KhR2rFjR9D3duzYocLC\nwrgFBQBIjpr66pja0112zZfHlbVnd+Bn66e7gy6TLD377B2T9GVc0cYUajkjjmkwknUcmdJfAGIT\n9R2g008/XYsXL5bFYlFFRUWg/b333tM999yjyZMnJyRAAEDiuNuCT2ATqj3dWd1fHpfF4wn6c+9l\nkqVnn73j6P062phCLWfEMQ1Gso4jU/oLQGyivgN0yy23qLCwUNdee61OOukktbW16Zvf/KZmzJih\nvLw8vgMEAGnIkVcUU3u68zk/Oyq9AAAbxklEQVS+PC5/r79j1/vn3sskS88+/Yf8bb2e19HGFGo5\nI45pMJJ1HJnSXwBiE3UBVFhYqA0bNujBBx/UnDlzdNlll+mqq67SQw89pHXr1mn48OGJjBMAkADl\nxRUxtae7rvIvj8s7emzgZ9+YsUGXSZaeffaOSfoyrmhjCrWcEcc0GMk6jkzpLwCxifoRuM7OTq1Y\nsUJXXXWVpk2bJpfLpR//+Md66qmnNGnSJK1YsUJDhw5NZKwAgDjrmejALLPAeUvHq0NffMfDarWq\na8SIL94weBa43nHJahnwLHB9ji+NZzVL1nFkSn8BiI3F7/f7o1lw+fLleu6557R27VqVlJToxhtv\nVHV1tS666CJt2rRJU6dO1W233ZboeCVJLldLxGWczoKoloN5kSMIh/xAOOQHwiE/EAk5khxOZ0HQ\n9qgfgXvxxRd11113qaSkRB6PR6+88ooWLFigO+64Q0uWLNGmTZviFiwAAAAAJELUBZDL5VJpaakk\n6Z///Ke8Xq+mTp0qSTrmmGPU2NiYmAgBAAAAIE6iLoAKCgrU3NwsSdq8ebNKS0t12GGHSZIOHDig\nnJycxEQIAAAAAHES9SQIJ510kpYuXapzzjlH69ev17x58yRJfr9ff/zjH3X88ccnLEgAAAAAiIeo\n7wDdfPPN+vTTT3XfffdpwoQJmjVrliRp48aNWrduna699tqEBQkAAAAA8RD1HaBjjjlGmzZt0r59\n+zSiZ9pQSeXl5aqsrNRXv/rVhAQIAAAAAPESdQHUo3fxI0mjRo3SqFGj4hYQAAAAACRK1I/AAQAA\nAEC6owACAAAAYBoUQAAAAABMgwIIAAAAgGlQAAEAAAAwDQogAAAAAKZBAQQAAADANFKiANq4caMu\nuugiXXrppfrb3/5mdDgAAAAAMlTMfwg13pqbm/WrX/1KGzZskMfj0apVq/T1r3/d6LAAw9W5a1VT\nXy13W5MceUUqL65QqWN8wvaXVVer7JpqWd1N8jmK1FVeIW9p4vaXjpLZR4M9/+HW730cWwu7tLlE\ncg3PDrqfZOdhKkmVa+LQc3Bm2yhN3Now4LgGek6jzamemCSlRP/FQ6peB6mSo4mWqv0fL+l4HtMx\n5t6yFi9evNjIAF599VVZrVZNmzZNdrtdkydPjriOx9MZcRm7PSeq5WBeqZwjde5aVe2slKerVX75\n5elq1fbmbRqR65Az3xn3/WXV1SqnqlIWT6vk98viaZVt+zb5Rjjkd8Z/f+ng0PxIZh8N9vyHW//w\nj12B4/jA0qhnu99RW+Nu+fPy1JqtPvtJdh6mkkjnO1njx6HnoL3+I+341wsqavVrpC8/5jwc6DmN\nNqd6+mpI9WbZ3nlbFintx5SB9Fky8sMs43amjkM9OZKO5zGdYrbbc4K2G/4I3Keffqr29nZde+21\nuuKKK/TGG28YHRJguJr66pjaByu7Jvh2Q7WbUTL7aLDnP9z6veOttu4O/Gz9dHef5eIRRzpLlWvi\n0L7uOU+9z50UfVwDPafR5lQP66e7lbVnd7/2dBxTUvU6SJUcTbRU7f94ScfzmI4xH8rwR+Akaf/+\n/XrooYf02Wef6aqrrtKrr74qi8UScvnCwnzZbFkRt+t0FsQzTGSgVM2RdmtL0P+1aLe0JCbm9hYp\n2P+StLeoIEX7KBn69HUS+2iw5z/c+gXtChzHQXVqSM/HQHencv7b3rOfpOdhKonifCejD/qdg+5O\nKdumg+qUPadXe5R5ONBzGm1O9YlTCuRUrHGmkoH2WcLzwyTjdiaPQ05nQXqex3SM+RCGF0AOh0Nl\nZWWy2Ww68sgjZbfbtW/fPjkcjpDrNDd7Im7X6SyQy9USz1CRYVI5R3J9BXJ5Gvu1O/NHJiTm3NwC\nWV399+dzjlR7ivZRoh2aH8nso8Ge/3Drt+QqcBzDbEPUaG2VJPnz7epu7eizn2TnYSqJdL6TNX4c\neg5stiGyeFo10mdXa3dHv7hi3V6PSOc02pzqHaekQE7FGmcqGUifJSM/zDJuZ+o41JMj6Xge0ynm\nUEWy4Y/AnXXWWaqpqZHP51Nzc7M8Ho8KCwuNDgswVHlxRUztg9XzheVo280omX002PMfbv3e8Vb4\nxgZ+9o0Z22e5eMSRzlLlmji0r3vOU+9zJ0Uf10DPabQ51cM3Zqy8o8f2a0/HMSVVr4NUydFES9X+\nj5d0PI/pGPOhDJ8EYejQoWpvb9fPfvYzVVVV6eabb9YxxxwTdh0mQUA8pHKOOPOdGpHrUHPHPrV3\nt6ko36nzjpyasFlv/E6nfCMcsjbvk6W9Tb4ipzrPm5pWM7rE26H5kcw+Guz5D7d+7+MY1WbViOGj\n5TpurDwjhvXbT7LzMJVEOt/JGj8OPQeOoqN0XvG5mvB5/oDycKDnNNqc6omp48KL1V32Pxkxpgyk\nz5KRH2YZtzN1HOrJkXQ8j+kUc6hJECx+v9+f5FgGLZpbnqn8eBNSAzmCcMgPhEN+IBzyA5GQI8mR\nso/AAQAAAECyUAABAAAAMA0KIAAAAACmQQEEAAAAwDQogAAAAACYBgUQAAAAANOgAAIAAABgGhRA\nAAAAAEyDAggAAACAaVAAAQAAADANCiAAAAAApkEBBAAAAMA0KIAAAAAAmAYFEAAAAADTsBkdQLqq\nc9eqpr5a7rYmOfKKVF5coVLHeKPDAmIWKpfTIcez6mqVXVMtq7tJPkeRusor5C1Nj9jTVag+T2WJ\njDkZ/ZGOfQ4AqYwCaADq3LWq2lkZeO3yNAZe80sW0kmoXN61f6c+cL/fr11KnRzPqqtVTtWXsVtd\njcqpqtS/Pt+pKn9qx56uQvV5h5Syv5AnMuZk9Ec69jkApDoegRuAmvrqmNqBVBUqZ5/d+XRMyxsh\nuyZ4LG/+K/VjT1eh+jxUeypIZMzJ6I907HMASHUUQAPgbmuKqR1IVaFy1uVpjGl5I1jdIa7D1tSP\nPV2F6vNQ7akgkTEnoz/Ssc8BINVRAA2AI68opnYgVYXKWWf+yJiWN4LPEeI6tKd+7OkqVJ+Hak8F\niYw5Gf2Rjn0OAKmOAmgAyosrYmoHUlWonP3OMZfGtLwRusqDx3L6yakfe7oK1eeh2lNBImNORn+k\nY58DQKrLWrx48WKjg4iVx9MZcRm7PSeq5QbCme/UiFyHmjv2qb27TUX5Tp135FS+YJ1mEpkj6SJU\nLleMPivlc9zvdMo3wiFr8z5Z2tvkK3Kq87ypGnFSfGInP/oL1eep/GX8RMVst+fo8/zhCe+PdOxz\nMH4gMnIkOez2nKDtFr/f709yLIPmcrVEXMbpLIhqOZgXOYJwyA+EQ34gHPIDkZAjyeF0FgRt5xE4\nAAAAAKZBAQQAAADANCiAAAAAAJgGBRAAAAAA06AAAgAAAGAaFEAAAAAATIMCCAAAAIBpUAABAAAA\nMA0KIAAAAACmQQEEAAAAwDQogAAAAACYBgUQAAAAANOgAAIAAABgGhRAAAAAAEyDAggAAACAadiM\nDqBHe3u7vv3tb2vu3Lm69NJLjQ4HSAl17lrV1FfL3dYkR16RyosrVOoYb3RYQERZdbXKrqmW1d0k\nn6NIXeUV8pb2zd1E5TfXjXlEk2cAcKiUKYB+/etfa/jw4UaHAaSMOnetqnZWBl67PI2B1/wyh1SW\nVVernKovc9fqalROVaU6pMAvp4nKb64b84gmzwAgmJR4BG7nzp3asWOHvv71rxsdCpAyauqrY2oH\nUkV2TfAc7d2eqPzmujGPaPIMAIJJiTtA9957r+68805VVlZGXlhSYWG+bLasiMs5nQWDDQ0ZLpVz\npN3aIrs9p3+7pSWl484k9PMAtbdIQXJX7S0q+G+fJiq/k3ndkB8GiyLPjER+IBJyxDiGF0CVlZU6\n+eSTNXbs2KjXaW72RFzG6SyQy9UymNCQ4VI9R3J9BXJ5Gvu1O/NHpnTcmSLV8yOV5eYWyOrqn7s+\n50i1/7dPE5XfybpuyA/jRZNnRiE/EAk5khyhikzDH4H729/+pldeeUUzZszQunXr9PDDD6u6mtvX\nQHlxRUztQKroKg+eo73bE5XfXDfmEU2eAUAwht8BWrlyZeDnVatWafTo0aqoYPACer6wzWxWSDfe\n0vHqkMLOzpWo/Oa6MY9o8gwAgjG8AAIQWqljPL+4IS15S8dH/EU0UfnNdWMe0eQZABwqpQqg66+/\n3ugQAAAAAGQww78DBAAAAADJQgEEAAAAwDQogAAAAACYBgUQAAAAANOgAAIAAABgGhRAAAAAAEyD\nAggAAACAaVAAAQAAADANCiAAAAAApkEBBAAAAMA0KIAAAAAAmAYFEAAAAADToAACAAAAYBoUQAAA\nAABMw2Z0AIBRtjZu1fMfvCR3W5MceUUqL65QqWO80WEBSVXnrlVNfTXXARBEVl2tsmuqZXU3yeco\nUld5hbylXB9AuqMAginVuWu1ae8LavV0SJJcnkZV7ayUJH75g2nUuWsDeS9xHQC9ZdXVKqfqy+vD\n6mpUTlWlOiSKICDN8QgcTKmmvjqmdiATcR0AoWXXBL8OQrUDSB8UQDAld1tTTO1AJuI6AEKzuoNf\nB6HaAaQPCiCYkiOvKKZ2IBNxHQCh+RzBr4NQ7QDSBwUQTKm8uCKmdiATcR0AoXWVB78OQrUDSB9M\nggBTKnWM14gRdmaBg6n15DuzwAH9eUvHq0NiFjggA1EAwbQmjJygUZYSo8MADFXqGE/BA4TgLR1P\nwQNkIB6BAwAAAGAaFEAAAAAATIMCCAAAAIBpUAABAAAAMA0KIAAAAACmQQEEAAAAwDQogAAAAACY\nBgUQAAAAANOgAAIAAABgGhRAAAAAAEyDAggAAACAaVAAAQAAADANCiAAAAAApkEBBAAAAMA0KIAA\nAAAAmIbN6AB6LF++XG+99Za6u7t1zTXXaNq0aUaHBAAZq85dq5r6arnbmuTIK1J5cYVKHeONDgsA\ngIRLiQKopqZGH374of785z+rublZl1xyCQUQACRInbtWVTsrA69dnsbAa4ogAECmS4kC6LTTTtPE\niRMlScOGDVNbW5u8Xq+ysrIMjgwAMk9NfXXIdgogAECmS4kCKCsrS/n5+ZKk9evX65xzzglb/BQW\n5stmi1wcOZ0FcYsRmYkcQTiZmh/t1hbZ7Tn92y0tGXvMiUBfIRzyA5GQI8ZJiQKox8svv6z169fr\nd7/7Xdjlmps9EbfldBbI5WqJV2jIQOQIwsnk/Mj1FcjlaezX7swfmbHHHG+ZnB8YPPIDkZAjyRGq\nyEyZWeA2b96s3/zmN1q9erUKCqiIASBRyosrYmoHACCTpMQdoJaWFi1fvly///3vddhhhxkdDgBk\ntJ7v+TALHADAjFKiAHr++efV3NysG2+8MdB277336ogjjjAwKgDIXKWO8RQ8AABTSokCaObMmZo5\nc6bRYQAAAADIcCnzHSAAAAAASDQKIAAAAACmQQEEAAAAwDQogAAAAACYBgUQAAAAANOgAAIAAABg\nGhRAAAAAAEyDAggAAACAaVAAAQAAADANi9/v9xsdBAAAAAAkA3eAAAAAAJgGBRAAAAAA06AAAgAA\nAGAaFEAAAAAATIMCCAAAAIBpUAABAAAAMI2MK4DWrVun2bNnB/6VlZUZHRJSSGtrq+bNm6fZs2fr\n8ssv1+bNm40OCSnG5/Ppzjvv1OWXX67Zs2dr586dRoeEFLB9+3ZNmTJFTzzxhCSpvr5es2fP1hVX\nXKH58+ers7PT4AhhtENzRJL+8Ic/6IQTTlBra6uBkSEVBBtD5syZo1mzZmnOnDlyuVwGR2guGVcA\nTZ8+XWvXrtXatWt1/fXX6+KLLzY6JKSQZ555RkcffbTWrl2rBx54QHfffbfRISHFvPLKK2ppadGf\n/vQn3X333Vq+fLnRIcFgHo9HS5Ys0aRJkwJtDz74oK644go99dRTKikp0fr16w2MEEYLliOVlZVy\nu90aOXKkgZEhFQTLj5UrV2rGjBl64oknNHXqVK1Zs8bACM0n4wqg3n71q19p7ty5RoeBFFJYWKj9\n+/dLkg4ePKjCwkKDI0Kq+eijjzRx4kRJ0pFHHqnPPvtMXq/X4KhgpCFDhmj16tV9fpHdsmWLzjvv\nPEnSueeeqzfeeMOo8JACguXIlClTdNNNN8lisRgYGVJBsPz46U9/qvPPP19S399NkBwZWwC99957\nKi4ultPpNDoUpJBvfetb+uyzzzR16lTNmjVLt956q9EhIcUcd9xxev311+X1erVr1y7t3r1bzc3N\nRocFA9lsNuXm5vZpa2tr05AhQyRJDoeDx1dMLliODB061KBokGqC5Ud+fr6ysrLk9Xr11FNP6cIL\nLzQoOnPK2AJo/fr1uuSSS4wOAynm2Wef1RFHHKGXXnpJjz/+uO666y6jQ0KK+drXvqYTTzxRV155\npR5//HF95Stfkd/vNzospDDyA8BAeL1eLVy4UOXl5X0ej0Pi2YwOIFG2bNmiO+64w+gwkGLefvtt\nnXXWWZKkcePGqbGxUV6vV1lZWQZHhlRy0003BX6eMmWKHA6HgdEgFeXn56u9vV25ublqaGjgex4A\nYnbbbbeppKRE8+bNMzoU08nIO0ANDQ2y2+2BxxOAHiUlJXr33XclSXv27JHdbqf4QR/btm3Tbbfd\nJkn6+9//rvHjx8tqzcihEoNQUVGhF198UZL017/+VWeffbbBEQFIJxs3blR2drZuuOEGo0MxJYs/\nA+/db926VStXrtRjjz1mdChIMa2trbr99tvldrvV3d2t+fPnc9sZffh8Pt1+++3asWOHcnJydN99\n96m4uNjosGCgrVu36t5779WePXtks9k0atQo3XfffVq0aJE6Ojp0xBFHaNmyZcrOzjY6VBgkWI5U\nVFSourpa//rXv3TiiSfq5JNP1sKFC40OFQYIlh9ut1s5OTmB74odc8wxWrx4sbGBmkhGFkAAAAAA\nEAzPdQAAAAAwDQogAAAAAKZBAQQAAADANCiAAAAAAJgGBRAAAAAA06AAAgAAAGAaFEAAgJhMnjxZ\nc+fOTci26+rqdPzxx+vpp59OyPbDmT17tr7zne8kfb8AgOSyGR0AAACpYNWqVeJP4wFA5qMAAgBA\n0mGHHWZ0CACAJOAROADAgGzYsEFTp07VhAkTdP755+uFF14IvLdr1y7NmzdPp59+uiZMmKBp06bp\n4YcfltfrDSzT2tqqRYsW6dRTT1VZWZnmzp2rxsbGwPsvv/yyjj/+eP3jH//ot+9vfOMbuuaaa2KK\n97HHHtM3v/lNnXTSSTrttNN01VVX6Z///Gfg/d6PwG3ZskXHH3980H+LFi0KrPPpp5/qpptu0rnn\nnqsTTzxR3/rWt/TnP/85prgAAMnFHSAAQMxqa2vl9Xp1//33Kzs7W/fff79uvvlmHXPMMSoqKtKV\nV16p0aNH6+GHH5bT6dTmzZt1zz33qKWlRbfeeqskacmSJXrhhRd011136eSTT9Zbb72lpUuXBvZx\n7rnn6ogjjtCGDRt05pln9tn3f/7zHy1YsCDqeDds2KD77rtPP//5z1VeXq62tjatXbtW3//+9/XK\nK69o5MiRfZYvKyvT66+/3qdt3bp1WrVqlb75zW9Kkg4ePKgrr7xSdrtdS5cu1eGHH66//vWv+ulP\nfyq/36/LL7885n4FACQeBRAAIGYHDhzQL3/5S9ntdknSsmXLVFFRoaqqKhUUFKi5uVnr1q3TmDFj\nJEklJSXatm2b/vSnP+mmm26S1+vV//3f/2nmzJmBuy4lJSVqaGjQypUrJUlZWVmaMWOGHn74Ye3f\nvz/wiNpzzz2nUaNG6etf/3rU8b7//vvKz8/XxRdfLJvti4++O++8U5dcckngGHobMmSInE5n4PVb\nb72lhx9+WDfeeKO+9rWvSZLWr1+vvXv36tlnn9W4ceMkSddcc40++OADPfzww5o5c6YsFkss3QoA\nSAIegQMAxKy0tLRP4VBYWKji4mLt2LFD77//voqLiwPFT4+ysjJ5PB7t2rVLH3/8sTo7OzVhwoQ+\ny5xyyil9Xk+fPl1+v19VVVWSJL/fr+eff16XXXaZsrKyoo73vPPOU3t7u2bOnKknn3xSO3bsUHZ2\ntsrKyoIWQL3V19frhhtu0JQpU/o8dvfOO+/I6XQGip8eZ555phoaGuRyuaKODwCQPNwBAgDEbNiw\nYf3a7Ha72traJElDhw7t935PW2tra2C2tfz8/H7b6K2oqEjTpk3T+vXrNXv2bL311ltqbGzU9OnT\nY4r37LPP1pNPPqk//OEPevDBB7V//36NHTtWc+fO1aWXXhpyvfb2dl133XUqKirSsmXL+rzX0tKi\npqYmlZWV9Wnv+Z5TQ0NDv0frAADGowACAMSstbU1aFtJSYmsVqv+85//9Hu/paVF0hfFU2dnpyQF\nCqZDl+ntiiuu0JVXXqm6ujo999xzOvvss1VcXBxzzGVlZSorK5PP59N7772n1atX67bbbtPo0aN1\nxhlnBF3nJz/5ifbs2aMNGzYoLy+vz3vDhg3T4YcfrscffzzouqNGjYo5RgBA4vEIHAAgZh988IE+\n//zzwOumpiZ99tlnOu644zRx4kTt3btXn3zySZ913nrrLRUUFOioo45SSUmJbDab3nvvvT7L9J6V\nrcepp56q4447TpWVlXrhhRc0c+bMmOPdvHmz3n//fUmS1WrVySefrJ///OeSpHfffTfoOo8++qj+\n8pe/aOXKlf0e55O+eFzP5XIpNzdXJSUlgX95eXkqKChQbm5uzHECABKPAggAELOhQ4dq4cKFqqur\n07Zt23TrrbfKZrPpwgsv1He/+105nU7ddNNNeuedd/TRRx/pt7/9rTZu3Kgf/OAHys7O1tChQzV1\n6lQ9/fTTev755/Xxxx9r3bp1ev7554Pu73vf+56eeOIJ5ebmBiYhiMWGDRs0d+5cvfzyy9qzZ492\n7dql3/zmN7LZbDr99NP7Lf/aa6/p/vvv1w9+8AMde+yxcrlcgX/79u2TJF166aVyOp26/vrrtWXL\nFu3Zs0evvfaaZs2apZtvvjnmGAEAyWHx82evAQAxmDx5ssrKylReXq5HHnlEe/fu1dixY7VgwQJN\nnjxZkvTRRx9pxYoVevPNN9XW1qYjjzxSl19+ua666qrAdg4cOKCf/exnevXVVyVJZ5xxhn70ox/p\niiuu0LJly/p8N+fAgQM644wzNG/ePM2bNy/mmD0ej375y1/qlVdekcvlkt1u11e/+lX98Ic/DMwm\nN3v2bB08eFDPPvusFi1apGeeeSbotkaPHq1NmzZJkvbs2aNf/OIX+sc//qGWlhYVFRXpG9/4hm64\n4Yag34MCABiPAggAkPLWrVunn//853r11Vc1YsQIo8MBAKQxJkEAAKSsvXv36p133tGyZct07bXX\nUvwAAAaNO0AAgJR1wgknqLCwUJdddpmuv/76Pn/75//9v/+nH/3oRxG3ceGFF+quu+5KZJgAgDRC\nAQQASEvt7e1qaGiIuNzQoUPlcDiSEBEAIB1QAAEAAAAwDabBBgAAAGAaFEAAAAAATIMCCAAAAIBp\nUAABAAAAMI3/D51pMQ+0jSP8AAAAAElFTkSuQmCC\n",
"text/plain": [
"<Figure size 1008x432 with 1 Axes>"
]
},
"metadata": {
"tags": []
}
}
]
},
{
"metadata": {
"id": "FWEHq3i_igNi",
"colab_type": "text"
},
"cell_type": "markdown",
"source": [
"とりあえずPlotまで行った、次はGLMにかける前に公式のテストデータで予習する"
]
},
{
"metadata": {
"id": "tYSv_yBbh4L5",
"colab_type": "text"
},
"cell_type": "markdown",
"source": [
"# いきなりやる前に公式にそって使い方を学ぶ\n",
"- http://www.statsmodels.org/dev/glm.html\n",
"- https://qiita.com/TomokIshii/items/374ac7d4231adf6a39f4\n",
"\n",
"\n",
"### 説明変数psi, tuce, gpa から応答変数post gradeを予測する問題\n",
"\n",
"1. 学期の通知表(GPA)\n",
"2. 夏季全国模試の成績(TUCE)\n",
"3. 夏期講習への参加状況(PSI)\n",
"\n",
"- 学期の成績が上がる/下がる(post grade)\n",
"\n",
"を推定するという問題"
]
},
{
"metadata": {
"id": "7BUt_Cs6fpaS",
"colab_type": "code",
"colab": {}
},
"cell_type": "code",
"source": [
"import statsmodels.api as sm\n",
"# Official dataset\n",
"spector_data = sm.datasets.spector.load()\n",
"spector_data.exog = sm.add_constant(spector_data.exog)\n",
"# add_constantが切片？"
],
"execution_count": 0,
"outputs": []
},
{
"metadata": {
"id": "Zk8jt1l4jq12",
"colab_type": "code",
"outputId": "c9f45907-2107-47c4-8ddb-ca94d05e999c",
"colab": {
"base_uri": "https://localhost:8080/",
"height": 306
}
},
"cell_type": "code",
"source": [
"print(pd.DataFrame(spector_data.data).head())\n",
"print('\\n')\n",
"print(pd.DataFrame(spector_data.data).describe())"
],
"execution_count": 150,
"outputs": [
{
"output_type": "stream",
"text": [
" GPA TUCE PSI GRADE\n",
"0 2.66 20.0 0.0 0.0\n",
"1 2.89 22.0 0.0 0.0\n",
"2 3.28 24.0 0.0 0.0\n",
"3 2.92 12.0 0.0 0.0\n",
"4 4.00 21.0 0.0 1.0\n",
"\n",
"\n",
" GPA TUCE PSI GRADE\n",
"count 32.000000 32.000000 32.000000 32.000000\n",
"mean 3.117188 21.937500 0.437500 0.343750\n",
"std 0.466713 3.901509 0.504016 0.482559\n",
"min 2.060000 12.000000 0.000000 0.000000\n",
"25% 2.812500 19.750000 0.000000 0.000000\n",
"50% 3.065000 22.500000 0.000000 0.000000\n",
"75% 3.515000 25.000000 1.000000 1.000000\n",
"max 4.000000 29.000000 1.000000 1.000000\n"
],
"name": "stdout"
}
]
},
{
"metadata": {
"id": "nJfRbk-EkD8n",
"colab_type": "code",
"colab": {
"base_uri": "https://localhost:8080/",
"height": 170
},
"outputId": "b692b6f7-4dc9-4a25-caeb-fedefb2a5877"
},
"cell_type": "code",
"source": [
"print(pd.DataFrame(spector_data.exog).describe())"
],
"execution_count": 152,
"outputs": [
{
"output_type": "stream",
"text": [
" 0 1 2 3\n",
"count 32.0 32.000000 32.000000 32.000000\n",
"mean 1.0 3.117188 21.937500 0.437500\n",
"std 0.0 0.466713 3.901509 0.504016\n",
"min 1.0 2.060000 12.000000 0.000000\n",
"25% 1.0 2.812500 19.750000 0.000000\n",
"50% 1.0 3.065000 22.500000 0.000000\n",
"75% 1.0 3.515000 25.000000 1.000000\n",
"max 1.0 4.000000 29.000000 1.000000\n"
],
"name": "stdout"
}
]
},
{
"metadata": {
"id": "lARfO8BThJ12",
"colab_type": "code",
"outputId": "66ec09cc-51f4-4d35-8035-a8be7ecca319",
"colab": {
"base_uri": "https://localhost:8080/",
"height": 374
}
},
"cell_type": "code",
"source": [
"logit_mod = sm.Logit(spector_data.endog, spector_data.exog)\n",
"logit_res = logit_mod.fit(disp=0)\n",
"print('Parameters: ', logit_res.params)\n",
"print('\\n')\n",
"print(logit_res.summary2())\n",
"# https://stackoverflow.com/questions/49814258/statsmodel-attributeerror-module-scipy-stats-has-no-attribute-chisqprob"
],
"execution_count": 0,
"outputs": [
{
"output_type": "stream",
"text": [
"Parameters: [-13.02134686 2.82611259 0.09515766 2.37868766]\n",
"\n",
"\n",
" Results: Logit\n",
"==============================================================\n",
"Model: Logit No. Iterations: 7.0000 \n",
"Dependent Variable: y Pseudo R-squared: 0.374 \n",
"Date: 2019-02-28 00:42 AIC: 33.7793\n",
"No. Observations: 32 BIC: 39.6422\n",
"Df Model: 3 Log-Likelihood: -12.890\n",
"Df Residuals: 28 LL-Null: -20.592\n",
"Converged: 1.0000 Scale: 1.0000 \n",
"---------------------------------------------------------------\n",
" Coef. Std.Err. z P>|z| [0.025 0.975]\n",
"---------------------------------------------------------------\n",
"const -13.0213 4.9313 -2.6405 0.0083 -22.6866 -3.3561\n",
"x1 2.8261 1.2629 2.2377 0.0252 0.3508 5.3014\n",
"x2 0.0952 0.1416 0.6722 0.5014 -0.1823 0.3726\n",
"x3 2.3787 1.0646 2.2344 0.0255 0.2922 4.4652\n",
"==============================================================\n",
"\n"
],
"name": "stdout"
}
]
},
{
"metadata": {
"id": "0CgcSueNl4Ur",
"colab_type": "text"
},
"cell_type": "markdown",
"source": [
"上ではいきなりLogitを読んでいたので下の書き方が正攻法だろう"
]
},
{
"metadata": {
"id": "30a2XlmWiAcy",
"colab_type": "code",
"colab": {
"base_uri": "https://localhost:8080/",
"height": 391
},
"outputId": "e35591ca-4274-44a1-e256-bb90f4345344"
},
"cell_type": "code",
"source": [
"glm_model = sm.GLM(spector_data.endog, spector_data.exog, family=sm.families.Binomial())\n",
"glm_reslt = glm_model.fit()\n",
"print('Parameters: ', glm_reslt.params)\n",
"print('\\n')\n",
"print(glm_reslt.summary2())"
],
"execution_count": 156,
"outputs": [
{
"output_type": "stream",
"text": [
"Parameters: [-13.02134686 2.82611259 0.09515766 2.37868765]\n",
"\n",
"\n",
" Results: Generalized linear model\n",
"=============================================================\n",
"Model: GLM AIC: 33.7793 \n",
"Link Function: logit BIC: -71.2613\n",
"Dependent Variable: y Log-Likelihood: -12.890 \n",
"Date: 2019-02-28 00:53 LL-Null: -20.592 \n",
"No. Observations: 32 Deviance: 25.779 \n",
"Df Model: 3 Pearson chi2: 27.3 \n",
"Df Residuals: 28 Scale: 1.0000 \n",
"Method: IRLS \n",
"--------------------------------------------------------------\n",
" Coef. Std.Err. z P>|z| [0.025 0.975]\n",
"--------------------------------------------------------------\n",
"const -13.0213 4.9313 -2.6406 0.0083 -22.6865 -3.3562\n",
"x1 2.8261 1.2629 2.2377 0.0252 0.3508 5.3014\n",
"x2 0.0952 0.1416 0.6722 0.5014 -0.1823 0.3726\n",
"x3 2.3787 1.0646 2.2344 0.0255 0.2922 4.4652\n",
"=============================================================\n",
"\n"
],
"name": "stdout"
}
]
},
{
"metadata": {
"id": "J4RN9lRsmVVa",
"colab_type": "text"
},
"cell_type": "markdown",
"source": [
"ちなみに残差も見てみる"
]
},
{
"metadata": {
"id": "3aLW3KCziHNW",
"colab_type": "code",
"colab": {
"base_uri": "https://localhost:8080/",
"height": 289
},
"outputId": "aea896bb-2c80-477a-973c-44a60a957cee"
},
"cell_type": "code",
"source": [
"resid1 = logit_res.resid_dev\n",
"resid2 = glm_reslt.resid_deviance\n",
"\n",
"print('Residual of Logit: \\n', resid1)\n",
"print('\\n')\n",
"print('Residual of GLM: \\n', resid2)"
],
"execution_count": 159,
"outputs": [
{
"output_type": "stream",
"text": [
"Residual of Logit: \n",
" [-0.23211021 -0.35027122 -0.64396264 -0.22909819 1.06047795 -0.26638437\n",
" -0.23178275 -0.32537884 -0.48538752 0.85555565 -0.22259715 -0.64918082\n",
" -0.88199929 1.81326864 -0.94639849 -0.24758297 -0.3320177 -0.28054444\n",
" -1.33513084 0.91030269 -0.35592175 0.44718924 -0.74400503 -1.95507406\n",
" 0.59395382 1.20963752 0.95233204 -0.85678568 0.58707192 0.33529199\n",
" -1.22731092 2.09663887]\n",
"\n",
"\n",
"Residual of GLM: \n",
" [-0.23211021 -0.35027122 -0.64396264 -0.22909819 1.06047795 -0.26638437\n",
" -0.23178275 -0.32537884 -0.48538752 0.85555565 -0.22259715 -0.64918082\n",
" -0.88199929 1.81326864 -0.94639849 -0.24758297 -0.3320177 -0.28054444\n",
" -1.33513084 0.91030269 -0.35592175 0.44718924 -0.74400503 -1.95507406\n",
" 0.59395382 1.20963752 0.95233204 -0.85678568 0.58707192 0.33529199\n",
" -1.22731092 2.09663887]\n"
],
"name": "stdout"
}
]
},
{
"metadata": {
"id": "ZlWQ2vkNm4AQ",
"colab_type": "text"
},
"cell_type": "markdown",
"source": [
"ちょっと使い方が怪しいのでもう少し公式を読みつつ確認する\n",
"- http://www.statsmodels.org/stable/index.html\n",
"- http://www.statsmodels.org/stable/examples/notebooks/generated/glm.html"
]
},
{
"metadata": {
"id": "QTBhga_pmX_v",
"colab_type": "code",
"colab": {
"base_uri": "https://localhost:8080/",
"height": 680
},
"outputId": "0a5384d1-1abd-4558-baa1-d60b70a019a4"
},
"cell_type": "code",
"source": [
"from __future__ import print_function\n",
"\n",
"print(sm.datasets.star98.NOTE)"
],
"execution_count": 162,
"outputs": [
{
"output_type": "stream",
"text": [
"::\n",
"\n",
" Number of Observations - 303 (counties in California).\n",
"\n",
" Number of Variables - 13 and 8 interaction terms.\n",
"\n",
" Definition of variables names::\n",
"\n",
" NABOVE - Total number of students above the national median for the\n",
" math section.\n",
" NBELOW - Total number of students below the national median for the\n",
" math section.\n",
" LOWINC - Percentage of low income students\n",
" PERASIAN - Percentage of Asian student\n",
" PERBLACK - Percentage of black students\n",
" PERHISP - Percentage of Hispanic students\n",
" PERMINTE - Percentage of minority teachers\n",
" AVYRSEXP - Sum of teachers' years in educational service divided by the\n",
" number of teachers.\n",
" AVSALK - Total salary budget including benefits divided by the number\n",
" of full-time teachers (in thousands)\n",
" PERSPENK - Per-pupil spending (in thousands)\n",
" PTRATIO - Pupil-teacher ratio.\n",
" PCTAF - Percentage of students taking UC/CSU prep courses\n",
" PCTCHRT - Percentage of charter schools\n",
" PCTYRRND - Percentage of year-round schools\n",
"\n",
" The below variables are interaction terms of the variables defined\n",
" above.\n",
"\n",
" PERMINTE_AVYRSEXP\n",
" PEMINTE_AVSAL\n",
" AVYRSEXP_AVSAL\n",
" PERSPEN_PTRATIO\n",
" PERSPEN_PCTAF\n",
" PTRATIO_PCTAF\n",
" PERMINTE_AVTRSEXP_AVSAL\n",
" PERSPEN_PTRATIO_PCTAF\n",
"\n"
],
"name": "stdout"
}
]
},
{
"metadata": {
"id": "owdmIC1QnD4B",
"colab_type": "code",
"colab": {
"base_uri": "https://localhost:8080/",
"height": 612
},
"outputId": "f0f75d3d-521f-405c-a774-8f4007a8bee3"
},
"cell_type": "code",
"source": [
"data = sm.datasets.star98.load()\n",
"# Load the data and add a constant to the exogenous (independent) variables:\n",
"data.exog = sm.add_constant(data.exog, prepend=False)\n",
"\n",
"glm_binom = sm.GLM(data.endog, data.exog, family=sm.families.Binomial())\n",
"res = glm_binom.fit()\n",
"print(res.summary())"
],
"execution_count": 163,
"outputs": [
{
"output_type": "stream",
"text": [
" Generalized Linear Model Regression Results \n",
"==============================================================================\n",
"Dep. Variable: ['y1', 'y2'] No. Observations: 303\n",
"Model: GLM Df Residuals: 282\n",
"Model Family: Binomial Df Model: 20\n",
"Link Function: logit Scale: 1.0\n",
"Method: IRLS Log-Likelihood: -2998.6\n",
"Date: Thu, 28 Feb 2019 Deviance: 4078.8\n",
"Time: 01:04:13 Pearson chi2: 9.60\n",
"No. Iterations: 5 \n",
"==============================================================================\n",
" coef std err z P>|z| [0.025 0.975]\n",
"------------------------------------------------------------------------------\n",
"x1 -0.0168 0.000 -38.749 0.000 -0.018 -0.016\n",
"x2 0.0099 0.001 16.505 0.000 0.009 0.011\n",
"x3 -0.0187 0.001 -25.182 0.000 -0.020 -0.017\n",
"x4 -0.0142 0.000 -32.818 0.000 -0.015 -0.013\n",
"x5 0.2545 0.030 8.498 0.000 0.196 0.313\n",
"x6 0.2407 0.057 4.212 0.000 0.129 0.353\n",
"x7 0.0804 0.014 5.775 0.000 0.053 0.108\n",
"x8 -1.9522 0.317 -6.162 0.000 -2.573 -1.331\n",
"x9 -0.3341 0.061 -5.453 0.000 -0.454 -0.214\n",
"x10 -0.1690 0.033 -5.169 0.000 -0.233 -0.105\n",
"x11 0.0049 0.001 3.921 0.000 0.002 0.007\n",
"x12 -0.0036 0.000 -15.878 0.000 -0.004 -0.003\n",
"x13 -0.0141 0.002 -7.391 0.000 -0.018 -0.010\n",
"x14 -0.0040 0.000 -8.450 0.000 -0.005 -0.003\n",
"x15 -0.0039 0.001 -4.059 0.000 -0.006 -0.002\n",
"x16 0.0917 0.015 6.321 0.000 0.063 0.120\n",
"x17 0.0490 0.007 6.574 0.000 0.034 0.064\n",
"x18 0.0080 0.001 5.362 0.000 0.005 0.011\n",
"x19 0.0002 2.99e-05 7.428 0.000 0.000 0.000\n",
"x20 -0.0022 0.000 -6.445 0.000 -0.003 -0.002\n",
"const 2.9589 1.547 1.913 0.056 -0.073 5.990\n",
"==============================================================================\n"
],
"name": "stdout"
}
]
},
{
"metadata": {
"id": "CGbXEDdZt0oB",
"colab_type": "text"
},
"cell_type": "markdown",
"source": [
"### OK、理解できたでは始めよう\n",
"\n",
"# Poisson Regression"
]
},
{
"metadata": {
"id": "zI_AvaSZo4zC",
"colab_type": "code",
"colab": {
"base_uri": "https://localhost:8080/",
"height": 204
},
"outputId": "c0cef1b7-6ce2-4beb-abb7-7ac458b2eff0"
},
"cell_type": "code",
"source": [
"X = sm.add_constant(df.iloc[:, 1], prepend=False)\n",
"X.head()"
],
"execution_count": 165,
"outputs": [
{
"output_type": "execute_result",
"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>body_size</th>\n",
" <th>const</th>\n",
" </tr>\n",
" </thead>\n",
" <tbody>\n",
" <tr>\n",
" <th>0</th>\n",
" <td>8.31</td>\n",
" <td>1.0</td>\n",
" </tr>\n",
" <tr>\n",
" <th>1</th>\n",
" <td>9.44</td>\n",
" <td>1.0</td>\n",
" </tr>\n",
" <tr>\n",
" <th>2</th>\n",
" <td>9.50</td>\n",
" <td>1.0</td>\n",
" </tr>\n",
" <tr>\n",
" <th>3</th>\n",
" <td>9.07</td>\n",
" <td>1.0</td>\n",
" </tr>\n",
" <tr>\n",
" <th>4</th>\n",
" <td>10.16</td>\n",
" <td>1.0</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"</div>"
],
"text/plain": [
" body_size const\n",
"0 8.31 1.0\n",
"1 9.44 1.0\n",
"2 9.50 1.0\n",
"3 9.07 1.0\n",
"4 10.16 1.0"
]
},
"metadata": {
"tags": []
},
"execution_count": 165
}
]
},
{
"metadata": {
"id": "vVwm0ttZoWDX",
"colab_type": "code",
"colab": {
"base_uri": "https://localhost:8080/",
"height": 544
},
"outputId": "99b970fc-bf0f-46f6-e5ed-cee8fdf0ef18"
},
"cell_type": "code",
"source": [
"poi_log = sm.GLM(df.iloc[:, 0], X, family=sm.families.Poisson(sm.families.links.log))\n",
"poi_log_results = poi_log.fit()\n",
"poi_resid = glm_reslt.resid_deviance\n",
"\n",
"print('Residual of GLM: \\n', poi_resid)\n",
"print('\\n')\n",
"print('Parameters: \\n', poi_log_results.params)\n",
"print('\\n')\n",
"print(poi_log_results.summary())"
],
"execution_count": 169,
"outputs": [
{
"output_type": "stream",
"text": [
"Residual of GLM: \n",
" [-0.23211021 -0.35027122 -0.64396264 -0.22909819 1.06047795 -0.26638437\n",
" -0.23178275 -0.32537884 -0.48538752 0.85555565 -0.22259715 -0.64918082\n",
" -0.88199929 1.81326864 -0.94639849 -0.24758297 -0.3320177 -0.28054444\n",
" -1.33513084 0.91030269 -0.35592175 0.44718924 -0.74400503 -1.95507406\n",
" 0.59395382 1.20963752 0.95233204 -0.85678568 0.58707192 0.33529199\n",
" -1.22731092 2.09663887]\n",
"\n",
"\n",
"Parameters: \n",
" body_size 0.075662\n",
"const 1.291721\n",
"dtype: float64\n",
"\n",
"\n",
" Generalized Linear Model Regression Results \n",
"==============================================================================\n",
"Dep. Variable: seeds_num No. Observations: 100\n",
"Model: GLM Df Residuals: 98\n",
"Model Family: Poisson Df Model: 1\n",
"Link Function: log Scale: 1.0\n",
"Method: IRLS Log-Likelihood: -235.39\n",
"Date: Thu, 28 Feb 2019 Deviance: 84.993\n",
"Time: 01:09:54 Pearson chi2: 83.8\n",
"No. Iterations: 4 \n",
"==============================================================================\n",
" coef std err z P>|z| [0.025 0.975]\n",
"------------------------------------------------------------------------------\n",
"body_size 0.0757 0.036 2.125 0.034 0.006 0.145\n",
"const 1.2917 0.364 3.552 0.000 0.579 2.005\n",
"==============================================================================\n"
],
"name": "stdout"
}
]
},
{
"metadata": {
"id": "w10pA1rGp3Px",
"colab_type": "text"
},
"cell_type": "markdown",
"source": [
"おおお、ちゃんと同じ値になった\n",
"- http://hosho.ees.hokudai.ac.jp/~kubo/stat/2017/Ees/c/HOkubostat2017c.pdf"
]
},
{
"metadata": {
"id": "oggICrxIpJ6X",
"colab_type": "code",
"colab": {
"base_uri": "https://localhost:8080/",
"height": 410
},
"outputId": "9430309f-59de-4557-b801-ea3f7f4d79aa"
},
"cell_type": "code",
"source": [
"plt.plot(df[df.columns[1]][mask_T], df[df.columns[0]][mask_T], 'o', color='r', label='T', alpha=0.5)\n",
"plt.plot(df[df.columns[1]][mask_C], df[df.columns[0]][mask_C], 'o', color='g', label='C', alpha=0.5)\n",
"\n",
"# Regression result\n",
"# Do not forget to apply data to Exponential function\n",
"preict = np.exp(poi_log_results.params[1] * np.ones_like(df.iloc[:, 1]) + poi_log_results.params[0] * df.iloc[:, 1])\n",
"plt.plot(df.iloc[:, 1], preict, color='black', label='Poisson Regression')\n",
"\n",
"plt.xlabel(\"{}\".format(df.columns[1]), size=17)\n",
"plt.ylabel(\"{}\".format(df.columns[0]), size=17)\n",
"plt.legend()"
],
"execution_count": 175,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"<matplotlib.legend.Legend at 0x7fbbac963f98>"
]
},
"metadata": {
"tags": []
},
"execution_count": 175
},
{
"output_type": "display_data",
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAA0AAAAF4CAYAAAB91ZmrAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzs3X18U/Xd//F30pSW1oIlDRUQUFBo\nK6I43Wp1TkRwl9MNb3BOQd38TZ3DyaaCKG4ol3fAvHA6p6Cig7l5ebNKvZyI4A2uFp2AAi333nDT\nmzQUKL1vkt8f2NLStEnaJCfJeT0fDx82356c8znf8zmn+XC+5xuL1+v1CgAAAABMwGp0AAAAAAAQ\nKRRAAAAAAEyDAggAAACAaVAAAQAAADANCiAAAAAApmEzOoDucDqr/S6Tnp6iqqraCESDWEWOoCvk\nB7pCfqAr5Af8IUciw+FI89ket3eAbLYEo0NAlCNH0BXyA10hP9AV8gP+kCPGitsCCAAAAACORgEE\nAAAAwDQogAAAAACYBgUQAAAAANOgAAIAAABgGhRAAAAAAEyDAggAAACAaVAAAQAAADANm9EBAAAA\nAPHkiSf+R1u2lGjfPpfq6+s1cOAg9enTVw89NM/o0CAKIAAAAJhcQkmxEosKZXVVymPPUFNuntzZ\nOd1e3223/VaS9NZbBdq5c4emTp0WqlARAhRAAADDhfrDBwAEKqGkWEkF+a2vrc4KJRXkq0HiOhSn\neAYIAGColg8fVmeF5PG0fvhIKCk2OjQAJpBYVBhUO2IfBRAAwFB8+ABgJKurMqh2xD4KIACAofjw\nAcBIHntGUO2IfRRAAABD8eEDgJGacvOCakfsowACABiKDx8AjOTOzlHDpRPlcfSXrFZ5HP3VcOlE\nJkCIY8wCBwAwlDs7Rw0Ss8ABMIw7Oycs15yLL7405OtEz1EAAQAMF64PHwAAHI0hcAAAAABMgwII\nAAAAgGlQAAEAAAAwDQogAAAAAKZBAQQAAADANJgFDgAAAAiDXbu+0Z/+9Eft318lt9ujU08drV//\neprRYZkeBRAAAABMrcRVrKLSQrnqKmXvnaHcAXnKtvdsan63261Zs6Zr2rS7NGbMd+T1erVgwTwt\nXrxIs2bdHaLI0R0UQAAAADCtElexCnbkt7521la0vu5JEfTpp2s0ZMgJGjPmO5Iki8WiW2/9jSwW\nnkAxGkcAAAAAplVUWhhUe6C++eYrnXzyiHZtSUnJ6tWrV4/Wi56jAAIAAIBpueoqg2oPnEUej6eH\n60A4UAABAADAtOy9M4JqD9TQoSeouHhTu7bGxkbt3Lm9R+tFz1EAAQAAwLRyB+QF1R6os876nsrL\nS/XRRx9Kkjwej/7ylye0cuWKHq0XPcckCAAAADCtlokOQj0LnNVq1R//+KTmzn1QixcvUmJios46\n63v6+c9/GYqw0QMUQAAAADC1bHtOjwseXzIyMjR37v+EfL3oGYbAAQAAADCNiBdAW7du1YUXXqil\nS5e2a1+9erVGjhwZ6XAAADBEQkmxkhc/q5T5jyh58bNKKCk2OiQAMIWIDoGrra3VnDlzdPbZZ7dr\nb2ho0MKFC+VwOCIZDgAAhkgoKVZSwZEvXrQ6K5RUkK8GSe7s0A/DAQAcEdE7QL169dKiRYvUv3//\ndu1PP/20rrnmGr4YCgBgColFvr9gsbN2AEDoRPQOkM1mk83WfpNffvmlNm/erNtvv13z5s0LaD3p\n6Smy2RL8LudwpHUrTpgHOYKukB/oSo/yo75aSk3y2Z5G3sUFrh/whxwxjuGzwD388MOaNWtWUO+p\nqqr1u4zDkSans7q7YcEEyBF0hfxAV3qaH8nJabI6Kzq0exz9VU/exTyuH/CHHImMzopMQwug8vJy\n7dy5U3feeackqaKiQpMnT+4wQQIAAPGkKTev3TNAbdsBxL7S0r267rqrNXJkliSpsbFR1157vX7w\ng7E+l3e5KvXcc89o+vR7wx7blVdeqv79M2W1WuX1epWUlKx77vm9MjKMexb/D3+YqXvu+YOSkpIj\nsj1DC6DMzEy9++67ra8vuOACih8AQNxzZ+eoQYef+bG6KuWxZ6gpN48JEIA4MmTIUD355EJJ0sGD\nB/Tzn1+r3NyzfX7It9szIlL8tJg//09KSUmRJL31VoEWLfqLZs78fcS2f7T77384otuLaAG0ceNG\nPfroo9qzZ49sNpuWL1+uJ554Qscee2wkwwAAwHDu7BwKHsAk+vTpK7s9Qy6XSzabTTNm3K6amjpZ\nrVbdffd9slgsmjVrhp57bomWLn1BH3zwnqxWq8455/u67rpf+Gxbu/Y/WrjwKdlsNjkc/TVz5u/1\n7rvL9cUX67V/f5W++eZrXXPNFF1yycQuY8vJGaX/+79lkqTPP1+nZ575s2w2m/r3z9SMGbNksVj0\nwAP3qaysVKeeOlqrVr2rf/7zLU2depOGDRsuSbrllql66KH7VV1dLbfbrWnT7tJJJ50c8L5ceeWl\n+utfX9ahQ9V6+OEH1NTU1K5vHnxwtgYOHKTt27dpxIiRuvvu+3p0PCJaAI0aNUpLlizp9PerVq2K\nYDQAAACId7Nnz1KBjyGnPXHppRM1e/Z/B7x8aeleHTx4QP37Z2ru3Ad15ZVX6qyzvq/33ntXzz+/\nUDfeeHPrsv/4x1Ll57+thIQE5ee/1mnb/PkP63/+58/KzDxOjz32qFaseFsWi0U7dmzX008/r927\nd+kPf7jHbwH0/vsrNWLE4aF6CxbM0+OP/0V9+vTVU089rvfee1cpKalqbGzQwoUv6N//Xq3//d+/\nt7532LDhmjjxSr3wwrP63vfydOmlE/Xllzv1+OPztWDBUwHvS4tnn31al1zyE40bN6Fd32zZUqL7\n739I6en9dNllF6u6ulppad2fRMLwSRAAAACAePPNN19r6tSbJB3+KphZs+6XzWbTli0luvfeu+X1\nSmeccaZeeOHZdu87//xxmjbtVo0f/0NNmPBDn20HDx6QxWJRZuZxkg6vZ/36tRoxIkujRo1WQkKC\nHI7+qqk55DO2O+/8jaxWq/bu3aPRo0/X9On3aN8+l3bv3qV77rlLklRfX6++fQ+P0jr11NMkSWef\nfY4SEo7MxJydPUqStGHDF9q/v0rLl78lSWpoqA94X9rasqVEt9wytXWfWvpm0KDBstszJEkZGQ7V\n1ByiAAIAAAB8mT37v4O6WxMqbZ8Bas8ir9crSWpqapbF0v5rOe+8c6a+/vorrVq1QrfddrMWLnyx\nQ9sf//hk6zoOr6epdT1tC5S2y7TV8gzQa6+9rF27diklJVXNzW5lZDg6xLx06QuyWg+v02KxyGKx\ntP4uMdHW+v/f/vYujRo1Ouh9WbjwRb9903afutqvQEX0i1ABAAAAM8vOztGaNWskSevXf6asrOzW\n3x06dEiLFy/S0KEn6Oc//6XS0vqqstLZoS0hwSqLxaKysrJv17O23XoC9ZOfXKF16z7Ttm1b1adP\nH0nSl1/ulCS9+uo/tH37Ng0adLy2bCmWJH3ySZHcbneH9eTkjNKHH77f+v5//GNpwPtSW1vTrm/W\nrv2Pz74JJe4AAQAAABHy//7fLfrjHx/S3/72d9lsiZo58z41NzdLko455hjt31+lX/7yOvXunaJR\no0bruOMGdGjr06evpk+fpfvvv1cJCQkaNOh4jRs3Qe+886+gYrHZbPr1r2/XY489oqeeek533/17\nPfTQ/UpMTFRGhkM//vHlGjx4iP7v/5bpV7+6UWPGfEd9+vTtsJ4rr/ypHnxwtm699f/J4/Fo2rQ7\ng9qXtn3z8MNzVFCQ36FvQsni7ek9JAME8sVRfMEU/CFH0BXyA10hP9AV8gP+xFKOHDx4QGvX/kfn\nnz9OTmeFbr/9V3rppdf8vzEKROUXoQIAAACIXikpqVq16l299NISeb0e3Xbb74wOqccogAAAAAD4\nZLPZ9MADkf2i0nBjEgQAAAAApkEBBAAAAMA0KIAAAAAAmAYFEAAAAADToAACAAAAYBoUQAAAAABM\ngwIIAAAAgGnwPUAAEAcSSoqVWFQoq6tSHnuGmnLz5M7OMTosAIgLXGPjCwUQAMS4hJJiJRXkt762\nOiuUVJCvBok/0ADQQ1xj4w9D4AAgxiUWFQbVDgAIHNfY+EMBBAAxzuqqDKodABA4rrHxhwIIAGKc\nx54RVDsAIHBcY+MPBRAAxLim3Lyg2gEAgeMaG3+YBAEAYpw7O0cNEjMUAUAYcI2NPxRAABAH3Nk5\n/DEGgDDhGhtfGAIHAAAAwDQogAAAAACYBgUQAAAAANOgAAIAAABgGhRAAAAAAEyDAggAAACAaVAA\nAQAAADANCiAAAAAApkEBBAAAAMA0KIAAAAAAmAYFEAAAAADToAACAAAAYBoUQAAAAABMgwIIAAAA\ngGlEvADaunWrLrzwQi1dulSSVFpaqhtuuEGTJ0/WDTfcIKfTGemQAAAAAJiELZIbq62t1Zw5c3T2\n2We3ti1YsEBXXXWVLr74Yv3tb3/T4sWLNX369EiGBQBBKXEVq6i0UK66Stl7Zyh3QJ6y7TlGhwUg\nghJKipVYVCirq1Iee4aacvPkzo6t6wDXMvrArCJaAPXq1UuLFi3SokWLWtv+8Ic/KCkpSZKUnp6u\nTZs2RTIkAAhKiatYBTvyW187aytaX/NHEzCHhJJiJRUcuQ5YnRVKKshXgxQzRRDXMvrAzCI6BM5m\nsyk5ObldW0pKihISEuR2u/XSSy/p0ksvjWRIABCUotLCoNoBxJ/EIt/ne2ft0YhrGX1gZhG9A9QZ\nt9ut6dOnKzc3t93wuM6kp6fIZkvwu5zDkRaK8BDHyBF0xVd+1FurlZqa1LHdUk0+mQzH28TqqyUf\n1wHVVyvt27yI9vzgWmZ8H5iln6NRVBRAM2fO1NChQzV16tSAlq+qqvW7jMORJqezuqehIY6RI+hK\nZ/mR7EmTs7ai4/Ip/cknE+H6YW7JyWmyOjteBzyO/qp3VsdEfnAtM7YPYiFH4kFnRabh02AvW7ZM\niYmJ+s1vfmN0KADgV+6AvKDaAcSfplzf53tn7dGIaxl9YGYRvQO0ceNGPfroo9qzZ49sNpuWL18u\nl8ulpKQkTZkyRZI0fPhwzZ49O5JhAUDAWh6MZdYgwLzc2TlqkGJ6FjiuZfSBmVm8Xq/X6CCCFcgt\nQ24twh9yBF0hP9AV8gNdIT/gDzkSGVE7BA4AAAAAIoUCCAAAAIBpUAABAAAAMA0KIAAAAACmQQEE\nAAAAwDQogAAAAACYBgUQAAAAANOgAAIAAABgGhRAAAAAAEyDAggAAACAaVAAAQAAADANCiAAAAAA\npkEBBAAAAMA0KIAAAAAAmAYFEAAAAADTsBkdAAAgckpcxcrf/po2Vn4hSRqVMVoTT7pC2fYcgyPr\nKKGkWIlFhbK6KuWxZ6gpN0/u7OiLM16VuIpVVFooV12l7L0zlDsgLyrzxJeW3CnZV6x/96lSxaB+\nSh+SHVP7EGqcT8ARFEAAYBIlrmIt3rhIW/Ztbm0r2luoffX79ItRv4yqD4YJJcVKKshvfW11Viip\nIF8NEh/aIqDEVayCHUf631lb0fo6mvLEl5bc2WR1Kt+2WaqRLFvL5JJHBbUVkqJ/H0KN8wlojyFw\nAGASRaWF2l29q0P7nupdKiotNCCiziUW+Y6ns3aEVmf5EG154ktLjhRa2+e6dffh17GwD6HG+QS0\nRwEEACbhqqtUbXNth/ba5lq56ioNiKhzVpfveDprR2h1lg/Rlie+tORIpaV9rlvqDr+OhX0INc4n\noD0KIAAwCXvvDKXYUjq0p9hSZO+dYUBEnfPYfcfTWTtCq7N8iLY88aUlRzK87XPd2/vw61jYh1Dj\nfALaowACAJPIHZCn49MGd2gflDZYuQPyDIioc025vuPprB2h1Vk+RFue+NKSI3me9rnuOf7w61jY\nh1DjfALaYxIEADCJbHuOfj7qlzExC5w7O0cNErNWGaQlH2JxFriW3MkuKpT2WVpngbObeBY4zieg\nPYvX6/UaHUSwnM5qv8s4HGkBLQfzIkfQFfIDXSE/0BXyA/6QI5HhcKT5bGcIHAAAAADToAACAAAA\nYBoUQAAAAABMgwIIAAAAgGlQAAEAAAAwDQogAAAAAKZBAQQAAADANCiAAAAAAJgGBRAAAAAA06AA\nAgAAAGAatmAWbmho0Ndff62DBw/K6/V2+P1ZZ50VssAAAAAAINQCLoDeeecdzZo1S9XV1R1+5/V6\nZbFYVFJSEtLgAAAAACCUAi6A5s+fr+zsbF177bXq27evLBZLOOMCAAAAgJALuAByOp1auHChTjjh\nhDCGAwAAAADhE3ABdMopp6isrKzHBdDWrVt166236oYbbtDkyZNVWlqq6dOny+12y+FwaN68eerV\nq1ePtgHAHEpcxSoqLZSrrlL23hnKHZCnbHtOQO/L3/6aNlZ+IUkalTFaE0+6IqD3hkpCSbESiwpl\ndVXKY89QU26e3NmR234s6e5xjgaBxB7ufAxH//lbJ/mNrpAfMFrC7NmzZwey4GmnndZanCQkJKi+\nvl7V1dXt/ktLS+tyHbW1tbrrrrt06qmnKiMjQ6NHj9ZDDz2kSy65RHfffbdKSkr0zTff6NRTT/Wz\nnka/8aamJgW0HMyLHIltJa5iFezIV21TjbzyqrapRlurNqtfsl2OFEeX71u8cZE+r1inRnejGt2N\n2l29S99Uf6PBaUNa3xvO/EgoKVZSQb4stTWS1ytLbY1sWzfL088ur6Pz2M2ou8c53ALJj0BiDzQf\nuysc/edvneQ3f1+6Qn4cRo5ERmpqks/2gKfBLi8v165du3TXXXfpkksu0bhx4zr850+vXr20aNEi\n9e/fv7VtzZo1re8dO3asPv7440BDAmBiRaWFQbW3/f3u6l0d2vdU7/L73lBJLPK9nc7azay7xzka\nBBJ7uPMxHP3nb53kN7pCfiAaBDwEbs6cOUpLS9P111/f7UkQbDabbLb2m6yrq2sd8ma32+V0Ov2u\nJz09RTZbgt/lHI6u70gB5EjsqrdW+/yXnXpLdZfHtd5arWZLo3r1an8talaj6q3t3xu2/Kivlnz9\nq1R9tdLIyXa6e5wjwd/2A4k9mHzsjnD0n991kt+S+PvSKfKjFTlinIALoL179+qVV17RiBEjwhaM\nr+8W8qWqqtbvMg5HmpzOjlN2Ay3IkdiW7EmTs7aiQ7sjpX+XxzXZkyabt5dqmmratacmph5e57fv\nDWd+JCenyersGLvH0V/15GQ73T3O4RZIfgQSe6D52F3h6D9/6yS/+fvSFfLjMHIkMjorMgMeAnfi\niSeqqakpZAG1SElJUX19vaTDw+zaDo8DgM7kDsgLqr3t749PG9yhfVDaYL/vDZWmXN/b6azdzLp7\nnKNBILGHOx/D0X/+1kl+oyvkB6JBwJMgDBs2TE8++aSys7Nlt9t7tNFPPvlEvXv31ujRo7V9+3bV\n1dUpKytLixcv1hlnnKFTTjmly/czCQJCgRyJbY4Uh/ol21XVsE/1zXXKSHFo3JDxfme3cqQ4dHza\nEB1qOqT9DVXqldBLZ2SeqZ9lTW733nDmh9fhkKefXdaqfbLU18mT4VDjuPHMguRDd49zuAWSH4HE\nHmg+dlc4+s/fOslv/r50hfw4jByJjM4mQbB4Axx3NmHCBB08eFAHDhxQYmKikpOT26/IYtGaNWu6\nXMfGjRv16KOPas+ePbLZbMrMzNT8+fN19913q6GhQQMHDtTDDz+sxMTELtcTyC1Dbi3CH3IEXSE/\n0BXyA10hP+APORIZnQ2BC/gZoDPOOKNbEx+0NWrUKC1ZsqRD++LFi3u0XgAAAAAIRMAF0COPPBLO\nOAAAAAAg7AIugA4dOuR3mWOOOaZHwQAAAABAOAVcAJ155pl+h8CVlJT0OCAAAAAACJeAC6Bf//rX\nHQqgmpoaffbZZzp48KCuv/76kAcHAAAAAKEUcAF02223dfq7uXPnyul0hiQgAAAAAAiXgL8ItSuT\nJk3SK6+8EopVAQAAAEDYhKQA2r9/f0CTJAAAAACAkQIeAvfXv/61Q5vX69W+fftUUFCgk046KaSB\nAQAAAECoBVwAPfTQQ53+7sQTT9Ts2bNDEQ8AAAAAhE3ABdDKlSs7tFksFvXp04fv/wEAAABMpKam\nRhs2fK61az/Tf/7ziT777FP97GeTNWPGvX6/OsdoARdAgwYNCmccAAAAAKJAZWWlXnppiRYt+ovK\ny8sCft9jj83Vr341VX37HhvG6Hou4AJIklasWKF169Zp//798nq97X5nsVi6HCYHAAAAwDher1c7\nd27XSy8t1cKFT6mhoSEk6x02bLi+852zdMMNN0Z98SMFUQA99thjWrhwoSwWi1JTUzvc2qIAAgAA\nAIzR2Nio9evX6eWXX9KSJYtDvv6TTx6hm266VVdccVXMP/4ScAH0xhtv6MYbb9Ttt9+uXr16hTMm\nAAAAADp816aqap/WrCnSa6/9r5Yt+2fIt3HWWd/TTTf9ShMm/Jd69+4d8vVHm4ALoIMHD+pnP/sZ\nxQ+AkCpxFauotFCuuko5DjTp+19Lo6oS5bFnqCk3T+7sHKNDxLcSSoqVWFQoq6sy6o9Pd2Jtm4v2\n3hnKHZCnURUKaj2+1pFt999HwcYb6HZ8LSfJ53t9xbCxv+9lQ6G7fQUYJRzXQLfbrW+++VqFhR/p\n9ddf0erVH4Qo2iMmTPihfvGLXyo39xylpKSEfP2xyOI9+mGeTvzsZz/TtGnT9L3vfS/cMfnldFb7\nXcbhSAtoOZgXOWK8ElexCnbkS5KslU4lbN0sSZrYnKVTPA5JUsOlEw35kE1+tJdQUqykgvwO7UYd\nn650J9a2udjCWunUFcVqzcW26+l33vc65IevdUjSpcMndvnBPth4A92Or+WcdU5ZJGX0br9PEy2n\n6vT3NrRr22R16rUcyZPRfll/+xOI7vZVrOD6EX+6c105dOiQtm3botWrP1B+/uvauPGLkMd15ZU/\n1TXXTNHpp4/RMcekhXz9sc7h8N0nAd8Bmjlzph599FHNnDlTo0aNCllgAMyrqLSw9Wfr7l2tPxda\nd7V+6EwsKoy6D9hmlFhU2Gl7tB2f7sTaNhdbWHfvUqG1YwGUWFQondfxHwN9raOlvasP9cHGG+h2\nfC23p/rweXZ0AfTJ+td1uk5u11Zo3SXr7o4FkL/9CUR3+wowStvz1OP1yllbow0VFXrn3un6+6YN\nqqqqCvk2b7zxJl166USdcsoo9enTN+qnlo4lARdAM2bM0MGDBzVp0iQlJiYqOTm53e8tFovWrFkT\n8gABxC9XXWXrz5a62tafKy1Hfra6KgXjdXYcovH4dCfWtrnYwlJXq0ofnzc6W4+vdXTV7m99Pd2O\nr+Vqm2s7tEmSq6ZCOqoAqrTUylIX+PaD0d2+AsKtvr5ee/bs0uefr9dbb70ZludtWgwePERXX32t\nxo4dp1GjRnf4bI3wCbgAOu2006g8AYSUvXeGnLUVkiRv7xRZamskSRneI2OUPfYMQ2JDex57hqzO\nCp/t0aY7sbbNxRbe3inKOOR7/YGuo6U9lPEGuh1fy6XYfI//t6f2lxrbt2V4U1Tu41lof/sTiO72\nFdBdXq9X+/bt0zfffKVPPinS0qUvasuWzWHb3jnnfF+TJl2tceMmqH///h0+QzNM0lgBF0CPPPJI\nwCt98803dcEFF/CgFYAu5Q7Ia30OwHP84NZngPI8g1uXacrNMyQ2tNeUm+dz/Hs0Hp/uxNo2F1t4\njh+svGLf6w90HS3toYw30O34Wm5Q2mD5+qfM755+uXTUM0B5nsF67Xjf2++p7vYV4EtjY6PKykq1\nfftWvfnmMi1d+mLYt3nJSSfrV2ecqdE33qykMWeEfXsIraC+CDVQv//973XaaadRAAHoUstY/6LS\nQrkcVvXr1U/f/1o6pSpRHkd0zzJmNu7sHDUouBnRjNKdWNvlYsusZN+bqJNOlDwBrsfnOgKY2SzY\neAPdjq/lLh0+0ed7R9hz1HDM8HYxnJQ7UT8K0yxw3e0rmIvX69WBA/u1d+9effDBe3riif9RZaUz\nItu+9trrNHXq7Ro+/MjQ0FiaCRNdC3gWuGCMGTNGy5Yt0+DBg/0v3A3MAodQIEfQFfIDXSE/0BXy\nw7/m5maVl5dp3bq1evbZp1VY+FHEtn3yySP07LN/VVZWtmGPd5AjkdHjWeAAAAAAfw4dqtb69ev0\n6qsv66WXlkR029/5zpl67LEnlc2dGXSBAggAAAB+ud1ubdjwud54459auPApNTU1RXT73/1urpYu\nfVnHHpse0e0i/lAAAQAAmJjH49GXX+5QQcEbevHF57Vnz+6Ix/CTn1yumTPv04knDmPWYYQdBRAA\nAEAcam5uVmnpXq1YsVwvvvi8Sko2GRbLkiUva8KEH1LcICpQAAEAAMSQlmmfi4oKtWTJC1qz5mPD\nYrngggv11FOL1K+f3bAYgGBRAAEAAESBuro6lZWV6osv1uuVV/6hd95529B4Fiz4sy6/fJKSk5MN\njQMINQogAACAMDp0qFplZWXatGmD/vWvN/X6668aHZLuvnuWrrlmio47boDRoQAR1+MCaN++ferX\nr1+7tltuuUV9+/bt6aoBAACiktfr1f79VSovL9fWrZv1zjtv67XX/ldut9vo0PT008/pkkt+ol69\nehkdChCVAi6AqqurNWPGDP3ud7/TSSedpJ07d+qmm27Snj17NHLkSD3zzDPKzMyUJN18881hCxgA\nACBcPB6PXC6XysvLtGPHNq1YsVxvvfWmDh0y/ksrhw0brqeffk6nnTaGyQSAHgi4AJo3b56++uor\nHXPMMZKkhx56SM3NzZo+fbqWLVumJ554Qv/93/8dtkABAAC6q7m5WZWVTpWXl2nnzh1auXKF3nnn\nX9q/f7/RoUmSLr/8St1772wNHjzE6FCAuBdwAfTBBx9o3rx5Ou6443TgwAF9/PHHmjdvni6++GKN\nHj1ad911VzjjBAAA6KCxsVEVFeUqKyvVN998rVWr3tWKFW+rqqrK6NAkSaNHn6758xdw1waIIgEX\nQPv27dPQoUMlSZ988oksFot+8IMfSJKOP/54OZ3O8EQIAABMp7a2VuXlZSovL9fu3d/o/fdXRdUd\nmxbLl7+n008/o11x43CkyenehvtwAAAgAElEQVQ0fsgcAN8CLoDS09NVVlamzMxMrVq1SqNHj1Zq\naqokyel0tv4MAADgi9fr1aFD1SovL1d5eZn27NmtDz98PyoLm/T0dK1Y8aGGDBlqdCgAQizgAui7\n3/2uZs2apTPPPFPLli3TfffdJ+nwv9AsXLhQp556atiCBAAA0avtjGjl5WXau3eP/v3v1VE1FK2t\np59+TpdddiVD0gCTCrgAuuuuu3THHXfon//8py6++GJNmjRJkvTuu+9q9erVWrJkSdiCBAAAkdd2\nRrSKijLt2bNHn3xSpHfe+VdUFjYZGQ6tWbNOaWl9jA4FQBQLuADKzMzU0qVLO7Sfd955WrVqldLT\n07sVQE1NjWbMmKEDBw6oqalJv/71r/X973+/W+sCEHtKXMUqKi2Uq65S9t4Zyh2Qp2x7jtFhhUxC\nSbESiwpldVXKY89QU26e3Nnxs3+IvLY5tTG9SauHSs6+iUGdP21nRCsrO3zHZu3a/0RtYSNJ7777\noUaPPt3QGLq6XrU9LjrheCWcckZQ53qkroXRcE0yKoZwbjfe/5b5Ew15FQyL1+v1GhnA0qVLVV5e\nrjvuuEPl5eW6/vrr9fbbb3f5nkAeLOQBRPhDjhivxFWsgh35HdovHT7R8D8cociPhJJiJRV03L+G\nSydG9R8G+GfU9aNtTm2yOpVv2yxJco/IkifDoeamZp2ddo7SGtNUXl6usrK9WrdubdQORZOkOXMe\n1k033Rr1w9G6ul6NqlC7cz01NUk1NQ0Bn+uRuhZGwzXJqBjCud3uHL94+gwSDXnVGYcjzWd7l3eA\nsrKygroglZSUBBeVDj9kuGXLFknSwYMHu30nCUDsKSot7LTd6AIoFBKLfO9fYlGh4X8UEFtaZkTb\n//elKt+zW6WHqvW/FcXatNOp+rpmSf9uXXaB5hsXaBs333yrZsyY1fr9gbGuq+vVmE99vyfQcz1S\n18JouCYZFUM4txvvf8v8iYa8ClaXBdCUKVPaFUD/+te/lJqaqjFjxig1NVUHDx7U2rVr1dzcrCuv\nvLJbAfzoRz/S66+/rvHjx+vgwYN65pln/L4nPT1FNluC3+U6q/qAFuSIseqt1UpNTerYbqmOimPT\n4xjqqyUf+6f6aqVFwf6hZ3qaH16vV9XV1SotLW39b/fu3Vq/fr3efvtt7du3L0SRhlZ6erq2bdsm\nu91udCgR1dX1Kq1eHc711NSkgM/1iF0Lo+GaZFQMYdxud49fNPydC4loyKsgdVkA3Xvvva0/P/PM\nMzrvvPP04IMPtlvG6/Xq7rvvls0W8ONE7bzxxhsaOHCgnnvuOW3evFn33HOPXn/99S7fU1VV63e9\n8XRrEeFBjhgv2ZMmZ21Fh3ZHSn/Dj00o8iM5OU1WZ8f98zj6q57ci2ld5cfRM6KVlZWqrKxUxcUb\n9f77q6J2KJokvfzyPzV27Di/y3k8gQ1HjyddXa+qk9XuXG8ZAhfouR6pa2E0XJOMiiGc2+3O8Yun\nzyDRkFed6dYQuLb+/ve/a9GiRR3aLRaLfvGLX+imm27SLbfcEnRga9eu1bnnnivp8JC7iooKud1u\nJST4v8MDILblDsjzOW46d0CeAdGEXlNuns9x0U258bF/ZtN2RrT6+gPauvXLbwubTfroow+i9o6N\nJP3xj3/SlCk3GB1GTOvqetWUqx6d65G6FkbDNcmoGMK53Xj/W+ZPNORVsAIugPbt26fGxkafv2tu\nbu72hX/o0KH6/PPPddFFF2nPnj1KTU2l+AFMomVsdLzOnOPOzlGDFFMz45hRy4xoZWWl7e7alJQU\nq7BwdVTfsfn1D8/TCdd9L+hZ4BC8rq5XbrvanevKzFRDELPARepaGA3XJKNiCOd24/1vmT/RkFfB\nCngWuCuuuEJWq1Vz5sxRVlZWa3txcbHmzJmjhoYGv0PXfKmpqdE999wjl8ul5uZm3X777Tr77LO7\nfA+zwCEUyBF0hfyIfQ0NDaqoOFzQtBQ25eWHC5s1az6O6sLm8cef0s9+NtnoMNBNXD/gDzkSGT0e\nAnfffffp5ptv1mWXXaaEhAQlJyeroaFBzc3N6t27d0CTF/iSmpqqxx9/vFvvBQCYT8uMaOXl5aqo\nKFN5eZlKS0u1ZUuJ/vOfT6K6sAn0GRsAQPgEXACdfvrpWrFihZYvX67t27fr0KFD6t27t0466SRN\nmDBB/fr1C2ecAIA45vV6dehQdbshaIe/x+ZwYfPpp5/o0KHo/dfS//xng4YMGWp0GACAAAQ1dVuf\nPn00adKkcMUCAIgzbWdEO1zUHLlzs3nzZn344XtGh9ip008fozffXKFevXp1+B3DVwAgdgVVADmd\nTv3tb3/Txo0bVVFRoWeeeUYZGRlavny5LrnkknDFCACIMi0zopWVlX47DO3InZvi4k1as+Zjo0Ps\n1JQpP9fcuY8x4Q4AmFTABdC2bds0efJkNTQ0KCsrSzt37lRzc7N27dqlGTNmyGq16uKLLw5nrACA\nMGtubpbTWXHUxAFlKisr02effari4o1Gh9ipNWvW68QThxkdBgAgygVcAM2dO1fZ2dl67LHH1K9f\nP40ZM0aSNGzYME2bNk3PP/88BRAARKnOZkTbs2ePCgs/0u7du4wO0adjjz1Wn376hfr2PdboUAAA\ncSLgAmjdunV64YUXfE52MGHCBD355JMhDQwA4N/RM6K1TB6wffs2/fvfq3Xw4AGjQ/TpmWee18SJ\nV8hisRgdCgDAZAIugCwWi1JTU33+rrGxkT9iABAibWdEaztxwN69u7Vx4wYVFn5kdIid2rr1ax17\nbLrRYQAA0KmAC6CTTjpJixcv1gMPPNDhd8uWLdPIkSNDGhgAxJuWGdHKyspan60pLy/Xjh3b9Pnn\n66P2+Zo777xbd901k3/oAgDEhYALoOuvv17Tpk1TcXGxzjnnHLndbv3tb3/Tzp079dFHH+lPf/pT\nOOMEgKjl8XhUWVmp8vKy1hnRdu/epU2bNmrDhs+j8vmaH/7wYj3xxNM8WwMAMB2L1+v1Brrw8uXL\n9eSTT2rbtm2H32yxaMSIEbr11lt10UUXhS3IowXy3Qt8RwP8IUfQFYcjTaWlVe1mRCsrK9X27du0\nYcPn2rDhC1VXHzQ6zHYyM4/TG2/8S8OGDTc6lLjH9QNdIT/gDzkSGQ5Hms/2oL4H6KKLLtJFF12k\nmpoaHTp0SGlpaUpJSQlJgAAQKUfPiLZnzy5t3LhBGzZ8oU2bNhgdXgf/+MfrGjt2HEPQAAAIgaAK\nIOnwjEPbt29XeXm5zj33XEmHh39YrdaQBwcAwTgyI9rh/7Zv36Yvvvhcn332qSoqyo0Or53nnvur\nLrnkJxQ1AABEWMAFkNfr1YIFC/TCCy+ooaFBFotF77zzjqqrq3XjjTfqhRdeUEZGRjhjBWBCLTOi\ntUwcsGfPkZnQNmz43Ojw2rn++hv14IOPqlevXkaHAgAAOhFwAbRw4UK98MILuvbaa5WXl6epU6dK\nkpKTk5WSkqLHH39cc+bMCVugAOKL1+tVVdW+1i/l3LKlRB988J5WrlyhIB5NDLt58xbommumKDEx\n0ehQAABACARcAL322mu69957ddVVV0lS67CNvn376s4779Tvfve78EQIIKa0nRFtx45tev/9VVq5\ncoXKy8uMDq3VI4/8UVOm3NBlUWOmB1RLXMUqKi2Uq65S9t4Zyh2Qp2x7jtFhhU1CSbESiwpldVVK\nTU2HGxMT5bFnqCk3T+5sY/a9Ja6EzcWyVlXJk95P7qzsoGNqu39G71NPRGo/4qW/AAQu4AKorKxM\nZ599ts/fDRo0SPv37w9ZUACiT3NzsyoqyrV5c4mKigq1atW7+uKL9UaH1eo3v/md7rhjhnr37m10\nKDGlxFWsgh35ra+dtRWtr+OxCEooKVZSweH9szidsm3dLElyj8ySPB4lFeSrQYr4B+CWuNrGZC0v\nk8XrkdVZEXBMbfdPkqzOCsP2qScitR/x0l8AghNwAZSZmant27dr8ODBHX63fft2pafzzd9ALDp0\n6JCKizfp88/X6r33VmrVqnfl8XiMDkuS9Je/PKvLLruSSVbCqKi0sNP2eCyAEouO7G/CniPfz2Td\nvUueDEfrMpH+8NsSV9uY2sYVaExt9+/o9lj6QB+p/YiX/gIQnIALoO9+97uaPXu2LBaL8vLyWtu/\n+OILPfLII7rgggvCEiCA4Hi9XtXUHNLXX3+tjRu/aL1bU1q619C4xo0br9tu+63OOONMJScnGxoL\njnDVVQbVHuusriP7Zamt9flz22UipWWbbeNo+zrQmDpbzoh96olI7Ue89BeA4ARcAN15553atGmT\nbrnlFlksFnm9Xv3Xf/2XmpublZOTwzNAQBi53W7t27dPTmeFvvxypz7/fJ3ef3+l1q9fZ2hcJ544\nTA888LDGjh3HzGcxyt47Q87aCp/t8chjz5DVeXh/vSkpstTUtP7cdhmj4mobU9u4Ao2p7f4d3R5L\nIrUf8dJfAIITcAGUnp6u1157TStXrtS6detUXV2tPn366IwzztDYsWOVkJAQzjiBuFNfX6/KSqec\nzgo5nRXaunWrior+rVWr3lVzc7NhcZ1yyqn6xS9+qUmTruZOjQnkDshr9wxQ2/Z41JSb1/rMh3vQ\n4NbnbTzHD263jFFxtY2pbVyBxtR2/45ujyWR2o946S8AwQm4AGpsbNS8efN03XXXacKECXI6nfrd\n736nl156SWeffbbmzZunY445JpyxAlHN6/Xq4MEDcjqdrYVNRUW5Nm7coJUrV6isrDTiMfXte6yu\nuGKSrrjiKo0Z8x3ZbEF/9zHiXMtzPmaZBc6dnaMGHX7Gw2q1qqlfv8O/MHgWuLZxyWrp9ixw7fYv\nhmc1i9R+xEt/AQiOxRvgF27MnTtXb775ppYsWaKhQ4dq2rRpKiws1I9//GOtWrVK48eP18yZM8Md\nryQFNDWtmaawRfcEkiNNTU3at8+lioqK1qKmsrJSe/fu1ttvv6Vvvvk6QtEe4XD016RJV+vyy6/U\nqFGjmSAgTLiGoCvkB7pCfsAfciQyHI40n+0B/3Pw8uXL9cADD2jo0KGqra3VypUr9fvf/16TJk3S\n2LFjNXv27IgVQEBP1NTUyOms0I4dtdq69StVVjpVXl6mbdu2aO3az/T111+FPYYhQ07Q2LHjdM45\n52rEiCw5HP3Vr18/hpICAACEWcAFkNPpVHZ2tiTp008/ldvt1vjx4yVJw4cPV0VFx4cIgUjweDyq\nqqpqc4empaDZpi1bSrR5c4kOHAj/91TZ7Xadf/44jR07Trm5eRo4cBBDzgAAAKJMwJ/O0tLSVFVV\npczMTK1evVrZ2dk69thjJUkHDhxQUlJS2IKE+TQ0NMjlqmwtaCoqKrRt21Zt3lyszZtLtHv3Lv8r\nCYHk5GSdf/44nX/+BTr//LE64YRhDDkDAACIYQEXQKeddpoeeughnXfeeXr11Vc1depUSYcf/P77\n3/+ukSNHhi1IxD6v16tDh6q/nfGsUhUV5dqxY5s2by5WSUmJNm8uVoCPowWsX79+ysrKUVZWtoYN\nG66MDIcyMhyy2zOUkeHQyJFDtX9/verr65ntDAAAwCQCLoDuuOMO/fKXv9T8+fN15plnavLkyZKk\nZcuW6ZVXXtHChQvDFiSik9vtlsvlap3GeceObSopKQnbsLPExESNHJmtrKxsZWXl6Pjjj5fdniG7\nPUMOh0P9+tmVmJgY1Pokih8AAAAzCbgAGj58uFatWqV9+/apX8u0oZJyc3OVn5+vk08+OSwBIrLq\n6uraFDTbtXnzkYJm165vQr69oUNPUHZ2jkaOzNbJJ49Q//6Z396pyVC/fna+XBMAAAAhFfQT2m2L\nH0nKzMxUZmZmyAJCaHm9Xu3fXyWn06mdO3d8+wxNceudGo/HE9LttR12lpWVo6FDT2hX0PCsGAAA\nAIzEFFUxqKmpSZWVTn355U6VlByeFKBlcoBQDzuz2WztCpqRI7M0cODA1qFnFDQAAACIJRRAUcDr\n9aqm5pC+/PLL1kKm5f/hGHY2ZMgJys7Obi1sTjrp5NbJAXgeBgAAAPGMAihMPB6Pdu36prWYaXun\nJlzDzkaOzFJWVo6ys3M0ePAQChoAAADgKBRAPbB9+zb9+Mc/VGWlM6TrPXrYWVZWlkaMyFJm5nHq\n3bt3SLcFAAAAmAkFUA889NADARc/Rw87y8rK0QknnKjU1NQwRwkAAACgBQVQDzz88HwNGzZcgwYd\n/+1UzllKT+/n/40AAAAADEEB1AOZmZmaNWu20WEAAAAACJDV6AAkadmyZfrxj3+syy+/XO+//77R\n4QAAAACIU4bfAaqqqtKf//xnvfbaa6qtrdUTTzyh888/3+iwAMOVuIpVVFooV12l7L0zlDsgT9n2\nnLBtL6GkWIlFhbK6KuWxZ6gpN0/u7PBtLxZFso96evy7en/b/diY3qTVQyVn30Sf24l0HkaTaDkn\njj4G59RlavTG8m7H1d1jGmhOtcQkKSr6LxSi9TyIlhwNt2jt/1CJxeMYizG3lTB79uzZRgbw3nvv\nyWq1asKECUpNTdUFF1zg9z21tY1+l0lNTQpoOZhXNOdIiatYBTvyVdtUI6+8qm2q0daqzeqXbJcj\nxRHy7SWUFCupIF+W2hrJ65Wltka2rZvl6WeX1xH67cWCo/Mjkn3U0+Pf1fuP+9rZuh+bLBV6o3md\n6ip2ydu7t2oS1W47kc7DaOLveEfq+nH0Magv/Urb1/9LGTVe9fekBJ2H3T2mgeZUS1/1Klwt27q1\nskgxf03pTp9FIj/Mct2O1+tQS47E4nGMpZhTU5N8ths+BG737t2qr6/XLbfcomuuuUYff/yx0SEB\nhisqLQyqvacSi3yvt7N2M4pkH/X0+Hf1/rbxFlp3tf5s3b2r3XKhiCOWRcs5cXRftxyntsdOCjyu\n7h7TQHOqhXX3LiXs2dWhPRavKdF6HkRLjoZbtPZ/qMTicYzFmI9m+BA4Sdq/f7+efPJJ7d27V9dd\nd53ee+89WSyWTpdPT0+RzZbgd70OR1oow0QcitYcqbdW+/xXi3pLdXhirq+WfP0rSX210qK0jyKh\nXV9HsI96evy7en9avVr346Aa1avlz0Bzo5K+bW/ZTsTzMJoEcLwj0QcdjkFzo5Ro00E1KjWpTXuA\nedjdYxpoTrWLU2rNqWDjjCbd7bOw54dJrtvxfB1yONJi8zjGYsxHMbwAstvtGjNmjGw2m4YMGaLU\n1FTt27dPdru90/dUVdX6Xa/DkSanszqUoSLORHOOJHvS5Kyt6NDuSOkflpiTk9NkdXbcnsfRX/VR\n2kfhdnR+RLKPenr8u3p/dbJa96OPrZcqrDWSJG9KqpprGtptJ9J5GE38He9IXT+OPgY2Wy9ZamvU\n35OqmuaGDnEFu74W/o5poDnVNk5JrTkVbJzRpDt9Fon8MMt1O16vQy05EovHMZZi7qxINnwI3Lnn\nnquioiJ5PB5VVVWptrZW6enpRocFGCp3QF5Q7T3V8sByoO1mFMk+6unx7+r9bePN8wxu/dlz/OB2\ny4UijlgWLefE0X3dcpzaHjsp8Li6e0wDzakWnuMHyz1ocIf2WLymROt5EC05Gm7R2v+hEovHMRZj\nPprhkyAcc8wxqq+v1/3336+CggLdcccdGj58eJfvYRIEhEI054gjxaF+yXZVNexTfXOdMlIcGjdk\nfNhmvfE6HPL0s8tatU+W+jp5MhxqHDc+pmZ0CbWj8yOSfdTT49/V+9vuR2adVf36DpJzxGDV9uvT\nYTuRzsNo4u94R+r6cfQxsGecoHEDxmrUoZRu5WF3j2mgOdUSU8OlE9U85jtxcU3pTp9FIj/Mct2O\n1+tQS47E4nGMpZg7mwTB4vV6vRGOpccCueUZzcObEB3IEXSF/EBXyA90hfyAP+RIZETtEDgAAAAA\niBQKIAAAAACmQQEEAAAAwDQogAAAAACYBgUQAAAAANOgAAIAAABgGhRAAAAAAEyDAggAAACAaVAA\nAQAAADANCiAAAAAApkEBBAAAAMA0KIAAAAAAmAYFEAAAAADToAACAAAAYBo2owOIVSWuYhWVFspV\nVyl77wzlDshTtj3H6LCAoHWWy7GQ4wklxUosKpTVVSmPPUNNuXlyZ8dG7LGqsz6PZuGMORL9EYt9\nDgDRjAKoG0pcxSrYkd/62llb0fqaD1mIJZ3l8s79O7TJtaFDuxQ9OZ5QUqykgiOxW50VSirI1/pD\nO1Tgje7YY1Vnfd4gRe0H8nDGHIn+iMU+B4BoxxC4bigqLQyqHYhWneXsGzteD2p5IyQW+Y7lk/XR\nH3us6qzPO2uPBuGMORL9EYt9DgDRjgKoG1x1lUG1A9Gqs5x11lYEtbwRrK5OzsOa6I89VnXW5521\nR4NwxhyJ/ojFPgeAaEcB1A323hlBtQPRqrOcdaT0D2p5I3jsnZyHqdEfe6zqrM87a48G4Yw5Ev0R\ni30OANGOAqgbcgfkBdUORKvOcvYnwy8PankjNOX6juW7p0d/7LGqsz7vrD0ahDPmSPRHLPY5AES7\nhNmzZ882Oohg1dY2+l0mNTUpoOW6w5HiUL9ku6oa9qm+uU4ZKQ6NGzKeB6xjTDhzJFZ0lst5g86N\n+hz3Ohzy9LPLWrVPlvo6eTIcahw3Xv1OC03s5EdHnfV5ND+MH66YU1OTdCilb9j7Ixb7HFw/4B85\nEhmpqUk+2y1er9cb4Vh6zOms9ruMw5EW0HIwL3IEXSE/0BXyA10hP+APORIZDkeaz3aGwAEAAAAw\nDQogAAAAAKZBAQQAAADANCiAAAAAAJgGBRAAAAAA06AAAgAAAGAaFEAAAAAATIMCCAAAAIBpUAAB\nAAAAMA0KIAAAAACmQQEEAAAAwDQogAAAAACYBgUQAAAAANOgAAIAAABgGhRAAAAAAEzDZnQALerr\n63XJJZfo1ltv1eWXX250OEBUKHEVq6i0UK66Stl7Zyh3QJ6y7TlGhwX4lVBSrMSiQlldlfLYM9SU\nmyd3dvvcDVd+c96YRyB5BgBHi5oC6C9/+Yv69u1rdBhA1ChxFatgR37ra2dtRetrPswhmiWUFCup\n4EjuWp0VSirIV4PU+uE0XPnNeWMegeQZAPgSFUPgduzYoe3bt+v88883OhQgahSVFgbVDkSLxCLf\nOdq2PVz5zXljHoHkGQD4EhV3gB599FHdd999ys/P97+wpPT0FNlsCX6XczjSehoa4lw050i9tVqp\nqUkd2y3VUR13PKGfu6m+WvKRu6qvVtq3fRqu/I7keUN+GCyAPDMS+QF/yBHjGF4A5efn6/TTT9fg\nwYMDfk9VVa3fZRyONDmd1T0JDXEu2nMk2ZMmZ21Fh3ZHSv+ojjteRHt+RLPk5DRZnR1z1+Por/pv\n+zRc+R2p84b8MF4geWYU8gP+kCOR0VmRafgQuPfff18rV67UVVddpVdeeUVPPfWUCgu5fQ3kDsgL\nqh2IFk25vnO0bXu48pvzxjwCyTMA8MXwO0ALFixo/fmJJ57QoEGDlJfHxQtoeWCb2awQa9zZOWqQ\nupydK1z5zXljHoHkGQD4YngBBKBz2fYcPrghJrmzc/x+EA1XfnPemEcgeQYAR4uqAui2224zOgQA\nAAAAcczwZ4AAAAAAIFIogAAAAACYBgUQAAAAANOgAAIAAABgGhRAAAAAAEyDAggAAACAaVAAAQAA\nADANCiAAAAAApkEBBAAAAMA0KIAAAAAAmAYFEAAAAADToAACAAAAYBoUQAAAAABMgwIIAAAAgGnY\njA4AMMrGio16a9MKueoqZe+dodwBecq25xgdFhBRJa5iFZUWch4APiSUFCuxqFBWV6U89gw15ebJ\nnc35AcQ6CiCYUomrWKvK/qWa2gZJkrO2QgU78iWJD38wjRJXcWveS5wHQFsJJcVKKjhyflidFUoq\nyFeDRBEExDiGwMGUikoLg2oH4hHnAdC5xCLf50Fn7QBiBwUQTMlVVxlUOxCPOA+Azlldvs+DztoB\nxA4KIJiSvXdGUO1APOI8ADrnsfs+DzprBxA7KIBgSrkD8oJqB+IR5wHQuaZc3+dBZ+0AYgeTIMCU\nsu056tcvlVngYGot+c4scEBH7uwcNUjMAgfEIQogmNao/qOUaRlqdBiAobLtORQ8QCfc2TkUPEAc\nYggcAAAAANOgAAIAAABgGhRAAAAAAEyDAggAAACAaVAAAQAAADANCiAAAAAApkEBBAAAAMA0KIAA\nAAAAmAYFEAAAAADToAACAAAAYBoUQAAAAABMgwIIAAAAgGlQAAEAAAAwDQogAAAAAKZBAQQAAADA\nNGxGB9Bi7ty5+uyzz9Tc3Kybb75ZEyZMMDokAIhbJa5iFZUWylVXKXvvDOUOyFO2PcfosAAACLuo\nKICKioq0bds2vfzyy6qqqtJll11GAQQAYVLiKlbBjvzW187aitbXFEEAgHgXFQXQWWedpdGjR0uS\n+vTpo7q6OrndbiUkJBgcGQDEn6LSwk7bKYAAAPEuKgqghIQEpaSkSJJeffVVnXfeeV0WP+npKbLZ\n/BdHDkdayGJEfCJH0JV4zY96a7VSU5M6tluq43afw4G+QlfID/hDjhgnKgqgFu+++65effVVPf/8\n810uV1VV63ddDkeanM7qUIWGOESOoCvxnB/JnjQ5ays6tDtS+sftPodaPOcHeo78gD/kSGR0VmRG\nzSxwq1ev1tNPP61FixYpLY2KGADCJXdAXlDtAADEk6i4A1RdXa25c+fqhRde0LHHHmt0OAAQ11qe\n82EWOACAGUVFAfTWW2+pqqpK06ZNa2179NFHNXDgQAOjAoD4lW3PoeABAJhSVBRAP/3pT/XTn/7U\n6DAAAAAAxLmoeQYIAAAAAMKNAggAAACAaVAAAQAAADANCiAAAAAApkEBBAAAAMA0KIAAAAAAmAYF\nEAAAAADToAACAAAAYBoUQAAAAABMw+L1er1GBwEAAAAAkcAdIAAAAACmQQEEAAAAwDQogAAAAACY\nBgUQAAAAANOgAAIAAGgyegIAAAtsSURBVABgGhRAAAAAAEwj7gqgV155RVOmTGn9b8yYMUaHhChS\nU1OjqVOnasqUKbr66qu1evVqo0NClPF4PLrvvvt09dVXa8qUKdqxY4fRISEKbN26VRdeeKGWLl0q\nSSotLdWUKVN0zTXX6Pbbb1djY6PBEcJoR+eIJP31r3/VKaecopqaGgMjQzTwdQ254YYbNHnyZN1w\nww1yOp0GR2gucVcATZo0SUuWLNGSJUt02223aeLEiUaHhCjyz3/+UyeeeKKWLFmixx9/XA8++KDR\nISHKrFy5UtXV1frHP/6hBx98UHPnzjU6JBistrZWc+bM0dlnn93a9qc//UnXXHONXnrpJQ0dOlSv\nvvqqgRHCaL5yJD8/Xy6XS/379zcwMkQDX/mxYMECXXXVVVq6dKnGjx+vxYsXGxih+cRdAdTWn//8\nZ916661Gh4Eokp6erv3790uSDh48qPT0dIMjQrT56quvNHr0aEnSkCFDtHfvXrndboOjgpF69eql\nRYsWtfsgu2bNGo0bN06SNHbsWH388cdGhYco4CtHLrzwQv32t7+VxWIxMDJEA1/58Yc//EEXXXSR\npPafTRAZcVsAffHFFxowYIAcDofRoSCK/OhHP9LevXs1fvx4TZ48WTNmzDA6JESZESNG6KOPPpLb\n7dbOnTu1a9cuVVVVGR0WDGSz2ZScnNyura6uTr169ZIk2e12hq+YnK8cOeaYYwyKBtHGV36kpKQo\nISFBbrdbL730ki699FKDojOnuC2AXn31VV122WVGh4Eo88Ybb2jgwIFasWKFXnzxRT3wwANGh4Qo\n84Mf/ECnnnqqrr32Wr344osaNmyYvF6v0WEhipEfALrD7XZr+vTpys3NbTc8DuFnMzqAcFmzZo1m\nzZpldBiIMmvXrtW5554rScrKylJFRYXcbrcSEhIMjgzR5Le//W3rzxdeeKHsdruB0SAapaSkqL6+\nXsnJySovL+c5DwBBmzlzpoYOHaqpU6caHYrpxOUdoPLycqWmprYOTwBaDB06VJ9//rkkac+ePUpN\nTaX4QTubN2/WzJkzJUkffvihcnJyZLXG5aUSPZCXl6fly5dLkt555x19//vfNzgiALFk2bJlSkxM\n1G9+8xujQzElizcO791v3LhRCxYs0LPPPmt0KIgyNTU1uueee+RyudTc3Kzbb7+d285ox+Px6J57\n7tH27duVlJSk+fPna8CAAUaHBQNt3LhRjz76qPbs2SObzabMzEzNnz9fd999txoaGjRw4EA9/PDD\nSkxMNDpUGMRXjuTl5amwsFDr16/XqaeeqtNPP13Tp083OlQYwFd+uFwuJSUltT4rNnz4cM2ePdvY\nQE0kLgsgAAAAAPCFcR0AAAAATIMCCAAAAIBpUAABAAAAMA0KIAAAAACmQQEEAAAAwDQogAAAAACY\nBgUQACAoF1xwgW699dawrLukpEQjR47U66+/Hpb1d2XKlCn6yU9+EvHtAgAiy2Z0APj/7d1fSFN/\nHwfwdzpD3ZSsLTWzBZUWWXnCVBIpLU0JoUybaXqRiWEqiZhGRGTigkolYfT3oj9CMVeZphGaiHlh\nNSTF9EL8E602J6Uu/yDYnotwuEcffu15yOrp/YJzcb7/zmfnZnw4n/M9RET0OygvLwc/jUdE9P+P\nCRARERGAJUuW/OoQiIhoAbAEjoiI/isajQYRERHw8/PDnj17UFdXZ+nr7e1FZmYmAgMD4efnh8jI\nSKhUKkxPT1vGjI2NoaCgAAEBARAEARkZGRgcHLT019fXw9fXFy0tLXOuHRUVhfT0dJvivXnzJqKj\no7FlyxZs27YNKSkpeP36taV/dglca2srfH195z0KCgoscz58+ICcnByEhYVh06ZN2Lt3Lx48eGBT\nXEREtLD4BIiIiGz27t07TE9Po7S0FA4ODigtLUVubi7WrFkDqVSKpKQkeHl5QaVSQSaTobm5GRcu\nXIDJZEJ+fj4A4Pz586irq0NhYSH8/f2h1WpRXFxsuUZYWBhWrFgBjUaDkJAQq2v39fUhLy/vh+PV\naDS4dOkSioqKEBwcjImJCdy9exdHjhxBQ0MDli9fbjVeEAS8fPnSqk2tVqO8vBzR0dEAgNHRUSQl\nJUEsFqO4uBgeHh54/vw5zp49C7PZjISEBJvvKxER/XxMgIiIyGYjIyMoKSmBWCwGACiVSmzfvh3V\n1dVwcXHBly9foFarsXLlSgCAXC5Hd3c37t+/j5ycHExPT+Pp06dQKBSWpy5yuRwGgwFlZWUAAHt7\nexw8eBAqlQrDw8OWErWamhq4u7tj586dPxxvR0cHnJ2dsW/fPohE3//6zpw5g/3791t+w2yLFy+G\nTCaznGu1WqhUKpw4cQI7duwAAFRWVkKv16Oqqgrr168HAKSnp6OzsxMqlQoKhQKLFi2y5bYSEdEC\nYAkcERHZbMOGDVaJg5ubGzw9PdHT04OOjg54enpakp8ZgiBgfHwcvb29GBgYwNTUFPz8/KzGbN26\n1eo8Pj4eZrMZ1dXVAACz2Yza2lrExcXB3t7+h+PdtWsXJicnoVAoUFFRgZ6eHjg4OEAQhHkToNk+\nffqE7Oxs7N6926rsrq2tDTKZzJL8zAgJCYHBYIDRaPzh+IiIaOHwCRAREdnM1dV1TptYLMbExAQA\nQCKRzOmfaRsbG7Pstubs7DxnjdmkUikiIyNRWVmJ5ORkaLVaDA4OIj4+3qZ4Q0NDUVFRgTt37uDK\nlSsYHh6Gt7c3MjIyEBsb+x/nTU5O4vjx45BKpVAqlVZ9JpMJQ0NDEATBqn3mPSeDwTCntI6IiH49\nJkBERGSzsbGxedvkcjns7OzQ19c3p99kMgH4njxNTU0BgCVh+vcxsyUmJiIpKQldXV2oqalBaGgo\nPD09bY5ZEAQIgoBv376hvb0dN27cwKlTp+Dl5YWgoKB555w+fRo6nQ4ajQZOTk5Wfa6urvDw8MDt\n27fnnevu7m5zjERE9POxBI6IiGzW2dmJr1+/Ws6Hhobw8eNH+Pj4YPPmzdDr9Xj//r3VHK1WCxcX\nF6xevRpyuRwikQjt7e1WY2bvyjYjICAAPj4+ePz4Merq6qBQKGyOt7m5GR0dHQAAOzs7+Pv7o6io\nCADw9u3beedcv34dz549Q1lZ2ZxyPuB7uZ7RaISjoyPkcrnlcHJygouLCxwdHW2Ok4iIfj4mQERE\nZDOJRIKTJ0+iq6sL3d3dyM/Ph0gkQkxMDA4cOACZTIacnBy0tbWhv78ft27dwpMnT5CamgoHBwdI\nJBJERETg4cOHqK2txcDAANRqNWpra+e93qFDh3Dv3j04OjpaNiGwhUajQUZGBurr66HT6dDb24ur\nV69CJBIhMDBwzvimpiaUlpYiNTUVa9euhdFotByfP38GAMTGxkImkyErKwutra3Q6XRoamrC4cOH\nkZuba3OMRES0MBaZ+dlrIiKyQXh4OARBQHBwMK5duwa9Xg9vb2/k5eUhPDwcANDf34+LFy/i1atX\nmJiYwKpVq5CQkICUlBTLOiMjIzh37hwaGxsBAEFBQUhLS0NiYiKUSqXVuzkjIyMICgpCZmYmMjMz\nbY55fHwcJSUlaGhogNFohFgsxrp163D06FHLbnLJyckYHR1FVVUVCgoK8OjRo3nX8vLywosXLwAA\nOp0Oly9fRktLC0wmE6RSKaKiopCdnT3ve1BERPTrMQEiIqLfnlqtRlFRERobG7F06dJfHQ4REf3B\nuAkCERH9tvR6Pdra2qBUKnHs2DEmP0RE9D/jEyAiIvptbdy4EW5uboiLi0NWVpbVt3/evHmDtLS0\nf1wjJiYGhYWFPzNMIiL6gzABIiKiP9Lk5CQMBsM/jpNIJFi2bNkCRERERH8CJkBERERERPTX4DbY\nRERERET012ACREREREREfw0mQERERERE9NdgAkRERERERH+NfwFBdi2Z2Q1wAgAAAABJRU5ErkJg\ngg==\n",
"text/plain": [
"<Figure size 1008x432 with 1 Axes>"
]
},
"metadata": {
"tags": []
}
}
]
},
{
"metadata": {
"id": "BHSH040UsLQT",
"colab_type": "text"
},
"cell_type": "markdown",
"source": [
"北大のスライドではこの後に因子を説明変数に追加してやっていたので一応やってみた"
]
},
{
"metadata": {
"id": "xbWO2r3Sp2NZ",
"colab_type": "code",
"colab": {
"base_uri": "https://localhost:8080/",
"height": 204
},
"outputId": "8d2b59f6-d294-4bea-e928-dd653f3acc28"
},
"cell_type": "code",
"source": [
"X = sm.add_constant(df[['body_size','label']], prepend=False)\n",
"X.head()"
],
"execution_count": 182,
"outputs": [
{
"output_type": "execute_result",
"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>body_size</th>\n",
" <th>label</th>\n",
" <th>const</th>\n",
" </tr>\n",
" </thead>\n",
" <tbody>\n",
" <tr>\n",
" <th>0</th>\n",
" <td>8.31</td>\n",
" <td>0</td>\n",
" <td>1.0</td>\n",
" </tr>\n",
" <tr>\n",
" <th>1</th>\n",
" <td>9.44</td>\n",
" <td>0</td>\n",
" <td>1.0</td>\n",
" </tr>\n",
" <tr>\n",
" <th>2</th>\n",
" <td>9.50</td>\n",
" <td>0</td>\n",
" <td>1.0</td>\n",
" </tr>\n",
" <tr>\n",
" <th>3</th>\n",
" <td>9.07</td>\n",
" <td>0</td>\n",
" <td>1.0</td>\n",
" </tr>\n",
" <tr>\n",
" <th>4</th>\n",
" <td>10.16</td>\n",
" <td>0</td>\n",
" <td>1.0</td>\n",
" </tr>\n",
" </tbody>\n",
"</table>\n",
"</div>"
],
"text/plain": [
" body_size label const\n",
"0 8.31 0 1.0\n",
"1 9.44 0 1.0\n",
"2 9.50 0 1.0\n",
"3 9.07 0 1.0\n",
"4 10.16 0 1.0"
]
},
"metadata": {
"tags": []
},
"execution_count": 182
}
]
},
{
"metadata": {
"id": "Eb7iffvvsZ7l",
"colab_type": "code",
"colab": {
"base_uri": "https://localhost:8080/",
"height": 578
},
"outputId": "671f0fab-8988-4bec-b5e5-f4c397470d07"
},
"cell_type": "code",
"source": [
"poi_log = sm.GLM(df.iloc[:, 0], X, family=sm.families.Poisson(sm.families.links.log))\n",
"poi_log_results = poi_log.fit()\n",
"poi_resid = glm_reslt.resid_deviance\n",
"\n",
"print('Residual of GLM: \\n', poi_resid)\n",
"print('\\n')\n",
"print('Parameters: \\n', poi_log_results.params)\n",
"print('\\n')\n",
"print(poi_log_results.summary())"
],
"execution_count": 183,
"outputs": [
{
"output_type": "stream",
"text": [
"Residual of GLM: \n",
" [-0.23211021 -0.35027122 -0.64396264 -0.22909819 1.06047795 -0.26638437\n",
" -0.23178275 -0.32537884 -0.48538752 0.85555565 -0.22259715 -0.64918082\n",
" -0.88199929 1.81326864 -0.94639849 -0.24758297 -0.3320177 -0.28054444\n",
" -1.33513084 0.91030269 -0.35592175 0.44718924 -0.74400503 -1.95507406\n",
" 0.59395382 1.20963752 0.95233204 -0.85678568 0.58707192 0.33529199\n",
" -1.22731092 2.09663887]\n",
"\n",
"\n",
"Parameters: \n",
" body_size 0.080073\n",
"label -0.031999\n",
"const 1.263105\n",
"dtype: float64\n",
"\n",
"\n",
" Generalized Linear Model Regression Results \n",
"==============================================================================\n",
"Dep. Variable: seeds_num No. Observations: 100\n",
"Model: GLM Df Residuals: 97\n",
"Model Family: Poisson Df Model: 2\n",
"Link Function: log Scale: 1.0\n",
"Method: IRLS Log-Likelihood: -235.29\n",
"Date: Thu, 28 Feb 2019 Deviance: 84.808\n",
"Time: 01:24:16 Pearson chi2: 83.8\n",
"No. Iterations: 4 \n",
"==============================================================================\n",
" coef std err z P>|z| [0.025 0.975]\n",
"------------------------------------------------------------------------------\n",
"body_size 0.0801 0.037 2.162 0.031 0.007 0.153\n",
"label -0.0320 0.074 -0.430 0.667 -0.178 0.114\n",
"const 1.2631 0.370 3.417 0.001 0.539 1.988\n",
"==============================================================================\n"
],
"name": "stdout"
}
]
},
{
"metadata": {
"id": "eCv1HXUBtBHB",
"colab_type": "text"
},
"cell_type": "markdown",
"source": [
"やったあ、こっちも同じ結果になった"
]
},
{
"metadata": {
"id": "ddDng62gs7x_",
"colab_type": "code",
"colab": {}
},
"cell_type": "code",
"source": [
""
],
"execution_count": 0,
"outputs": []
}
]
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment