Created
February 19, 2019 22:27
-
-
Save shotahorii/c4f793a200f6b56683e383eae4352ea8 to your computer and use it in GitHub Desktop.
PoissonRegression
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| { | |
| "cells": [ | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "### Poisson Regression\n", | |
| "- http://hosho.ees.hokudai.ac.jp/~kubo/stat/2013/ou1/kubostat2013ou1.pdf\n", | |
| "- \"データ解析のための統計モデリング入門\" (久保, 2012) - 第2章\n", | |
| "- https://risalc.info/src/st-maximum-likelihood-Poisson.html\n", | |
| "\n", | |
| "ポアソン回帰は、誤差構造がポアソン分布、リンク関数がlog(y)であるGLMである。 \n", | |
| "最尤推定によってポアソン回帰のパラメータ推定を行う。 \n", | |
| "**まず、ポアソン回帰って?** \n", | |
| "データ（の誤差）がポアソン分布に従うと仮定した場合の一般化線形モデル。通常の線形モデルは、誤差を正規分布で仮定したモデルであるため、例えばyが0以下にならないカウントデータの場合などには整合しない場合が多い。 " | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 42, | |
| "metadata": { | |
| "collapsed": false | |
| }, | |
| "outputs": [ | |
| { | |
| "data": { | |
| "text/latex": [ | |
| "ポアソン回帰では、応答変数の誤差構造にポアソン分布を仮定している。これはつまり、データ（の誤差）がyである確率は\n", | |
| "パラメータλによって以下のように決められる、という仮定と等しい。<br><br>\n", | |
| "$$p(y|\\lambda) = \\frac{\\lambda^y exp(-\\lambda)}{y!}$$ <br>\n", | |
| "ここで、λはポアソン分布の唯一のパラメータでλ=平均=分散である。<br>\n", | |
| "また、yは0以上の整数で、全てのyを足し合わせると1になる。つまり以下。<br><br>\n", | |
| "$$\\sum_{y=0}^\\infty p(y|\\lambda)=1$$" | |
| ], | |
| "text/plain": [ | |
| "<IPython.core.display.Latex object>" | |
| ] | |
| }, | |
| "metadata": {}, | |
| "output_type": "display_data" | |
| } | |
| ], | |
| "source": [ | |
| "%%latex\n", | |
| "ポアソン回帰では、応答変数の誤差構造にポアソン分布を仮定している。これはつまり、データ（の誤差）がyである確率は\n", | |
| "パラメータλによって以下のように決められる、という仮定と等しい。<br><br>\n", | |
| "$$p(y|\\lambda) = \\frac{\\lambda^y exp(-\\lambda)}{y!}$$ <br>\n", | |
| "ここで、λはポアソン分布の唯一のパラメータでλ=平均=分散である。<br>\n", | |
| "また、yは0以上の整数で、全てのyを足し合わせると1になる。つまり以下。<br><br>\n", | |
| "$$\\sum_{y=0}^\\infty p(y|\\lambda)=1$$" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 64, | |
| "metadata": { | |
| "collapsed": false | |
| }, | |
| "outputs": [ | |
| { | |
| "data": { | |
| "text/latex": [ | |
| "ポアソン<b>回帰</b>に入る前に、ポアソン分布のパラメータの最尤推定について考える。<br>\n", | |
| "例えばn個のデータがあるとき、パラメータλの尤度は以下のように表される。<br><br>\n", | |
| "$$L(\\lambda)=p(y_1|\\lambda)*p(y_2|\\lambda)*...*p(y_{n}|\\lambda)$$<br>\n", | |
| "$$=\\prod_{i=1}^{n} p(y_i|\\lambda)=\\prod_{i=1}^{n} \\frac{\\lambda^{y_i} exp(-\\lambda)}{y_i!}$$<br>\n", | |
| "$$=exp(-n\\lambda) \\frac{\\lambda^{(y_1+y_2+...+y_n)}}{y_1!y_2!...y_n!}$$<br><br>\n", | |
| "\n", | |
| "この尤度を最大化するパラメータλを求める。ここで、この尤度をこのまま計算するのは大変なので、対数をとる。\n", | |
| "対数尤度最大化は尤度最大化となる。<br><br>\n", | |
| "$$logL(\\lambda) = log(exp(-n\\lambda) \\frac{\\lambda^{(y_1+y_2+...+y_n)}}{y_1!y_2!...y_n!})$$<br>\n", | |
| "$$= -n\\lambda + (y_1+y_2+...+y_n)log\\lambda - log(y_1!y_2!...y_n!)$$<br><br>\n", | |
| "\n", | |
| "この対数尤度を最大化するのは微分して0になる点である。すなわち以下。<br><br>\n", | |
| "$$\\frac{d}{d\\lambda}logL = -n + \\frac{(y_1+y_2+...+y_n)}{\\lambda}$$<br>\n", | |
| "$$0 = -n + \\frac{(y_1+y_2+...+y_n)}{\\lambda}$$<br>\n", | |
| "$$\\lambda = \\frac{(y_1+y_2+...+y_n)}{n}$$<br><br>\n", | |
| "\n", | |
| "つまりポアソン分布におけるλの最尤推定値は標本平均と同様となる。" | |
| ], | |
| "text/plain": [ | |
| "<IPython.core.display.Latex object>" | |
| ] | |
| }, | |
| "metadata": {}, | |
| "output_type": "display_data" | |
| } | |
| ], | |
| "source": [ | |
| "%%latex\n", | |
| "ポアソン<b>回帰</b>に入る前に、ポアソン分布のパラメータの最尤推定について考える。<br>\n", | |
| "例えばn個のデータがあるとき、パラメータλの尤度は以下のように表される。<br><br>\n", | |
| "$$L(\\lambda)=p(y_1|\\lambda)*p(y_2|\\lambda)*...*p(y_{n}|\\lambda)$$<br>\n", | |
| "$$=\\prod_{i=1}^{n} p(y_i|\\lambda)=\\prod_{i=1}^{n} \\frac{\\lambda^{y_i} exp(-\\lambda)}{y_i!}$$<br>\n", | |
| "$$=exp(-n\\lambda) \\frac{\\lambda^{(y_1+y_2+...+y_n)}}{y_1!y_2!...y_n!}$$<br><br>\n", | |
| "\n", | |
| "この尤度を最大化するパラメータλを求める。ここで、この尤度をこのまま計算するのは大変なので、対数をとる。\n", | |
| "対数尤度最大化は尤度最大化となる。<br><br>\n", | |
| "$$logL(\\lambda) = log(exp(-n\\lambda) \\frac{\\lambda^{(y_1+y_2+...+y_n)}}{y_1!y_2!...y_n!})$$<br>\n", | |
| "$$= -n\\lambda + (y_1+y_2+...+y_n)log\\lambda - log(y_1!y_2!...y_n!)$$<br><br>\n", | |
| "\n", | |
| "この対数尤度を最大化するのは微分して0になる点である。すなわち以下。<br><br>\n", | |
| "$$\\frac{d}{d\\lambda}logL = -n + \\frac{(y_1+y_2+...+y_n)}{\\lambda}$$<br>\n", | |
| "$$0 = -n + \\frac{(y_1+y_2+...+y_n)}{\\lambda}$$<br>\n", | |
| "$$\\lambda = \\frac{(y_1+y_2+...+y_n)}{n}$$<br><br>\n", | |
| "\n", | |
| "つまりポアソン分布におけるλの最尤推定値は標本平均と同様となる。" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "ここからポアソン**回帰**の話に入る。ここでは、久保緑本のデータを使う。 \n", | |
| "https://github.com/tushuhei/statisticalDataModeling/blob/master/data3a.csv \n", | |
| "yが応答変数で種子数、xが体サイズ、そしてfが施肥処理の有無を表す(Cが肥料なし、Tが施肥処理)。" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 68, | |
| "metadata": { | |
| "collapsed": false | |
| }, | |
| "outputs": [ | |
| { | |
| "data": { | |
| "text/html": [ | |
| "<div>\n", | |
| "<table border=\"1\" class=\"dataframe\">\n", | |
| " <thead>\n", | |
| " <tr style=\"text-align: right;\">\n", | |
| " <th></th>\n", | |
| " <th>y</th>\n", | |
| " <th>x</th>\n", | |
| " <th>f</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", | |
| " </tbody>\n", | |
| "</table>\n", | |
| "</div>" | |
| ], | |
| "text/plain": [ | |
| " y x f\n", | |
| "0 6 8.31 C\n", | |
| "1 6 9.44 C" | |
| ] | |
| }, | |
| "execution_count": 68, | |
| "metadata": {}, | |
| "output_type": "execute_result" | |
| } | |
| ], | |
| "source": [ | |
| "import pandas as pd\n", | |
| "df = pd.read_csv('./data3a.csv')\n", | |
| "df.head(2)" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 70, | |
| "metadata": { | |
| "collapsed": false | |
| }, | |
| "outputs": [ | |
| { | |
| "data": { | |
| "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAEACAYAAACj0I2EAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAD8BJREFUeJzt3X+sZHdZx/H3Z1ktlIZa0e4VFrvUBCQNWBqtaKMZfoVG\nA200QcHEFowhRmgDBlvqH3s1xuBGMUTlD6Q0C6EqVKHFgF2aMiHVtCB06QJrxWC3rbCXX0VsmmyA\n+/jHTOt22bvze87O975fyWTPPfec8zx778znfud7Zs6kqpAkrb4dXTcgSZoPA12SGmGgS1IjDHRJ\naoSBLkmNMNAlqREjAz3JGUnuSnJ3kkNJ9g7Xn5PkQJJ7k9ya5OzFtytJ2krGeR16kjOr6pEkTwD+\nBbgK+FXgG1W1L8k1wDlVde1i25UkbWWsKZeqemS4eAawEyjgMmD/cP1+4PK5dydJGttYgZ5kR5K7\ngaPAx6rqU8CuqtoAqKqjwLmLa1OSNMq4I/TNqno+sBu4OMkFDEbpj9ts3s1Jksa3c5KNq+rbSfrA\npcBGkl1VtZFkDfjqyfZJYtBL0hSqKpNsP86rXH7k0VewJHkS8FLgMHALcOVwsyuAm0/R1Mre9u7d\n23kP27X/Ve7d/ru/rXr/0xhnhP5jwP4kOxj8Afj7qvpIkjuB9yd5LXAEeOVUHUiS5mJkoFfVIeCi\nk6z/JvCSRTQlSZqc7xQdodfrdd3CTFa5/1XuHey/a6ve/zTGemPRTAWSWnQNSWpNEmreJ0UlSavB\nQJekRhjoktQIA12SGmGgS1IjDHRJaoSBLkmNMNAlqREGeoPW1vaQpJPb2tqerv/70rblO0UblITu\nLk+fqa8UJ+n/+U5RSdrGDHRJaoSBLkmNMNAlqREGuiQ1wkCXpEYY6JLUCANdkhphoEtSIwx0SWqE\ngS5JjTDQJakRBrokNcJAl6RGGOiS1AgDXZIaMTLQk+xOcnuSzyc5lOQNw/V7kzyY5DPD26WLb1eS\ntJWRn1iUZA1Yq6qDSc4CPg1cBvwa8L9V9bYR+/uJRUvmJxZJq2+aTyzaOWqDqjoKHB0uP5zkMPD0\nR2tO3KUkaSEmmkNPsge4ELhruOr1SQ4meVeSs+fcmyRpAmMH+nC65Sbg6qp6GHgHcH5VXchgBH/K\nqRdJ0mKNnHIBSLKTQZi/t6puBqiqrx23yd8AH95q//X19ceWe70evV5vilYlqV39fp9+vz/TMUae\nFAVI8h7g61X1puPWrQ3n10nyRuBnqurVJ9nXk6JL5klRafVNc1J0nFe5XAJ8AjjEICUKuA54NYP5\n9E3gPuB1VbVxkv0N9CUz0KXVt5BAn5WBvnwGurT6pgl03ykqSY0w0CWpEQa6JDXCQJekRhjoktQI\nA12SGmGgS1IjDHRJaoSBLkmNMNAlqREGuiQ1wkCXpEYY6JLUCANdkhphoEtSIwx0SWqEgS5JjTDQ\nJakRBrokNcJAl6RGGOiS1AgDXZIaYaBLUiMMdElqhIEuSY0w0CWpEQa6JDXCQJekRhjoktSIkYGe\nZHeS25N8PsmhJFcN15+T5ECSe5PcmuTsxbcrSdpKqurUGyRrwFpVHUxyFvBp4DLgNcA3qmpfkmuA\nc6rq2pPsX6NqaL6SAF39zIO/b2l2SaiqTLLPyBF6VR2tqoPD5YeBw8BuBqG+f7jZfuDyydqVJM3T\nRHPoSfYAFwJ3AruqagMGoQ+cO+/mJEnj2znuhsPplpuAq6vq4SQnPq/e8nn2+vr6Y8u9Xo9erzdZ\nl5LUuH6/T7/fn+kYI+fQAZLsBP4J+GhVvX247jDQq6qN4Tz7x6vqOSfZ1zn0JXMOXVp9C5lDH3o3\n8IVHw3zoFuDK4fIVwM2TFJYkzdc4r3K5BPgEcIjBsK+A64BPAu8HngEcAV5ZVd86yf6O0JfMEbq0\n+qYZoY815TILA335DHRp9S1yykWSdJoz0CWpEQa6NAdra3tI0tltbW1P1z8CnQacQ2+Qc+jL1+3P\nHLbrz71lzqFL0jZmoEtSIwx0SWqEgS5JjTDQJakRBrokNcJAl6RGGOiS1AgDXZIaYaBLUiMMdElq\nhIEuSY0w0CWpEQa6JDXCQJekRhjoktQIA12SGmGgS1IjDHRJaoSBLkmNMNAlqREGuiQ1wkCXpEYY\n6JLUiJGBnuT6JBtJ7jlu3d4kDyb5zPB26WLblCSNMs4I/QbgZSdZ/7aqumh4++c59yVJmtDIQK+q\nO4CHTvKtzL8dSdK0ZplDf32Sg0neleTsuXUkSZrKzin3ewfwR1VVSf4YeBvwW1ttvL6+/thyr9ej\n1+tNWVba2traHjY2jnTdhjSVfr9Pv9+f6RipqtEbJecBH66q503yveH3a5wamp8kQFc/89DV77vr\n/3d3tQf1fZy1JQlVNdHU9rhTLuG4OfMka8d971eAz01SVJI0fyOnXJLcCPSApya5H9gLvDDJhcAm\ncB/wugX2KEkaw1hTLjMVcMpl6bqeenDKpZv6Ps7assgpF0nSac5Al6RGGOiS1AgDXZIaYaBLUiMM\ndElqhIEuSY0w0CWpEQa6JDXCQJekRhjoktQIA12SGmGgS1IjDHRJaoSBLkmNMNAlqREGuiQ1wkCX\npEaM/ExRaTJnDD8KTtKyGeias2N0+7me0vbllIskNcJAl6RGGOiS1AgDXZIaYaBLUiMMdElqhIEu\nSY0w0CWpESMDPcn1STaS3HPcunOSHEhyb5Jbk5y92DYlSaOMM0K/AXjZCeuuBW6rqmcDtwNvmXdj\nkqTJjAz0qroDeOiE1ZcB+4fL+4HL59yXJGlC086hn1tVGwBVdRQ4d34tSZKmMa+Lc53yakzr6+uP\nLfd6PXq93pzKSlIb+v0+/X5/pmOkavSV8ZKcB3y4qp43/Pow0KuqjSRrwMer6jlb7Fvj1ND8DC5f\n2+UVD63dRX0fZ21JQlVNdAnRcadcwuOvTXoLcOVw+Qrg5kmKSpLmb+QIPcmNQA94KrAB7AU+BHwA\neAZwBHhlVX1ri/0doS+ZI/TtVntQ38dZW6YZoY815TILA335DPTtVntQ38dZWxY55SJJOs0Z6JLU\nCANdkhrhh0QvyNraHjY2jnTdhqRtxJOiC+KJSWsvu/52fJy1zJOikrSNGeiS1AgDXZIaYaBLUiMM\ndElqhIEuSY0w0CWpEQa6JDXCQJekRhjoktQIA12SGmGgS1IjDHRJaoSBLkmNMNAlqREGuiQ1wkCX\npEYY6JLUCANdkhphoEtSIwx0SWqEgS5JjTDQJakRO2fZOcl9wP8Am8B3qurieTQlSZrcTIHOIMh7\nVfXQPJqRJE1v1imXzOEYkqQ5mDWMC/hYkk8l+e15NCRJms6sUy6XVNVXkvwog2A/XFV3zKMxSdJk\nZgr0qvrK8N+vJfkgcDHwfYG+vr7+2HKv16PX681Sdmxra3vY2DiylFqSNIt+v0+/35/pGKmq6XZM\nzgR2VNXDSZ4MHAD+sKoOnLBdTVtjVkkYzAp1Ut3a1l5q/a4eZ1qMJFRVJtlnlhH6LuCDSWp4nPed\nGOaSpOWZeoQ+dgFH6Na29lLqO0JvyzQjdF9yKEmNMNAlqRGzvmxxLPv27VtGGUna1pYyh75z55sX\nWuNkNjfvYnPzE2zPOVVrb6/aAE8EjnVSedeu8zh69L5Oardsmjn0pQR6N3f0fcA1bM8HuLW3V+2u\n63tCdhE8KSpJ25iBLkmNMNAlqREGuiQ1wkCXpEYY6JLUCANdkhphoEtSIwx0SWqEgS5JjTDQJakR\nBrokNcJAl6RGLOV66JJadsbw4x6Xb8eOM9ncfKST2qfjZYMNdEkzOkZXl+7d3OzussEbG938ETsV\np1wkqREGuiQ1wkCXpEYY6JLUCANdkhphoEtSIwx0SWqEgS5JjZgp0JNcmuTfk/xHkmvm1ZQkaXJT\nB3qSHcBfAS8DLgBeleQn59XY6aPfdQMz6nfdwAz6XTcwo37XDcyo33UDM+p33cDSzTJCvxj4YlUd\nqarvAH8HXDaftk4n/a4bmFG/6wZm0O+6gRn1u25gRv2uG5hRv+sGlm6WQH868MBxXz84XCdJ6sBS\nLs71lKe8fBllHufYsf/k2LGll5WkzqRquiuVJXkBsF5Vlw6/vhaoqvrTE7br5lJokrTiqmqiSzrO\nEuhPAO4FXgx8Bfgk8KqqOjzVASVJM5l6yqWqvpfk9cABBnPx1xvmktSdqUfokqTTy8LeKZpkd5Lb\nk3w+yaEkVy2q1qIk2ZHkM0lu6bqXSSU5O8kHkhwe/g5+tuueJpHkjUk+l+SeJO9L8oNd93QqSa5P\nspHknuPWnZPkQJJ7k9ya5OwuezyVLfrfN7z/HEzyD0me0mWPp3Ky/o/73u8l2Uzyw130No6t+k/y\nhuHv4FCSt446ziLf+v9d4E1VdQHwc8DvruAbj64GvtB1E1N6O/CRqnoO8FPAykyHJXka8Abgoqp6\nHoOpwV/vtquRbmDwJrvjXQvcVlXPBm4H3rL0rsZ3sv4PABdU1YXAF1m9/kmyG3gpcGTpHU3m+/pP\n0gNeDjy3qp4L/Nmogyws0KvqaFUdHC4/zCBQVuZ16sM7wi8B7+q6l0kNR1K/UFU3AFTVd6vq2x23\nNaknAE9OshM4E/hyx/2cUlXdATx0wurLgP3D5f3A5UttagIn67+qbquqzeGXdwK7l97YmLb4+QP8\nBfDmJbczsS36/x3grVX13eE2Xx91nKVcnCvJHuBC4K5l1JuTR+8Iq3iS4ZnA15PcMJwyemeSJ3Xd\n1Liq6svAnwP3A/8NfKuqbuu2q6mcW1UbMBjgAOd23M8sXgt8tOsmJpHkFcADVXWo616m9CzgF5Pc\nmeTjSX561A4LD/QkZwE3AVcPR+qnvSS/DGwMn2FkeFslO4GLgL+uqouARxg8/V8JSX6Iwej2POBp\nwFlJXt1tV3OxioMDkvwB8J2qurHrXsY1HMBcB+w9fnVH7UxrJ3BOVb0A+H3g/aN2WGigD58u3wS8\nt6puXmStObsEeEWSLwF/C7wwyXs67mkSDzIYmfzb8OubGAT8qngJ8KWq+mZVfQ/4R+DnO+5pGhtJ\ndgEkWQO+2nE/E0tyJYOpx1X7g/oTwB7gs0n+i8F00aeTrNKzpAcY3Pepqk8Bm0meeqodFj1Cfzfw\nhap6+4LrzFVVXVdVP15V5zM4GXd7Vf1m132Na/g0/4EkzxquejGrdXL3fuAFSZ6YJAz6X4WTuic+\nm7sFuHK4fAVwug9qHtd/kksZTDu+oqpW4UIaj/VfVZ+rqrWqOr+qnslgkPP8qjqd/6ieeP/5EPAi\ngOFj+Qeq6hunOsAiX7Z4CfAbwIuS3D2cy710UfX0fa4C3pfkIINXufxJx/2Mrao+yeBZxd3AZxnc\nyd/ZaVMjJLkR+FfgWUnuT/Ia4K3AS5M8+o7qkS8768oW/f8lcBbwseHj9x2dNnkKW/R/vOI0nnLZ\nov93A+cnOQTcCIwcVPrGIklqhB9BJ0mNMNAlqREGuiQ1wkCXpEYY6JLUCANdkhphoEtSIwx0SWrE\n/wHVo4k4FTCGRQAAAABJRU5ErkJggg==\n", | |
| "text/plain": [ | |
| "<matplotlib.figure.Figure at 0x10755a0b8>" | |
| ] | |
| }, | |
| "metadata": {}, | |
| "output_type": "display_data" | |
| } | |
| ], | |
| "source": [ | |
| "from matplotlib import pyplot as plt\n", | |
| "%matplotlib inline\n", | |
| "fig, ax = plt.subplots(1,1)\n", | |
| "chart = ax.hist(df['y'].tolist())" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 71, | |
| "metadata": { | |
| "collapsed": false | |
| }, | |
| "outputs": [ | |
| { | |
| "data": { | |
| "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAEACAYAAACj0I2EAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAGdNJREFUeJzt3X2QHPdd5/H3V6sddmRZa9mW7UtkdmIbY1fiJ10iUsfT\nGOKKeXKA4zhiHi65448cKE7lqFQCXKGtu6urOFUUR8HlD0AIh8rKwSGEpAjHksIDFShHIrGzTrR2\nUiSzliHRjh3HF/m2rI30vT9mRjs72zPT3dM9PfObz6tqSjs9/fDtntF3e3/b/Vlzd0REZPLtKroA\nERHJhhq6iEgg1NBFRAKhhi4iEgg1dBGRQKihi4gEYmBDN7NjZnbWzFa6pr/dzFbN7Ekze29+JYqI\nSBy7Y8xzHPgd4APtCWZWBX4MuM3dv2VmV+dTnoiIxDXwDN3dPwW80DX5PwPvdfdvteZ5LofaREQk\ngbRj6DcD32dmj5nZo2b22iyLEhGR5OIMufRabr+7v97MXgf8CXBDdmWJiEhSaRv6GeAjAO5+yswu\nmtlV7v5894xmprAYEZEU3N2SzB93yMVaj7aPAj8AYGY3A7NRzbyjqIl9HD16tPAaprX+Sa5d9Rf/\nmPT60xh4hm5mS0AVuMrMngGOAn8IHDezJ4GXgV9ItXUREcnMwIbu7vf3eOnnM65FRESGoDtFB6hW\nq0WXMJRJrn+SawfVX7RJrz8NSztWE3sDZp73NkREQmNmeE6/FBURkTGnhi4iEgg1dBGRQKihi4gE\nQg1dRCQQaugiPTQaDU6dOkWj0Si6FJFY1NBFIpw48SEWFm7hnnvexsLCLZw48aGiSxIZSNehi3Rp\nNBosLNzCxsajwO3ACuXy3aytPcWBAweKLk+mhK5DF8lAvV6nVKrQbOYAtzM7u0C9Xi+uKJEY1NBF\nulQqFc6frwPtP6O7wubmGpVKpbiiRGJQQxfpcuDAAY4dez/l8t3s23eIcvlujh17v4ZbZOxpDF2k\nh0ajQb1ep1KpqJnLyKUZQ1dDFxEZQ/qlqIjIFFNDFxEJhBq6iEgg1NBFRAIxsKGb2TEzO2tmKxGv\n/YqZXTSzK/MpT0RE4opzhn4ceGP3RDM7CNwDrGVdlIiIJDewobv7p4AXIl76LeBdmVckEjilOEpe\nUo2hm9l9wBl3fzLjekSCphRHyVOsG4vMbAH4uLvfbmZl4FHgHnf/ppl9BXituz/fY1ndWCSCUhwl\nmTQ3Fu1OsZ0bgQrwOTMz4CDwGTM77O7rUQssLi5e+rparVKtVlNsVmSytVMcNzZ2pjiqoUutVqNW\nqw21jrhn6BWaZ+i3Rbz2FeCQu0eNs+sMXaRFZ+iSRC63/pvZEvAPwM1m9oyZvbVrFgcSbVRkGinF\nUfKmcC6REVOKo8ShtEURkUAobVFEZIqpoYuIBEINXUQkEGroIiKBUEMXEQmEGroESyFY4dN7vJ0a\nugRJIVjh03u8k65Dl+DoFvvwTcN7rOvQRdgKwWr+R4fOECwJg97jaGroEpxKpcL583Wg/VcTV9jc\nXKNSqRRXlGRK73E0NXQJjkKwwqf3OJrG0CVYCsEKX8jvscK5REQCoV+KiohMMTV0EZFAqKGLiARC\nDV1EJBBq6CIigYjzR6KPmdlZM1vpmPY+M1s1syfM7E/NbF++ZYqIyCBxztCPA2/smrYMvNrd7wS+\nBPxq1oWJdFKq3mhkcZzzfK/0OehvYEN3908BL3RN+6S7X2w9fQw4mENtIoBS9UYli+Oc53ulz0EM\n7j7wASwAKz1e+xhwf59lXSSt9fV1L5evdPicgzt8zsvlK319fb3o0oKSxXHO872axs9Bq3fG6tHt\nx+5hvhmY2a8Dm+6+1G++xcXFS19Xq1Wq1eowm5Up0k7V29jYmaoX2q3eRcriOOf5Xk3D56BWq1Gr\n1YZbSZyuT8QZOvAW4O+BbxuwbO7fySRc03hmVgSdoY8fUpyhx71s0VqP5hOze4F3Afe5+8vDfUsR\n6U2peqORxXHO873S5yCegeFcZrYEVIGrgLPAUeDXgBLwfGu2x9z9l3os74O2ITJIyKl64ySL45zn\nezVNnwOlLYqIBEJpiyIiU0wNXUQkEGroIiKBUEMXEQmEGrqISCDU0EVEAqGGLlOh0WiwvLzM8vJy\n4qS+cUn4y6KO1dVVHnroIVZXVzOsrGlcjtNUS3pradIHuvVfCra09LCXSvMONzns8dnZvb609HDs\nZcvlK31+/pCXy1fGXi5rWdRx5Mg7HMoONzuU/ciRB8aqPtmOFLf+q6FL0KIyQGC/z81dMTAHZFzy\nQ7Ko4/Tp061m3nkcyn769OmxqE92StPQNeQiQavX6+zadT2wldIHFWZmrqFerw9ctlSqbFu2nfA3\nSlnUcfLkSaD7OBxsTS++PsmGGroErVKpcPHiGaD9FxRXgDoXLqxTqVQGLnv+fH3bspubawOXy1oW\ndRw+fBjoPg7PtqYXX59kJOkpfdIHGnKRgm2Nod+Yegx93767xmIMfZg6jhx5oDXs8h25jaEXfZxC\nQoohF4VzyVRoNBo8/vjjANx1112JkvrGJeEvizpWV1c5efIkhw8f5tZbbx27+mSL0hZFRAKhtEUR\nkSmmhi4iEgg1dBGRQKihi4gEQg1dRCQQAxu6mR0zs7NmttIxbb+ZLZvZ02b2V2Y2n2+ZIiIyyMDL\nFs3se4BzwAfc/fbWtAeB5939fWb2bmC/u7+nx/K6bFEiRV23nPR68VFd+5zFdtrr2Lt3L2fOnAF6\n72O/7XW+Bkz0td+6dr23NJctxr3bcwFY6Xj+FHBt6+vrgKf6LJvHTVQy4aLS+ZKmIo4q4S+L7bTX\nUS7f1rpb87qe+9hve52vzc5e7qXS/MQmHCqhsT/ySluMaOhf73r9632WzX3HZbJEpfPNzV3hc3P7\nY6cijirhL4vtRCc+Xunw6I597Le97a+tO+zPff/zooTGwdI09N1D/1zQOtHv9+Li4uKlr6vVKtVq\nNaPNyiRqp/NtbGyl883MXMPFiyV2piK+RL1e3/HjeNQ62gl/Wf7onsV2otbRPEe6jO597Lc9oOO1\nU8CriEo4nIShi1G9f5OkVqtRq9WGW0mcrs/OM/RVtg+5rPZZdgTfy2SS6AxdZ+g6Qx+MHIdcKsCT\nHc8fBN7d+vrdwHv7LDuCXZdJE5XOlzQVcVQJf1lsZ2sM/TWtMfRrB46hR22v87XZ2b1eKs1PbMKh\nEhr7S9PQ41zlsgRUgauAs8BR4KPAIzQT89eAn3b3b/RY3gdtQ6aTrnLRVS66yqU3pS2KiARCaYsi\nIlNMDV1EJBBq6CIigVBDFxEJhBq6iEgg1NBFRAKR1a3/In0lvb58FIa9Brp7+c7rzM+dO7djevf1\n9lHLdj5vH6/rr7+ec+fO7Vhv1vsziaZxn/tKeidS0ge6U3TqLS097LOzlzvscbjJS6X5wu8KHDbp\nr3v5I0ceaN0JeoND2cvl27ZN706VjFq28/nW8fpXDmUvlW5trfdVkfVOY3Jh6PtMXrf+D/NQQ59u\n6+vrPjd3xVjljgybIxKdzVJ2+LNWRkv39Ee3ZdZEL/to1/P51rToDJjOeqcxF2Ua9jlNQ9cYuuSq\nXq8zM3Mt3cmAu3YdvJQgWERNpVKFqKTCtMvDQeBFmrFH3dMvu/R8ZuYadu26vu888EqamXeXRayv\nmdLYWe+w+zOJpnGf41BDl1xVKhUuXDgLfAVo/xXDFS5efPZSFkkRNZ0/X99Wz+bmWux6opaHZ4F5\nIGr6S5eeX7iwzsWLZ/rOA/8MfK01rXt9a8BL2+oddn8m0TTucyxJT+mTPtCQy9RrjqHvbY0J3zhW\nY+hpk/66l2+Pg8/NVVpj3a/ZNr07VTJq2c7nW8fr2m1j6HNzlb5j6NOUXBj6PpNH2uKwFM4loKtc\ndJVLPkLeZ6UtiogEQmmLIiJTTA1dRCQQaugiIoFQQxcRCYQauohIIIZq6Gb2TjP7vJmtmNkHzayU\nVWEiIpJM6oZuZq8A3g4ccvfbaSY3/kxWhUlYGo0Gp06dotFojOX6itZoNFheXmZ5ebmwfeo8plHH\nN84xT/q+xJ0/j23HNVGftaR3IvnWHaCvoHkf8n6azfzjwBsi5svvViqZCFmn4oWWsre09LCXSvMO\nNzns8dnZvSPfp85jOjt7uZdK830TIqPqS/q+xJ0/j23HVeRnjVGnLQIPAN8EzgJ/3GOe3HdcxlfW\nqXihpexFJzfu97m5K0a2T9trWI9MxmwmZvY+5knfl7jzx5kvr89E0Z+1NA099R+4MLMrgDfRjH97\nEfiwmd3v7kvd8y4uLl76ulqtUq1W025WJkw7FW9jY2cqXppbtbNeX9Hq9XpE+mKFmZmXRrZP24/p\nKaKSMeFlopIN2/UlfV/izh9nvrw+E6P+rNVqNWq12nArSfodoP0Afgr4/Y7nPw/8bsR8eX8jkzGm\nM/T+dIauM/ReGOWQC3AYeBKYAwz4I+CXI+Ybwa7LOMs6FS+0lL2tMfQbCx9D37fvLp+d3eul0nzf\nhMh+49hx35e48+ex7biK/KylaehDhXOZ2VGaV7ZsAo8Dv+jum13z+DDbkDBknYoXWsreOKRRdh5T\nYGBC5KB1xNmHuPPnse24ivqsKW1RRCQQSlsUEZliaugiIoFQQxcRCYQauohIINTQRUQCoYYuIhKI\n1Lf+i0yiPK4pTnMNeb9lsqixyOu2pUBJ70RK+kB3isqYyCM5b2npYZ+dvdxhj8NNXirNx7pDstcy\nWdRYZDqhZIdRpy3G2oAauoyBPHI51tfXWxknO7NP+mWY9Frm9OnTQ9dYZPaJZCtNQ9cYukyFdnJe\nVGLgMOucmbmWqHTCXuvtt8zJkyeHrjHOfuZxLGQ8qKHLVKhUKpw/XwdWWlNW2Nxcu5RbknadFy6c\nBb6ybb0XLz7bc739ljl8+PDQNcbZzzyOhYyJpKf0SR9oyEXGRB7Jec3x8L2t8fAbE4yhRy+TRY1F\nphNKdhh12mIcCueScaKrXJLNI8VR2qKISCCUtigiMsXU0EVEAqGGLiISCDV0EZFAqKGLiARiqIZu\nZvNm9oiZrZrZF8zsu7IqTEREkhk2bfG3gU+4+78zs93AngxqkoJl+ZfYR1VPr3k6pz/33HOcPHmS\nw4cPc/XVV4/FX4hvz793717OnTu3bbmodXVOA3oum3fdwy43KeubuGv1k96J5Ft3gO4D/inGfPnc\nRiW5iJvCN6q0vmGSAzunz8xc5lB2uNmh7DMzezKvPekxac9fLt/gUPZy+bZLy0Wtq3Pa7OzlXirN\ne7l8W2vZV+WazpjlcpOyvqITKRll2iJwB/Bp4DjwWeD3gHLEfCPYdclC3BS+UaX1DZMcuD25cH1H\numHz+XpmtSc9JlvzP+qwc7lmIuPWtLm5KwbsT3NdeaQzZrncpKxvHBIp0zT0YYZcdgOHgF929380\ns/8FvAc42j3j4uLipa+r1SrVanWIzUpe2il8Gxs7U/g6f9yMO98o6uk1Tzu5sDn9FN3phlAB6sDr\nMqk96THZmv+yVi3bkxfh5W3TZmauAcqtaVH7swBclnhf0r6XWX8Gxm19o/qMd6rVatRqteFWkvQ7\nQPsBXAt8ueP59wAfj5gv729kkhGdoedba/T8OkMfx/VN6hl66obe3B5/C9zc+voo8GDEPLnvuGQn\nbgrfqNL6hkkO7Jw+M7PHm2Po3+HtMfSsa096TNrzz81VvDkO/pod4+Wd6+qcNju7tzWG/hqHss/N\nVXJNZ8xyuUlZX9GJlGka+lDhXGZ2B/AHwCzwZeCt7v5i1zw+zDZk9HSVS761Rs2vq1zGc31FXuWi\ntEURkUAobVFEZIqpoYuIBEINXUQkEGroIiKBUEMXEQmEGrqISCCGTVucehOXxpahftdQj2L5LOR9\nnfKwr4skkvROpKQPAr5TtOg0tiJtJQWmS/sbdvks5J3GN+zrMt0Y9a3/sTYQaEMfh6yHokTte5Is\nkWGXz2sfssz6GPZ1kTQNXWPoKbXT2DoT79ppbKGL2vfutL88l8/CsO/foOWHfV0kDTX0lCqVCufP\n14GV1pQVNjfXLmVshCxq32ENeCnWMRh2+SwM+/4NWn7Y10VSSXpKn/RBoEMu7sWnsRVpaww8Xdrf\nsMtnIe80vmFfl+nGqNMW4wg9nGuar1LQVS66ykXyo7RFEZFAKG1RRGSKqaGLiARCDV1EJBBq6CIi\ngVBDFxEJxNAN3cx2mdlnzexjWRQkIiLpZHGG/g7gdAbrkTHSaDQ4deoUjUaj6FJGpnOfB+3/JB2f\nSapVhpT0TiTffhfoQeCvgSrwsR7z5HYnleRjGlMAO/e5VJr32dm9QaQkTlKtsh2jTlsEHgHuBL5f\nDT0M05gCGJ3+uN9hfaJTEiepVtkpTUNP/QcuzOxHgLPu/oSZVYGedzQtLi5e+rparVKtVtNuVnLW\nTgHc2NiZAhjqrelR+wwVoA68btv+T9LxmaRaBWq1GrVabah1pL7138z+J/BzwLeAMnA58BF3/4Wu\n+TztNmT0Go0GCwu3sLHxKM3GtkK5fDdra08F2wSi9rk5ivg08NVt+z9Jx2eSapWd0tz6n1WiooZc\nAjKNKYCd+9weQw8hJXGSapXtKCpt0cy+H/gVd78v4jXPYhsyWtOYAti5z0AwKYmTVKtsUdqiiEgg\nlLYoIjLF1NBFRAKhhi4iEgg1dBGRQKihi4gEQg1dRCQQauiBU9LeFh2L3nRswqCGHrATJz7EwsIt\n3HPP21hYuIUTJz5UdEmF0bHoTccmHLqxKFDK8diiY9Gbjs340o1Fckk7aa/5nxQ6k/amjY5Fbzo2\nYVFDD1SlUuH8+TrN5ECAFTY31y7llEwTHYvedGzCooYeqAMHDnDs2Pspl+9m375DlMt3c+zY+6fy\nx2gdi950bMKiMfTAKWlvi45Fbzo240dpiyIigdAvRUVEppgauohIINTQRUQCoYYuIhKI1A3dzA6a\n2d+Y2RfM7EkzeyDLwkREJJnUV7mY2XXAde7+hJntBT4DvMndn+qaT1e5iIgkNNKrXNz9a+7+ROvr\nc8Aq8Mq06xPJmxIFJXSZjKGbWQW4E/h0FusTyZoSBWUaDH1jUWu4pQb8d3f/84jXNeQihVKioEyi\nNEMuu4fc4G7gw8AfRzXztsXFxUtfV6tVqtXqMJsVSaSdKLixsTNRUA1dxkWtVqNWqw21jqHO0M3s\nA8Bz7v5f+syjM3QplM7QZRKN9JeiZvbdwM8CP2Bmj5vZZ83s3rTrE8mLEgVlWiicS6aGEgVlkiht\nUUQkEEpbFBGZYmroIiKBUEMXEQmEGrqISCDU0EVEAqGGLiISCDV0EZFAqKGLiARCDV1EJBBq6CIi\ngVBDFxEJhBq6iEgg1NBFRAKhhi4iEgg1dBGRQKihi4gEQg1dRCQQaugiIoEYqqGb2b1m9pSZfdHM\n3p1VUSIiklzqhm5mu4DfBd4IvBp4s5ndklVh46JWqxVdwlAmuf5Jrh1Uf9Emvf40hjlDPwx8yd3X\n3H0TeBh4UzZljY9J/1BMcv2TXDuo/qJNev1pDNPQXwmc6Xj+bGuaiIgUQL8UFREJhLl7ugXNXg8s\nuvu9refvAdzdH+yaL90GRESmnLtbkvmHaegzwNPADwJfBU4Cb3b31VQrFBGRoexOu6C7XzCzI8Ay\nzaGbY2rmIiLFSX2GLiIi4yXXX4qa2byZPWJmq2b2BTP7rjy3lxUzu9nMHjezz7b+fdHMHii6riTM\n7J1m9nkzWzGzD5pZqeiakjCzd5jZk63H2B97MztmZmfNbKVj2n4zWzazp83sr8xsvsga++lR/0+1\nPkMXzOxQkfUN0qP+97V6zxNm9qdmtq/IGvvpUf9/M7PPtXrQ/zGz6watJ++rXH4b+IS73wrcAUzE\nkIy7f9Hd73L3Q8C/Bl4C/qzgsmIzs1cAbwcOufvtNIfWfqbYquIzs1cD/wl4LXAn8KNmdkOxVQ10\nnOZNdp3eA3zS3b8T+BvgV0deVXxR9T8J/ATwt6MvJ7Go+peBV7v7ncCXmLzj/z53v8Pd7wL+Ajg6\naCW5NfTWd8PvdffjAO7+LXf/v3ltL0dvAP7J3c8MnHO8zACXmdluYA/wLwXXk8StwKfd/WV3vwD8\nHfCTBdfUl7t/Cniha/KbgIdaXz8E/PhIi0ogqn53f9rdvwQkutKiCD3q/6S7X2w9fQw4OPLCYupR\n/7mOp5cBFxkgzzP0VwHPmdnx1tDF75lZOcft5eXfAyeKLiIJd/8X4DeBZ4B/Br7h7p8stqpEPg98\nb2vIYg/ww8D1BdeUxjXufhbA3b8GXFNwPdPsPwJ/WXQRSZnZ/zCzZ4D7gd8YNH+eDX03cAj4362h\ni/9H80fQiWFms8B9wCNF15KEmV1B8+xwAXgFsNfM7i+2qvjc/SngQeCvgU8AjwMXCi0qG7oCoQBm\n9uvAprsvFV1LUu7+X93924EP0hxG7SvPhv4scMbd/7H1/MM0G/wk+SHgM+7eKLqQhN4AfNndv94a\nsvgI8G8KrikRdz/u7q919yrwDeCLBZeUxlkzuxag9Qut9YLrmTpm9haaP+FNzAlND0vAvx00U24N\nvfWj5hkzu7k16QeB03ltLydvZsKGW1qeAV5vZnNmZjSP/UT8QrrNzA60/v12mr+Ym4SzK2P7ePPH\ngLe0vv4PwJ+PuqCEuuvvfm3cbavfzO4F3gXc5+4vF1ZVfN3139Tx2o8T4/9wrtehm9kdwB8As8CX\ngbe6+4u5bTBDrbHbNeAGd/9m0fUkZWZHaV7ZsklzyOIXW6mYE8HM/g64kmb973T3WrEV9WdmS0AV\nuAo4S/OKhI/SHK67nuZn6afd/RtF1dhPj/pfAH4HuJrmT0lPuPsPFVVjPz3q/zWgBDzfmu0xd/+l\nQgocoEf9PwJ8J83hxjXgbe7+1b7r0Y1FIiJhUNqiiEgg1NBFRAKhhi4iEgg1dBGRQKihi4gEQg1d\nRCQQaugiIoFQQxcRCcT/B6L6GDTRUgUXAAAAAElFTkSuQmCC\n", | |
| "text/plain": [ | |
| "<matplotlib.figure.Figure at 0x10a401e10>" | |
| ] | |
| }, | |
| "metadata": {}, | |
| "output_type": "display_data" | |
| } | |
| ], | |
| "source": [ | |
| "fig, ax = plt.subplots(1,1)\n", | |
| "chart = ax.scatter(df['x'].tolist(), df['y'].tolist())" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "まずはfの値は無視して、xのみを考慮してモデリングする。 \n", | |
| " \n", | |
| "このデータ(y,x)に対して、ポアソン回帰を行うというのはどういうことだろうか。 \n", | |
| "yデータの分布(上のヒストグラム)がポアソン分布に従っている、と仮定してみる。 \n", | |
| "すると、下の散布図において、こう考えることができる。「yデータ(y軸)は全体としてポアソン分布に従って分布している。では、これをx方向も考慮して考えてみると、各x点において、それぞれ特有の平均λのポアソン分布に従ってyが分布しているのではないか。」 \n", | |
| "ということで、これを数式で表してみる。" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 82, | |
| "metadata": { | |
| "collapsed": false | |
| }, | |
| "outputs": [ | |
| { | |
| "data": { | |
| "text/latex": [ | |
| "まず、「yデータは全体としてポアソン分布に従って分布している」という部分。これは本記事の最初に書いたように、以下の式で表される。<br><br>\n", | |
| "$$p(y|\\lambda) = \\frac{\\lambda^y exp(-\\lambda)}{y!}$$<br><br>\n", | |
| "\n", | |
| "ここで、「各x点においてそれぞれ特有の平均λのポアソン分布に従ってyが分布している。」という点を考慮して式を書き直す。つまり各x点でのy_iは、\n", | |
| "そのx点でのλ_iに従う、というように書き直す。<br><br>\n", | |
| "$$p(y_i|\\lambda_i) = \\frac{\\lambda_i^{y_i} exp(-\\lambda_i)}{y_i!}$$<br><br>\n", | |
| "\n", | |
| "そしてここで、このように考えてみる。yは種子数で、これはポアソン分布に従ってばらつき分布している。この分布の平均の違いは、体サイズxによって\n", | |
| "表すことができる、と。これを以下のような式で表してみる。<br><br>\n", | |
| "\n", | |
| "$$\\lambda_i = exp(\\beta_1 + \\beta_2 x_i)$$<br><br>\n", | |
| "\n", | |
| "これはつまり、各個体iにおける種子数の平均λ_iは、その個体の体サイズx_iによって上のように算出できる、という仮定である。<br>\n", | |
| "この両辺を対数化すると以下となる。(このlogがリンク関数)<br><br>\n", | |
| "$$log(\\lambda_i) = \\beta_1 + \\beta_2 x_i$$" | |
| ], | |
| "text/plain": [ | |
| "<IPython.core.display.Latex object>" | |
| ] | |
| }, | |
| "metadata": {}, | |
| "output_type": "display_data" | |
| } | |
| ], | |
| "source": [ | |
| "%%latex\n", | |
| "まず、「yデータは全体としてポアソン分布に従って分布している」という部分。これは本記事の最初に書いたように、以下の式で表される。<br><br>\n", | |
| "$$p(y|\\lambda) = \\frac{\\lambda^y exp(-\\lambda)}{y!}$$<br><br>\n", | |
| "\n", | |
| "ここで、「各x点においてそれぞれ特有の平均λのポアソン分布に従ってyが分布している。」という点を考慮して式を書き直す。つまり各x点でのy_iは、\n", | |
| "そのx点でのλ_iに従う、というように書き直す。<br><br>\n", | |
| "$$p(y_i|\\lambda_i) = \\frac{\\lambda_i^{y_i} exp(-\\lambda_i)}{y_i!}$$<br><br>\n", | |
| "\n", | |
| "そしてここで、このように考えてみる。yは種子数で、これはポアソン分布に従ってばらつき分布している。この分布の平均の違いは、体サイズxによって\n", | |
| "表すことができる、と。これを以下のような式で表してみる。<br><br>\n", | |
| "\n", | |
| "$$\\lambda_i = exp(\\beta_1 + \\beta_2 x_i)$$<br><br>\n", | |
| "\n", | |
| "これはつまり、各個体iにおける種子数の平均λ_iは、その個体の体サイズx_iによって上のように算出できる、という仮定である。<br>\n", | |
| "この両辺を対数化すると以下となる。(このlogがリンク関数)<br><br>\n", | |
| "$$log(\\lambda_i) = \\beta_1 + \\beta_2 x_i$$" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "ここでもう一度、何がやりたかったんだっけ、というのを再確認する。 \n", | |
| "今、種子数yのデータがあり、これをxによって回帰するモデルを作りたい。つまり、各x点において、最も尤もらしいyの値を一点プロットしていく。 \n", | |
| "ここで、yの分布は各xにおいてそれぞれの平均λをもったポアソン分布であると仮定している。 \n", | |
| "各x点におけるxとλの関係性をある特定の式(上の式)で固定できると仮定すると、上の式に各xを入れた時の全てのλの尤度の合計が最も大きくなるような\n", | |
| "パラメータβ1,β2を見つけることで、データ全体を通して最適なλリストを算出できる。 \n", | |
| "ポアソン分布において、期待値はλそのものなので、各x点におけるλが算出できれば、それはyの値を算出するのと同様である。" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 107, | |
| "metadata": { | |
| "collapsed": true | |
| }, | |
| "outputs": [], | |
| "source": [ | |
| "import statsmodels.api as sm\n", | |
| "# add constant to the data\n", | |
| "x_c = sm.add_constant(df.x)\n", | |
| "# modeling with GLM (log link function is automatically selected for Poisson)\n", | |
| "model = sm.GLM(df.y, x_c, family=sm.families.Poisson())\n", | |
| "result = model.fit()" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 108, | |
| "metadata": { | |
| "collapsed": false | |
| }, | |
| "outputs": [ | |
| { | |
| "data": { | |
| "text/html": [ | |
| "<table class=\"simpletable\">\n", | |
| "<caption>Generalized Linear Model Regression Results</caption>\n", | |
| "<tr>\n", | |
| " <th>Dep. Variable:</th> <td>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>Poisson</td> <th> Df Model: </th> <td> 1</td> \n", | |
| "</tr>\n", | |
| "<tr>\n", | |
| " <th>Link Function:</th> <td>log</td> <th> Scale: </th> <td>1.0</td> \n", | |
| "</tr>\n", | |
| "<tr>\n", | |
| " <th>Method:</th> <td>IRLS</td> <th> Log-Likelihood: </th> <td> -235.39</td>\n", | |
| "</tr>\n", | |
| "<tr>\n", | |
| " <th>Date:</th> <td>Tue, 19 Feb 2019</td> <th> Deviance: </th> <td> 84.993</td>\n", | |
| "</tr>\n", | |
| "<tr>\n", | |
| " <th>Time:</th> <td>22:03:36</td> <th> Pearson chi2: </th> <td> 83.8</td> \n", | |
| "</tr>\n", | |
| "<tr>\n", | |
| " <th>No. Iterations:</th> <td>7</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>[95.0% Conf. Int.]</th> \n", | |
| "</tr>\n", | |
| "<tr>\n", | |
| " <th>const</th> <td> 1.2917</td> <td> 0.364</td> <td> 3.552</td> <td> 0.000</td> <td> 0.579 2.005</td>\n", | |
| "</tr>\n", | |
| "<tr>\n", | |
| " <th>x</th> <td> 0.0757</td> <td> 0.036</td> <td> 2.125</td> <td> 0.034</td> <td> 0.006 0.145</td>\n", | |
| "</tr>\n", | |
| "</table>" | |
| ], | |
| "text/plain": [ | |
| "<class 'statsmodels.iolib.summary.Summary'>\n", | |
| "\"\"\"\n", | |
| " Generalized Linear Model Regression Results \n", | |
| "==============================================================================\n", | |
| "Dep. Variable: y 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: Tue, 19 Feb 2019 Deviance: 84.993\n", | |
| "Time: 22:03:36 Pearson chi2: 83.8\n", | |
| "No. Iterations: 7 \n", | |
| "==============================================================================\n", | |
| " coef std err z P>|z| [95.0% Conf. Int.]\n", | |
| "------------------------------------------------------------------------------\n", | |
| "const 1.2917 0.364 3.552 0.000 0.579 2.005\n", | |
| "x 0.0757 0.036 2.125 0.034 0.006 0.145\n", | |
| "==============================================================================\n", | |
| "\"\"\"" | |
| ] | |
| }, | |
| "execution_count": 108, | |
| "metadata": {}, | |
| "output_type": "execute_result" | |
| } | |
| ], | |
| "source": [ | |
| "# view the summary of the result\n", | |
| "result.summary()" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 119, | |
| "metadata": { | |
| "collapsed": false | |
| }, | |
| "outputs": [ | |
| { | |
| "data": { | |
| "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAEACAYAAACj0I2EAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAHdhJREFUeJzt3X2UXHWd5/H3N53u7coDzVN4DJPiYRgYnjMYGVd3ikEE\nEdBdnVll1EHH9biKcHQOB0c9ps+uswc4s7O7R+UPd2IExgSQUZ58CsxQetSDCRLoQMLDDHYTBNKV\nQMBAn3TT/Z0/blW6urqqq+rWvXWrbn1e59Sh69bvd+/33iq+uf3r3+9b5u6IiEj3W5R0ACIiEg0l\ndBGRlFBCFxFJCSV0EZGUUEIXEUkJJXQRkZSom9DNbJ2Z7TKzkYrtnzWzHWa2zcyujy9EERFpxOIG\n2qwHvgbcUtpgZjngMuAMd3/TzA6PJzwREWlU3Tt0d/858ErF5v8OXO/ubxbb7I4hNhERaULYMfST\ngf9kZg+Z2YNmdm6UQYmISPMaGXKp1e8Qdz/PzN4C3AGcEF1YIiLSrLAJfSfwPQB332JmM2Z2mLvv\nqWxoZioWIyISgrtbM+0bHXKx4qPkLuBPAczsZKC/WjIvC6prH2vXrk08hl6Nv5tjV/zJP7o9/jDq\n3qGb2QYgBxxmZs8Ba4FvAevNbBuwH/hoqKOLiEhk6iZ0d7+ixksfiTgWERFpgVaK1pHL5ZIOoSXd\nHH83xw6KP2ndHn8YFnaspuEDmHncxxARSRszw2P6o6iIiHQ4JXQRkZRQQhcRSQkldBGRlFBCFxFJ\nCSV0kRoKhQJbtmyhUCgkHYpIQ5TQRarYuPF2Vq06hQsv/BSrVp3Cxo23Jx2SSF2ahy5SoVAosGrV\nKUxMPAicCYyQyZzP2NiTrFixIunwpEdoHrpIBEZHRxkYyBIkc4Az6e9fxejoaHJBiTRACV2kQjab\nZXJyFCh9je4IU1NjZLPZ5IISaYASukiFFStWsG7dTWQy53PQQavJZM5n3bqbNNwiHU9j6CI1FAoF\nRkdHyWazSubSdmHG0JXQRUQ6kP4oKiLSw5TQRURSQgldRCQllNBFRFKibkI3s3VmtsvMRqq89tdm\nNmNmh8YTnoiINKqRO/T1wEWVG81sJXAhMBZ1UCIi0ry6Cd3dfw68UuWl/wNcG3lEIimnKo4Sl1Bj\n6GZ2ObDT3bdFHI9IqqmKo8SpoYVFZrYKuNfdzzSzDPAgcKG7/87MfgOc6+57avTVwiIRVMVRmhNm\nYdHiEMc5EcgCj5mZASuBX5vZGncfr9ZheHj4wM+5XI5cLhfisCLdrVTFcWJifhVHJXTJ5/Pk8/mW\n9tHoHXqW4A79jCqv/QZY7e7Vxtl1hy5SpDt0aUYsS//NbAPwS+BkM3vOzD5W0cSBpg4q0otUxVHi\npuJcIm2mKo7SCFVbFBFJCVVbFBHpYUroIiIpoYQuIpISSugiIimhhC4ikhJK6JJaKoKVfnqP51JC\nl1RSEaz003s8n+ahS+poiX369cJ7rHnoIswWwQr+R4fyIliSDnqPq1NCl9TJZrNMTo4CpW9NHGFq\naoxsNptcUBIpvcfVKaFL6qgIVvrpPa5OY+iSWiqClX5pfo9VnEtEJCX0R1ERkR6mhC4ikhJK6CIi\nKaGELiKSEkroIiIp0ciXRK8zs11mNlK27UYz22Fmj5rZP5nZQfGGKSIi9TRyh74euKhi2ybgNHc/\nG3gG+JuoAxMpp6p67RHFdY7zvdLnYGF1E7q7/xx4pWLbA+4+U3z6ELAyhthEAFXVa5cornOc75U+\nBw1w97oPYBUwUuO1e4ArFujrImGNj497JnOow2MO7vCYZzKH+vj4eNKhpUoU1znO96oXPwfF3NlQ\nji49Frfyj4GZfQmYcvcNC7UbHh4+8HMulyOXy7VyWOkhpap6ExPzq+qlbal3kqK4znG+V73wOcjn\n8+Tz+dZ20kjWp8odOnAl8AvgP9TpG/u/ZJJevXhnlgTdoXceQtyhNzpt0YqP4InZxcC1wOXuvr+1\nf1JEalNVvfaI4jrH+V7pc9CYusW5zGwDkAMOA3YBa4EvAgPAnmKzh9z90zX6e71jiNST5qp6nSSK\n6xzne9VLnwNVWxQRSQlVWxQR6WFK6CIiKaGELiKSEkroIiIpoYQuIpISSugiIimhhC49oVAosGnT\nJjZt2tR0pb5OqfAXRRw7duzg5ptvZseOHRFGFuiU69TTml1a2uwDLf2XhG3YcJsPDAw5nOSwxPv7\nl/mGDbc13DeTOdSHhlZ7JnNow/2iFkUcV111jUPG4WSHjF911dUdFZ/MRYil/0rokmrVaoDAIT44\neHDdOiCdUj8kiji2b99eTObl1yHj27dv74j4ZL4wCV1DLpJqo6OjLFp0HDBbpQ+y9PUdwejoaN2+\nAwPZOX1LFf7aKYo4Nm/eDFReh5XF7cnHJ9FQQpdUy2azzMzsBErfoDgCjDI9PU42m63bd3JydE7f\nqamxuv2iFkUca9asASqvw/PF7cnHJxFp9pa+2QcacpGEzY6hnxh6DP2gg87piDH0VuK46qqri8Mu\nvx/bGHrS1ylNCDHkouJc0hMKhQJbt24F4JxzzmmqUl+nVPiLIo4dO3awefNm1qxZw6mnntpx8cks\nVVsUEUkJVVsUEelhSugiIimhhC4ikhJK6CIiKaGELiKSEnUTupmtM7NdZjZStu0QM9tkZk+Z2U/M\nbCjeMEVEpJ660xbN7O3APuAWdz+zuO0GYI+732hm1wGHuPsXavTXtEWpqtq85Wbni7dr7nMUxynt\nY9myZezcuROofY4LHa/8NaCr535r7nptYaYtNrracxUwUvb8SeDI4s9HAU8u0DeORVTS5apV52u2\nKmK7KvxFcZzSPjKZM4qrNY+qeY4LHa/8tf7+5T4wMNS1FQ5VoXFhxFVtsUpCf7ni9ZcX6Bv7iUt3\nqVadb3DwYB8cPKThqojtqvAXxXGqV3w81OHBeee40PHmvjbucEjs5x8XVWisL0xCX9zy7wXFG/2F\nXhweHj7wcy6XI5fLRXRY6Ual6nwTE7PV+fr6jmBmZoD5VRFfZ3R0dN6v49X2UarwF+Wv7lEcp9o+\ngnukpVSe40LHA8pe2wIcT7UKh90wdNGu96+b5PN58vl8aztpJOsz/w59B3OHXHYs0LcN/5ZJN9Ed\nuu7QdYdeHzEOuWSBbWXPbwCuK/58HXD9An3bcOrSbapV52u2KmK7KvxFcZzZMfTTi2PoR9YdQ692\nvPLX+vuX+cDAUNdWOFSFxoWFSeiNzHLZAOSAw4BdwFrgLuC7BBXzx4A/d/e9Nfp7vWNIb9IsF81y\n0SyX2lRtUUQkJVRtUUSkhymhi4ikhBK6iEhKKKGLiKSEErqISBXT03DzzXDFFfDUU0lH0xgldBHp\neXv2wOWXg9nsY/FiuPJK2LgRnn026QgbE9XSf5EFNTu/vB1anQNd2b98nvm+ffvmba+cb1+tb/nz\n0vU67rjj2Ldv37z9Rn0+3SjMOT/yCFx2GbzwwsLt7rgD/uzPIgiynZpdidTsA60U7XkbNtzm/f3L\nHZY4nOQDA0OJrwpstdJfZf+rrrq6uBL0BIeMZzJnzNleWVWyWt/y57PX62iHjA8MnFrc7/FV4+3F\nyoX1znlqyv3WW71YWqD24+ij3bduTegkFkBcS/9beSih97bx8XEfHDy4o+qOtFpHpHptlozD94s1\nWiq3PzinZk31vg9WPB8qbqteA6Y83l6sizL/nB/3vr5v1k3e73mP+wsvJB19Y8IkdI2hS6xGR0fp\n6zuSysqAixatPFBBMImYBgayVKtUGLY/rAReJSh7VLl96YHnfX1HsGjRcQu2gWMJat4trbK/oEpj\nebytnk83cYetW2HlysOYmNjD7DmfxvT0f5vX/vvfn5vS77sPjj66rSG3lRK6xCqbzTI9vQv4DVD6\nFsMRZmaeP1CLJImYJidH58QzNTXWcDzV+sPzwBBQbfvrB55PT48zM7NzwTbwW+Cl4rbK/Y0Br8+J\nt9Xz6VQTE/D1r8/9Q+WiRbB6NUxOzk9d/f1/xfh4YU4Cf9/7Egg8Sc3e0jf7QEMuPS8YQ19WHBM+\nsaPG0MNW+qvsXxoHHxzMFse6T5+zvbKqZLW+5c9nr9eRc8bQBwezC46hd2PlwpkZ92eecV++vP54\nd+lx003dfc6NII5qi61ScS4BzXLRLJfA5CTccAN85SuN91m3Dj7+8eqvdcM5h6VqiyLSMX7xC3j7\n25vv87a3xRNPt1G1RRFpu1274AMfmDvWbVY/mT/xxPzBFCXz1mhhkYg05PXX4Vvfgquvbq7fMcfA\nc89BX188ccks3aGLyByTk/DP/xwk4vI77mXL6ifzSy+FmZm5d92//a2SebvoDl2kR01Pw44d8NWv\nwu23N9//Zz+Dd7wj+rgkPCV0kZRzD4Y8br4Z1q4Nt4+JCRgcjDYuiV5LQy5m9jkze9zMRszsO2Y2\nEFVgItK88XG49VY49ti5i3Gy2caS+Sc+MX/IxF3JvFuETuhmdgzwWWC1u59JcLf/wagCk3QpFAps\n2bKFQqHQkftLWqFQYNOmTWzatKmhc9q7F+6/Hy64YO4495FHwkc/Wr+SIMBDD81N2uPjBT75yS3s\n3l2oen0buebNvi+Nto/j2I3qqs9asyuRSg/gGIJ1yIcQJPN7gXdWaRffUirpClFXAkxbZcENG27z\ngYEhh5Mclnh//7ID5/T66+6bN7t/8pONr6Ks9piZqR9D6Zr29y/3gYGhBStEVrvmzb4vjbaP49iN\nSvKzRrurLQJXA78DdgG31mgT+4lL54q6EmDaKgvOns9I8XxecJgOnbjDlIGde03HvVplzKBiZu1r\n3uz70mj7RtrF9ZlI+rMWJqGH/qOomR0MvJeg/NurwJ1mdoW7b6hsOzw8fODnXC5HLpcLe1jpMqVK\ngBMT8ysBhlmqHfX+2m16Ovj2m1/+Ev72b+GZZ1YAe8paNFYK8KKL4Mc/jiamudd0C9UqY8J+qlVz\nLF3zZt+XRts30i6uz0S7P2v5fJ58Pt/aTpr9F6D0AD4A/P+y5x8Bvl6lXdz/kEkH69U79JkZ97Ex\n97vvdr/kktaGS/btizdW3aE3cl264w69lYS+BtgGDAIGfBv4TJV2bTh16WRRV8XrtCp7u3a5P/CA\n+6c/3Vrihi965Rh6u5Rf0/7+ZT4wMLRghciFxrEbfV8abR/HsRuV5GctTEJvqTiXma0lmNkyBWwF\nPuHuUxVtvJVjSDpEXRUviSp7e/cG9Ufyefjyl8PvZ8mSYBl9uU6oRll+TYG6FSLr7aORc2i0fRzH\nblRSFR1VbVEkAm+8Adu3w8MPw5e+BC+/HH5f+uhLWGESulaKSs+anISnn4bHHoO/+zt49NHw+5qZ\nCeaBiyRJCV1SrzSzZNu24MsSfvjD1va1SCXtpEMpoUtquMPOnUHivuMOuOWW8Pt6801VCJTuo4Qu\nXWl8HEZG4K674BvfCL8f3XFLmiihS0fbuze4477vPrjxxvD70Ri39AIldOkIb7wBjz8O99wTrKAM\n48MfDioNivQqJXRpq8nJ4EsV/vEfg5klYdx5J7z//eH6xjGnOMwc8oX6RBFjkvO2JUHNrkRq9oFW\nivakN990f/zx1qoE7tkTbUxxVM7bsOE27+9f7rDE4SQfGBhqaIVkrT5RxJhkdUKJDu2uttjQAZTQ\nU21mxv3hh90/+MFwSfvii+uXdo1CHHU5xsfHizVO5tc+WaiGSa0+27dvbznGJGufSLTCJHT9fV8a\nMjUVLHm/5JK5X6iwaBGcey7cdlvtvu96F7z2WvWU/qMfteePlaXKedUqBrayz76+I6lWnbDWfhfq\ns3nz5pZjbOQ847gW0hmU0GWO3bvh7rvhvPPmJu6BATj//CAB13LttcE0wMqk/ZOfwPLl7TuHarLZ\nLJOTo8BIccsIU1NjB+qWhN3n9PQu4Ddz9jsz83zN/S7UZ82aNS3H2Mh5xnEtpEM0e0vf7AMNuXSc\n/fvdt293//a33YeGmhsi+fzn3V97LekzCCeOynnBePiy4nj4iU2MoVfvE0WMSVYnlOjQ7mqLjVBx\nrmS4B4tvnnoquKu+/vrG+n3kI/CZz8CKFcEXC6dt0Y1muTTXRpKjaos9aP9++Nd/Db4B57774N57\nG6vwd911cM01cHRjX5AjIm2maosp5Q4vvRSUc7333iBxv/hi7fYXXACXXQZ//MeweDGceCIMDbUv\nXhFJhhJ6B5mYCOqT/OAHQdIu/kZe1emnw6WXBt8t+Yd/GAyRaGm7SG9TQm8zdxgdDca177134S/6\nPfzw4E770kvhrW+FY45R0haR2pTQY/Laa3D//bPj2nv21G777ncHSftd74Ljj1fZVhEJRwm9RePj\n8Pd/HyTuJ56o3W716iBpX3opnH029Pe3L0YR6Q0tJXQzGwL+ATgdmAE+7u6/iiKwbvG1r8ENNwQ/\nH310MERy2WXBIpylS5ONTUR6S0vTFs3s28BP3X29mS0Glrj7axVtNG2xy0T5TeztiqdWm/Ltu3fv\nZvPmzaxZs4bDDz+8I74hvtR+2bJl7Nu3b06/avsq3wbU7Bt33K3265b9JTlXP8y0xVZWgB4E/FsD\n7aJeQCUxarQKX7uq9bVSObB8e1/fUoeMw8kOGe/rWxJ57M1ek1L7TOYEh4xnMmcc6FdtX+Xb+vuX\n+8DAkGcyZxT7Hh9rdcYo+3XL/pKuSEk7qy0CZwG/AtYDjwDfBDJV2rXh1CUKjVbha1e1vlYqB86t\nXDg+r7ph8Hw8stibvSaz7R90mN8vqMg4u21w8OA65xPsK47qjFH265b9dUJFyjAJvZUx9MXAauAz\n7v6wmf1f4AvA2sqGw8PDB37O5XLkcrkWDitxKVXhm5iYX4Wv/NfNRtu1I55abUqVC4PtW6isbghZ\nYBR4SySxN3tNZtsvLcYyt/Ii7J+zra/vCCBT3FbtfFYBS5s+l7DvZdSfgU7bX7s+4+Xy+Tz5fL61\nnTT7L0DpARwJPFv2/O3AvVXaxf0PmUREd+jxxlq9ve7QO3F/3XqHHjqhB8fjp8DJxZ/XAjdUaRP7\niUt0Gq3C165qfa1UDizf3te3xIMx9N/30hh61LE3e01K7QcHsx6Mg58+b7y8fF/l2/r7lxXH0E93\nyPjgYDbW6oxR9uuW/SVdkTJMQm91lstZBNMW+4FngY+5+6sVbbyVY0j7aZZLvLFWa69ZLp25v26b\n5aJqiyIiHShMQk9ZtWsRkd6lhC4ikhJK6CIiKaGELiKSEkroIiIpoYQuIpISqofeol7+5vSF5lC3\no38U4p6n3OrrIk1pdiVSsw9SvFI06WpsSZqtFBiu2l+r/aMQdzW+Vl+X3ka7l/43dICUJvROqPWQ\nlGrn3kwtkVb7x3UOUdb6aPV1kTAJXWPoIZWqsZVXvCtVY0u7audeWe0vzv5RaPX9q9e/1ddFwlBC\nDymbzTI5OQqMFLeMMDU1dqDGRppVO3cYA15v6Bq02j8Krb5/9fq3+rpIKM3e0jf7IKVDLu7JV2NL\n0uwYeLhqf632j0Lc1fhafV16G+2uttiItBfn6uVZCprlolkuEh9VWxQRSQlVWxQR6WFK6CIiKaGE\nLiKSEkroIiIpoYQuIpISLSd0M1tkZo+Y2T1RBCQiIuFEcYd+DbA9gv1IBykUCmzZsoVCoZB0KG1T\nfs71zr+brk83xSotanYlks9dBboSuB/IAffUaBPbSiqJRy9WASw/54GBIe/vX5aKKondFKvMRbur\nLQLfBc4G/kQJPR16sQpg9eqPhziMd3WVxG6KVeYLk9BDf8GFmb0H2OXuj5pZDqi5oml4ePjAz7lc\njlwuF/awErNSFcCJiflVANO6NL3aOUMWGAXeMuf8u+n6dFOsAvl8nnw+39I+Qi/9N7P/BXwYeBPI\nAMuB77n7RyvaedhjSPsVCgVWrTqFiYkHCRLbCJnM+YyNPZnaJFDtnINRxKeAF+ecfzddn26KVeYL\ns/Q/qoqKGnJJkV6sAlh+zqUx9DRUSeymWGUukqq2aGZ/Avy1u19e5TWP4hjSXr1YBbD8nIHUVEns\nplhllqotioikhKotioj0MCV0EZGUUEIXEUkJJXQRkZRQQhcRSQkldBGRlFBCTzlV2pula1Gbrk06\nKKGn2MaNt7Nq1SlceOGnWLXqFDZuvD3pkBKja1Gbrk16aGFRSqmOxyxdi9p0bTqXFhbJAaVKe8H/\npFBeaa/X6FrUpmuTLkroKZXNZpmcHCWoHAgwwtTU2IE6Jb1E16I2XZt0UUJPqRUrVrBu3U1kMudz\n0EGryWTOZ926m3ry12hdi9p0bdJFY+gpp0p7s3QtatO16TyqtigikhL6o6iISA9TQhcRSQkldBGR\nlFBCFxFJidAJ3cxWmtm/mNkTZrbNzK6OMjAREWlO6FkuZnYUcJS7P2pmy4BfA+919ycr2mmWi4hI\nk9o6y8XdX3L3R4s/7wN2AMeG3Z9I3FRRUNIukjF0M8sCZwO/imJ/IlFTRUHpBS0vLCoOt+SB/+nu\nd1d5XUMukihVFJRuFGbIZXGLB1wM3AncWi2ZlwwPDx/4OZfLkcvlWjmsSFNKFQUnJuZXFFRCl06R\nz+fJ5/Mt7aOlO3QzuwXY7e6fX6CN7tAlUbpDl27U1j+Kmtl/BP4C+FMz22pmj5jZxWH3JxIXVRSU\nXqHiXNIzVFFQuomqLYqIpISqLYqI9DAldBGRlFBCFxFJCSV0EZGUUEIXEUkJJXQRkZRQQhcRSQkl\ndBGRlFBCFxFJCSV0EZGUUEIXEUkJJXQRkZRQQhcRSQkldBGRlFBCFxFJCSV0EZGUUEIXEUkJJXQR\nkZRoKaGb2cVm9qSZPW1m10UVlIiINC90QjezRcDXgYuA04APmdkpUQXWKfL5fNIhtKSb4+/m2EHx\nJ63b4w+jlTv0NcAz7j7m7lPAbcB7owmrc3T7h6Kb4+/m2EHxJ63b4w+jlYR+LLCz7PnzxW0iIpIA\n/VFURCQlzN3DdTQ7Dxh294uLz78AuLvfUNEu3AFERHqcu1sz7VtJ6H3AU8AFwIvAZuBD7r4j1A5F\nRKQli8N2dPdpM7sK2EQwdLNOyVxEJDmh79BFRKSzxPpHUTMbMrPvmtkOM3vCzN4a5/GiYmYnm9lW\nM3uk+N9XzezqpONqhpl9zsweN7MRM/uOmQ0kHVMzzOwaM9tWfHT8tTezdWa2y8xGyrYdYmabzOwp\nM/uJmQ0lGeNCasT/geJnaNrMVicZXz014r+xmHseNbN/MrODkoxxITXi/x9m9lgxB/3YzI6qt5+4\nZ7n8P+CH7n4qcBbQFUMy7v60u5/j7quBPwJeB76fcFgNM7NjgM8Cq939TIKhtQ8mG1XjzOw04K+A\nc4GzgUvN7IRko6prPcEiu3JfAB5w9z8A/gX4m7ZH1bhq8W8D/jPw0/aH07Rq8W8CTnP3s4Fn6L7r\nf6O7n+Xu5wA/ANbW20lsCb34r+E73H09gLu/6e6vxXW8GL0T+Dd331m3ZWfpA5aa2WJgCfBCwvE0\n41TgV+6+392ngZ8B/yXhmBbk7j8HXqnY/F7g5uLPNwPva2tQTagWv7s/5e7PAE3NtEhCjfgfcPeZ\n4tOHgJVtD6xBNeLfV/Z0KTBDHXHeoR8P7Daz9cWhi2+aWSbG48XlvwIbkw6iGe7+AvC/geeA3wJ7\n3f2BZKNqyuPAO4pDFkuAS4DjEo4pjCPcfReAu78EHJFwPL3s48CPkg6iWWb2VTN7DrgC+Eq99nEm\n9MXAauAbxaGLNwh+Be0aZtYPXA58N+lYmmFmBxPcHa4CjgGWmdkVyUbVOHd/ErgBuB/4IbAVmE40\nqGhoBkICzOxLwJS7b0g6lma5+5fd/feA7xAMoy4ozoT+PLDT3R8uPr+TIMF3k3cDv3b3QtKBNOmd\nwLPu/nJxyOJ7wNsSjqkp7r7e3c919xywF3g64ZDC2GVmRwIU/6A1nnA8PcfMriT4Da9rbmhq2AC8\nv16j2BJ68VfNnWZ2cnHTBcD2uI4Xkw/RZcMtRc8B55nZoJkZwbXvij9Il5jZiuJ/f4/gD3PdcHdl\nzB1vvge4svjzXwJ3tzugJlXGX/lap5sTv5ldDFwLXO7u+xOLqnGV8Z9U9tr7aOD/4VjnoZvZWcA/\nAP3As8DH3P3V2A4YoeLY7Rhwgrv/Lul4mmVmawlmtkwRDFl8olgVsyuY2c+AQwni/5y755ONaGFm\ntgHIAYcBuwhmJNxFMFx3HMFn6c/dfW9SMS6kRvyvAF8DDif4LelRd393UjEupEb8XwQGgD3FZg+5\n+6cTCbCOGvG/B/gDguHGMeBT7v7igvvRwiIRkXRQtUURkZRQQhcRSQkldBGRlFBCFxFJCSV0EZGU\nUEIXEUkJJXQRkZRQQhcRSYl/Bx0FevqNj/+9AAAAAElFTkSuQmCC\n", | |
| "text/plain": [ | |
| "<matplotlib.figure.Figure at 0x1122a6ef0>" | |
| ] | |
| }, | |
| "metadata": {}, | |
| "output_type": "display_data" | |
| } | |
| ], | |
| "source": [ | |
| "# plot result\n", | |
| "import math\n", | |
| "# get params\n", | |
| "b1,b2 = result.params\n", | |
| "# plot\n", | |
| "fig, ax = plt.subplots(1,1)\n", | |
| "raw = ax.scatter(df.x, df.y) # plot original data\n", | |
| "fitted = ax.plot(df.x, df.x.map(lambda x:math.exp(b1+b2*x))) # plot the fitted line" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "今度は、施肥処理の有無を考慮したモデルを考えてみる。" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 122, | |
| "metadata": { | |
| "collapsed": false | |
| }, | |
| "outputs": [], | |
| "source": [ | |
| "# binarize f \n", | |
| "df['f_'] = df['f'].map(lambda x: 0 if x=='C' else 1) " | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 127, | |
| "metadata": { | |
| "collapsed": true | |
| }, | |
| "outputs": [], | |
| "source": [ | |
| "# add constant to the data\n", | |
| "x_f_c = sm.add_constant(df[['x','f_']])\n", | |
| "# modeling with GLM (log link function is automatically selected for Poisson)\n", | |
| "model = sm.GLM(df.y, x_f_c, family=sm.families.Poisson())\n", | |
| "result = model.fit()" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 128, | |
| "metadata": { | |
| "collapsed": false | |
| }, | |
| "outputs": [ | |
| { | |
| "data": { | |
| "text/html": [ | |
| "<table class=\"simpletable\">\n", | |
| "<caption>Generalized Linear Model Regression Results</caption>\n", | |
| "<tr>\n", | |
| " <th>Dep. Variable:</th> <td>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> 97</td> \n", | |
| "</tr>\n", | |
| "<tr>\n", | |
| " <th>Model Family:</th> <td>Poisson</td> <th> Df Model: </th> <td> 2</td> \n", | |
| "</tr>\n", | |
| "<tr>\n", | |
| " <th>Link Function:</th> <td>log</td> <th> Scale: </th> <td>1.0</td> \n", | |
| "</tr>\n", | |
| "<tr>\n", | |
| " <th>Method:</th> <td>IRLS</td> <th> Log-Likelihood: </th> <td> -235.29</td>\n", | |
| "</tr>\n", | |
| "<tr>\n", | |
| " <th>Date:</th> <td>Tue, 19 Feb 2019</td> <th> Deviance: </th> <td> 84.808</td>\n", | |
| "</tr>\n", | |
| "<tr>\n", | |
| " <th>Time:</th> <td>22:19:38</td> <th> Pearson chi2: </th> <td> 83.8</td> \n", | |
| "</tr>\n", | |
| "<tr>\n", | |
| " <th>No. Iterations:</th> <td>7</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>[95.0% Conf. Int.]</th> \n", | |
| "</tr>\n", | |
| "<tr>\n", | |
| " <th>const</th> <td> 1.2631</td> <td> 0.370</td> <td> 3.417</td> <td> 0.001</td> <td> 0.539 1.988</td>\n", | |
| "</tr>\n", | |
| "<tr>\n", | |
| " <th>x</th> <td> 0.0801</td> <td> 0.037</td> <td> 2.162</td> <td> 0.031</td> <td> 0.007 0.153</td>\n", | |
| "</tr>\n", | |
| "<tr>\n", | |
| " <th>f_</th> <td> -0.0320</td> <td> 0.074</td> <td> -0.430</td> <td> 0.667</td> <td> -0.178 0.114</td>\n", | |
| "</tr>\n", | |
| "</table>" | |
| ], | |
| "text/plain": [ | |
| "<class 'statsmodels.iolib.summary.Summary'>\n", | |
| "\"\"\"\n", | |
| " Generalized Linear Model Regression Results \n", | |
| "==============================================================================\n", | |
| "Dep. Variable: y 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: Tue, 19 Feb 2019 Deviance: 84.808\n", | |
| "Time: 22:19:38 Pearson chi2: 83.8\n", | |
| "No. Iterations: 7 \n", | |
| "==============================================================================\n", | |
| " coef std err z P>|z| [95.0% Conf. Int.]\n", | |
| "------------------------------------------------------------------------------\n", | |
| "const 1.2631 0.370 3.417 0.001 0.539 1.988\n", | |
| "x 0.0801 0.037 2.162 0.031 0.007 0.153\n", | |
| "f_ -0.0320 0.074 -0.430 0.667 -0.178 0.114\n", | |
| "==============================================================================\n", | |
| "\"\"\"" | |
| ] | |
| }, | |
| "execution_count": 128, | |
| "metadata": {}, | |
| "output_type": "execute_result" | |
| } | |
| ], | |
| "source": [ | |
| "# view the summary of the result\n", | |
| "result.summary()" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 131, | |
| "metadata": { | |
| "collapsed": false | |
| }, | |
| "outputs": [], | |
| "source": [ | |
| "df_c = df.query('f_==0')\n", | |
| "df_t = df.query('f_==1')" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 137, | |
| "metadata": { | |
| "collapsed": false | |
| }, | |
| "outputs": [ | |
| { | |
| "data": { | |
| "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAEACAYAAACj0I2EAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XuUXGWZ7/Hvk6SbriSkIRrCCsQuxGHCECKJGF0ePVMI\nKMoR1DMzjpzjDf0D5bZwlgOoY3o8jgrKqGuUuUgbYZxOIhhuDkLESeGAQEfITdIgErqJoOnmTqBJ\nd9LP+WNXdao7VV171712/T5r1UrXrne/+9l7V57e/da7nzJ3R0REmt+MegcgIiKVoYQuIhITSugi\nIjGhhC4iEhNK6CIiMaGELiISE0UTupn1mNluM9s2ZfmFZtZvZtvN7OvVC1FERMKYFaLNauCfgOuy\nC8wsBbwPONHd95nZa6sTnoiIhFX0Ct3d7waem7L408DX3X1fps3TVYhNREQiKHUM/Tjgf5rZfWa2\n0cxOrmRQIiISXZghl0LrHe7ubzWzNwM/Bl5fubBERCSqUhP6LmA9gLtvMrNxM3uNuz8ztaGZqViM\niEgJ3N2itA875GKZR9ZNwDsBzOw4oC1fMs8Jqmkfq1atqnsMrRp/M8eu+Ov/aPb4S1H0Ct3MeoEU\n8BozewJYBfwAWG1m24G9wEdL2rqIiFRM0YTu7ucUeOkjFY5FRETKoDtFi0ilUvUOoSzNHH8zxw6K\nv96aPf5SWKljNaE3YObV3oaISNyYGV6lD0VFRKTBKaGLiMSEErqISEwooYuIxIQSuohITCihixQw\nPDzMpk2bGB4erncoIqEooYvksW7NGpZ0dXHe6aezpKuLdWvW1DskkaI0D11kiuHhYZZ0dbFxZIRl\nwDbglESChwcHWbBgQb3DkxaheegiFTAwMECyvZ1lmefLgK62NgYGBuoYlUhxSugiUySTSQZGR8l+\nie42YHBsjGQyWceoRIpTQheZYsGCBVzd08MpiQQr5s3jlESCq3t6NNwiDU9j6CIFDA8PB8MvyaSS\nudRcKWPoSugiIg1IH4qKiLQwJXQRkZhQQhcRiQkldBGRmCia0M2sx8x2m9m2PK/9jZmNm9n86oQn\nIiJhhblCXw28e+pCMzsaOB0YrHRQIiISXdGE7u53A8/leelbwOcqHpFIzKmKo1RLSWPoZnYWsMvd\nt1c4HpFYUxVHqaZQNxaZWRdwq7svM7MEsBE43d1fMrPHgZPd/ZkC6+rGIhFUxVGiKeXGolklbOdY\nIAlsNTMDjgYeMLOV7j6Ub4Xu7u6Jn1OpFKlUqoTNijS3iSqOIyPA5CqOSuiSTqdJp9Nl9RH2Cj1J\ncIV+Yp7XHgdWuHu+cXZdoYtk6ApdoqjKrf9m1gv8CjjOzJ4ws09MaeJApI2KtCJVcZRqU3EukRpT\nFUcJQ9UWRURiQtUWRURamBK6iEhMKKGLiMSEErqISEwooYuIxIQSusSWimDFn87xZEroEksqghV/\nOscH0zx0iR3dYh9/rXCONQ9dhJwiWJnnuUWwJB50jvNTQpfYSSaTDIyOkv3OxG3A4NgYyWSyjlFJ\nJekc56eELrGjIljxp3Ocn8bQJbZUBCv+4nyOVZxLRCQm9KGoiEgLU0IXEYkJJXQRkZhQQhcRiQkl\ndBGRmAjzJdE9ZrbbzLblLLvSzPrNbIuZ/cTM5lU3TBERKSbMFfpq4N1Tlm0ATnD3k4BHgcsrHZhI\nLlXVq41KHOdqniu9D6ZXNKG7+93Ac1OW3enu45mn9wFHVyE2EQDWrFlHV9cSTj/9PLq6lrBmzbp6\nhxRLlaheWM1zpfdBCO5e9AF0AdsKvHYLcM4067pIqYaGhjyRmO+w1cEdtnoiMd+HhobqHVqsDA0N\n+fxEwrcGB9m3gs9PJCId52qeq1Z8H2RyZ6gcnX3MKueXgZl9ARhz997p2nV3d0/8nEqlSKVS5WxW\nWsjAwADt7UlGRg7U1Wtr62JgYCB2t3rX00T1wpERYHL1wrDHuZrnqhXeB+l0mnQ6XV4nYbI+ea7Q\ngY8D9wCHFFm36r/JJL5a8cqsHnSF3ngo4Qo97LRFyzyCJ2ZnAJ8DznL3veX9ShEpbMGCBfT0XE0i\ncQrz5q0gkTiFnp6rY3NV1igqUb2wmudK74NwihbnMrNeIAW8BtgNrAI+D7QDz2Sa3efunymwvhfb\nhkgxca6q10gqcZyrea5a6X2gaosiIjGhaosiIi1MCV1EJCaU0EVEYkIJXUQkJpTQRURiQgldRCQm\nlNClJQwPD7NhwwY2bNgQuVJfo1T4q0Qc/f39XHvttfT391cwskCjHKeWFvXW0qgPdOu/1Flv71pv\nb+90eIPDbG9rm+u9vWtDrbu2t9fnJxK+orPT5ycSvra3t8rR5tfbu9YTifne2bnCE4n5oePPddEF\nF3oC/A3M8AT4hRdcWLH4GuU4xQkl3PqvhC6xlq8GCBzuHR2HFa0DUon6JpVQiTomO3bs8ERmH7L7\nkgDfsWNHReJrhOMUN6UkdA25SKwNDAwwY8ZigvqBZP5NMnPmEQwMDBRdN9nePmnNbAXCWspWGszd\nh2ylwbD6+vo4ihmT9mURM+jr66tIfI1wnERj6BJzyWSS8fFdQPYbFLcBA+zfP0QymSy67sDo6KQ1\nB8fGiq5XaclkktHRAXL3YWxsMFIcK1eu5EnGJ+3LU4yzcuXKisTXCMdJ0JCLxN+BMfRjSx5DXz5v\nXkOMoc+bt7zkMfQLM2Pox1ZxDL3exylOKGHIRcW5pCUMDw+zefNmAJYvXx6pUl+jVPirRBz9/f30\n9fWxcuVKjj/++IaLTw5QtUURkZhQtUURkRamhC4iEhNK6CIiMaGELiISE0roIiIxUTShm1mPme02\ns205yw43sw1m9oiZ3WFmndUNU0REiik6bdHM3g7sAa5z92WZZVcAz7j7lWZ2KXC4u19WYH1NW5S8\n8s1bjjpfvFZznyuxnWwfc+fOZdeuXUDhfZxue7mvAU0991tz1wsrZdpi2Ls9u4BtOc8fBhZmfj4S\neHiadatxE5U0uXzVA6NWRaxVhb9KVDrM9pFInOiQcDiy4D5Ot73c1w5pm+Od7e1NW+FQFRqnR7Wq\nLeZJ6M9Oef3Zadat+o5Lc8lXPbCj4zDv6Dg8dFXEWlX4q0Slw/wVH+c7bDxoH6fb3uTXhryDjqat\ncKgKjcWVktBnlf13QeZCf7oXu7u7J35OpVKkUqkKbVaaUbZ64MjIgfp8M2cewfh4OwdXRXyZgYGB\ng/4cn6jwNzIy0Tpb4a+Sf7rnizVb6TDsdvL1EVwjzWHqPk63PSDntU0cxSEs49UDPVZh/6ulVuev\nmaTTadLpdHmdhMn6HHyF3s/kIZf+adatwe8yaSa6QtcVuq7Qi6OKQy5JYHvO8yuASzM/Xwp8fZp1\na7Dr0mzyVQ+MWhWxVhX+KlHp8MAY+tLMGPrComPo+baX+1p7Zgy9WSscqkLj9EpJ6GFmufQCKeA1\nwG5gFXATcD2wGBgE/srdny+wvhfbhrQmzXLRLBfNcilM1RZFRGJC1RZFRFqYErqISEwooYuIxIQS\nuohITCihi4jEhBK6iMhUe/fCpz8NZsHjZz+rd0ShVOrWf5FpRZ1fXgvlzoGeun7uPPM9e/YctHzq\nfPt86+Y+zx6vxYsXs2fPnoP6rfT+NJM9e+CGG+Caa0a55572ieX798OMqJepIyNw8cXw/e/nf33u\nXDjmmNKDraWodyJFfaA7RVteb+9aP6RtjicwPxa8s7297ncFllvpb+r6F15woScS831OIukJ8BMT\niUnLcysnFlo32+aCCy7ytrZDHWa7sdAT4Me3t3sCfHYimfdu1UpUhGxUv/2t+6c+5ZlSCAc/ZtrT\n3jHrW37VN35avLPnn3c/99zCnWUft95a/R0rgmrd+l/OQwm9tQ0NDQV1Whqo7ki5dUTyrZ8Ahxu9\ng3zLN06qWZN/3QNtgtIAnQ4bD+qvg4TDxkn1ZCpRb6YRjI253367+7veNX2uPeYY90ceGS5+Dnfv\ndv/oR4sn71tucR8fr9+OF1BKQtcYulTVwMAAM2cuzFQGDCwDjp4xY6KCYD1iSra3T4onW+mv1PUX\nMQN4gaPIt3zOxJKZM49g8YwZ07aBowhq3s3J018bMGdSBcZshcbcSpW5rzeip5+Gq68ORjKyw9Rt\nbXDGGbBhw+S211wDL798IAPv3AkvvPD4QefgprExFhxxxIEOFy6E666b3NnNNwfjMrkp/X3vC9rH\ngBK6VFUymWT//t08yV6y32G4Dfj9+PhELZJ6xDQwOjopnsGxsdDx5Fv/KcaBTp4k3/KXJ5bs3z/E\nrvHxadvAk8AfgZfz9DcGvMzY2OBEvMlkktHRgUyLoGXu6/U0Pg5btsBnP3sgz5rBggVw/vkw9XfO\n/Pnw6KOT8+0nPwmzZ2cauEN/Pyedey4PvPDCpHXfsW/f5M7+5E9g27bJnZ11VgmD7E0k6iV91Aca\ncml5vb1rvb1Bx9BLrfQ3df3sOHiio8sT4EunjKHnVk4stG62TTCGPtdhtsMRngBfkhlDT3R0TTuG\nXk5FyHLt2eP+k5+4n3VW8VEOcD/8cPfNm6fpcN8+91/9yj2RKNrZPTNm+M3f/W7N9rUWqEa1xXKp\nOJeAZrnEaZaLOwwOwvXXw7XXwkMPFV+nsxM2boTlyws0GB2F226DD3wgXBDnnw/f+Q7Dzz4b25k9\nqrYoIhW1dy9s2gQ/+lGQvF99Ndx63/kOXHRRgRdfegl6euCSS8J1duqpQbJvby/eNkZKSeiahy4i\nAAwNBR9I/vCH8ItfhF/ve987cA/OQXbtgq99Df75n8N1dtll8NWvxuZDylpTQhdpMfv2wcMPBzfm\nXHvtwR9MTueSS+Cb38zzuaI73Hd/8OnnvfeG73BgALq6wreXaSmhi8TYCy/AffcFQyY/+lG0df/+\n7+Hyy4PphJO8+irccgt86EPhO0sm4Te/gTlzijaV0imhi8TA+Dg8/jjceWdw1R3lIhng7LOD9To7\np7zw1FOwejV88YvhO/vYx+Bf/gU6OqIFIWVTQhdpMq+8EkyvXr8+SMJDQ9HWv/hiuOoqmDkzZ+G+\nfUGn3/gGrF0brcORESXvBlFWQjezS4BPAuPAduAT7j5aicBEWp17cIF8333Q2xsk8KhWrQoekz5j\nfP55+Old8LnPBXfxhHXRRfCtb8X7xpwmV/K0RTNbBNwNLHH3UTNbB/ynu183pZ2mLUrF50jHrbLg\nk08Oc8stv2PLlvncffex7NgR/Vrry18ORkYmknd2HGb9evjbv43U154vfpH+s86adu58PlHPS9j2\n1dh2WPV6r5UybbGcO0AXAYPA4QRX+rcCp+VpV42bqKSJVLoSYLmVEuttaMh9wwb3z3zG/ZBDwt1V\nOfXx4x9P6XTPHvd773X/yEeid3bJJZPO0SFtc7yzvX3S8Q1zDqOel7Dtw7SrVrXJer7XqHW1ReAi\n4CVgN/DvBdpUfcelcVW6EmC5lRJrad8+9x073K+91v3d7y4tcYP7unU5nY6Pu+/a5X7DDe5Ll0bv\n7PvfPyjOyedoKG9lzI6Ow6Y9h1HPS9j2YdpVq9pkvd9rpST0ksfQzeww4GygC3gBuMHMznH33qlt\nu7u7J35OpVKkUqlSNytNJlsJcGTk4EqApfz5OlHpcGQk09uBSon1HHp58cXgM8V0Opjg8eSTpfXz\nDPN5M/vZPedI0htWc/Ls2XDHHWCXRe/siSdg8eKizSafo02ZypjBLaHZypiPsZB81RyzxzzqeQnb\nPky7Sr/HosZYKel0mnQ6XV4nUX8DZB/AXwDfz3n+EeC7edpV+xeZNLC4XaGPj7vv3Om+fr37xz5W\n+lV3R0fQV3Z/fpN54dlSOmtvd9+/v+R90hV64ePSbFfo5ST0lQQzWzoAA34InJ+nXQ12XRpZpSsB\nllspMaxXXnHv63P/3vfcTzih9OS9dWtOp/v2uff3u/f0uC9eHL2zvr6q7GvuOWrPjKHnHt8w5zDq\neQnbPky7alWbrNV7LZ9SEnpZxbnMbBXw18AYsBn4lLuPTWnj5WxD4qGRZ7m4wx/+ENTt/sUv4B//\nsbR+3vMe+OlPc2b1Zcdhbr0Vrrwyeoc1nt+de0yBohUii/WhWS7lUbVFkSJGR4M6Jlu2wL/+K/zq\nV6X1MzgIr3td5ol7UJPkgQfg29+Ge+6J1tlXvxrcYy+SQwldJMczz8DWrXDXXcEc7VL82Z8FJUgm\n5naPjAQL7rgD/u7vone4d2/LlYGV0qh8rrSk/fvhd7+DBx8MboW/447S+plU+M8d/vjHoNNTrwq+\nnSGK224LxmBEakgJXZrKSy8Fw9J33glXXBFcMEe1dCls356zYGwsGIe5/o7gdvgourqCuzFVv1sa\ngBK6NCT3YJy6ry8o+3rrrdH7WLIk+Hq0SaVHnn02GOte+YXgq3iiUBEqaXAaQ5e6yw5L33wz/MM/\nlNbHiy/CoYfmLBgfD8Zh1q2DL30pWmf60gVpABpDl4aWHZbOfkjZ3x+9j4ceCj6onGTPHvjv/4b3\nvjdaZ0rcEjNK6FIVY2PB1MCvfQ1uvDH6+j//OZx22pSF7rBzJ7z3QvjZz8J39utfw5veBFRnTvHw\n8DCbN28GYPny5aHnXxdapxIx1nPettRR1DuRoj7QnaKx99RT7p//fGl3UT74YIFOX3nF/YoronV2\n//3TxlmNinxre3v90LY2nw3+BvDO9vZQd0gWWqcS1f3C9NHsFStbAbW89T/0BpTQY+Oll9x/8AP3\nRYui5dnPftb95ZcLdDo+7n7nndE6vPfeyLFXo97H0NCQH9bR4Ydn6nyErWFSaJ0dO3aUXTskbO2T\nZqlY2cpKSej66hE5yIsvHhjyMDvwOPRQOPfc4Ft08lm2LBgjn5qBr7oKZic8+Hact799cqczZuQZ\nWwHe9ragZGG+lP7Wt0bep2xFvnwVA0s1MDDAwpkzOWZSr0F1wkL9TrdOX19fUN0vZ3m2ul+UmIr1\nEaaNNCcl9Bb27LNw991w3nmTc2xnJ7zrXUFdk3wuvTT4HHJqnt26FRYuJPiw8ZvfhLe8ZXLiPu64\ng2+Lf/3rg09Hp3Z2zz2waFHF9jWZTDI6OgBsyyzZxtjY4ETdklL73L1/P49P6hV+Pz5esN/p1lm5\nciUDo6OTlg+OjUWKMZlMFu0jTBtpUlEv6aM+0JBLXY2Pu+/e7b5xo/vll4cb0Vi2zP2d73S/8EL3\nF18ssoE1a9xPPrlwZwsXup9/vns67T46WotdLqgaFfnW9vb63Mx4+LERxtALrVOJ6n5h+qhnFUEJ\nh1pXWwxD89Brwz0YCtm+HW6/PZhZ8sQT06+zcCF85Stwzjkwe3aJG87eIXnkkfCXfxk83va2KV8p\n3zg0yyVaG6kfFedqAePjQaLu64Obbgq+/3fv3vxtly4NbtgB6O6Giy+Gww6rWagiUgbdWBQj+/fD\nY48FH07eeGPh8WyAU0+F978/uDA+7jiYO7d2cYpI41BCr7PR0QO3va9ff+CKeqpDDoEPfCB4rFwZ\n1OKeoY+0RSSHEnqNvPoq/PKXQdJevx6Gh/O3W7wYPvhBOPtsWLEimHEiIhKGEnqFPfdc8DVkN944\n/S3vJ58cJO73vS+oCjhLZ0JEyqQ0Uqbbb5/+ewzOPDMYJjnzzGAiiIhItZSV0M2sE7gGWAqMA+e6\n+/2VCKxZvPJKMHMkO7592mmQSNQ7KhFpRWVNWzSzHwJ3uftqM5sFzHb3F6e00bTFJlPJb2KvVTyF\n2uQuf/rpp+nr62PlypW89rWvbYhviM+2nzt3Lnv27Jm0Xr6+cpcBBdetdtzlrtcs/dVzrn4p0xbL\nuQN0HvBYiHaVvoFKqihsRcJaVesLs51CMecub5uZ8AT4G5jhCfBZMxMVrboYNtZ8cc9JJD0BfmIi\nMbFevn3KXXZI2xzvbG/3ExPBfh2Ts25UpVahrPR7oNLVMMuNr94VKalltUXgjcD9wGrgQeDfgESe\ndjXYdamEsBUJa1WtL2zlwHwx79ixI2f5kHfQMamfDjochipSdTFsrFPbB/Ft9A4OXq+j47BJ+9TR\ncdi0+zMffGOJ1RlLqUJZ6fdApathlhtfI1SkLCWhlzOGPgtYAZzv7r82s28DlwGrpjbs7u6e+DmV\nSpFKpcrYrFRLtiLhyMjBFQlz/9ycqNaX+Ybm3Gp9lfyzNMx2CsXc19eXs3wTR3EIy3h1op9FtLOT\nAeDNefexGrFObR/EN4ejaGcZB9Y7esYMHmMhufUYZ848Akhklh28P13AnCLbLBR3mHNe7v5WK45q\nxVer93iudDpNOp0ur5OovwGyD2AhsDPn+duBW/O0q/YvMqkQXaFXN9ap7XWFXn4c1YqvWa/QS07o\nwfa4Czgu8/Mq4Io8baq+41I5YSsS1qpaX5jtFIo5d/mszBj6sTlj6JWsuhg21nxxJzq6PAG+NM8Y\nem6MucvaM2PoSzNj6MmOjrLH0KMej0q/BypdDbPc+OpdkbKUhF7uLJc3EkxbbAN2Ap9w9xemtPFy\ntiG1p1ku1Y01X3vNcmnM/pptlouqLYqINKBSErrKO4mIxIQSuohITCihi4jEhBK6iEhMKKGLiMSE\nErqISEyoHnqZWvmb06ebQ12L9Suh2vOUy31dJJKodyJFfRDjO0UrXR2umUzc5Zg40SHhicQxJd1l\nmK0WODuRrPkxrHY1vnJfl9ZGrW/9D7WBmCb0SteeaCb59h2CuiSl1gHpIBF6/UrtQzVrfZT7ukgp\nCV1j6CXKVofLrYiXrQ4Xd/n2PVvvL8wxmKhkl7P2ItpCr18J+WLIVtOrxPrlvi5SCiX0EiWTSUZH\nB4BtmSXbGBsbnKixEWf59h0GgZdDHYNkMsnA6OiktZ9iLPT6lZAvhsGxsdDbLrZ+ua+LlCTqJX3U\nBzEdcnGvfHW4ZnJgDH2pQ8I7OqKNgWfHj7PVAhMdXXUbQ69WNb5yX5fWRq2rLYYR9+JcrTxLQbNc\nNMtFqkfVFkVEYkLVFkVEWpgSuohITCihi4jEhBK6iEhMKKGLiMRE2QndzGaY2YNmdkslAhIRkdJU\n4gr9YmBHBfqRBjI8PMymTZsYHh6udyg1k7vPxfa/mY5PM8UqZYp6J5JPvgv0aODnQAq4pUCbqt1J\nJdXRilUkcysfdra3e3vbnIL730zHRxUdmxe1rrYIXA+cBPy5Eno8tGIVyfzVHzschg7a/2Y6Pqro\n2NxKSeglf8GFmZ0J7Hb3LWaWAgre0dTd3T3xcyqVIpVKlbpZqbJsJcWRkYOrSMb11vSJyocjI0C2\n+mM7OxkA3jxp/5vp+OTbr2xFx0aLVSCdTpNOp8vqo+Rb/83sq8D/BfYBCeBQYL27f3RKOy91G1J7\nw8PDdHUtYWRkI0EK2EYicQqDgw/HNgkMDw+zpKuLjSMjmT2Gt9DBqzwB/GHS/jfT8cm3X6ckEjw8\nONhwscrBSrn1v1IVFTXkEiOtWEUyt/Jhdgy90P430/FRRcfmRb2qLZrZnwN/4+5n5XnNK7ENqa1W\nrAKYu89AbKokNlOscoCqLYqIxISqLYqItDAldBGRmFBCFxGJCSV0EZGYUEIXEYkJJXQRkZhQQo85\nVdo7QMeiMB2beFBCj7E1a9bR1bWE008/j66uJaxZs67eIdWNjkVh69asYUlXF+edfjpLurpYt2ZN\nvUOSEunGophqppoj1aZjUZjqvTQu3VgkE7JVAYMEBrlVAVuNjkVhExUZM89zKzJK81FCj6lkMsno\n6ADBNRfANsbGBifqlLQSHYvCkskkA6OjOUcGBsfGdGyalBJ6TC1YsICenqtJJE5h3rwVJBKn0NNz\ndUv+Ga1jUdiCBQu4uqeHUxIJVsybxymJBFf39OjYNCmNocecKu0doGNRmI5N41G1RRGRmNCHoiIi\nLUwJXUQkJpTQRURiQgldRCQmSk7oZna0mf2XmT1kZtvN7KJKBiYiItGUPMvFzI4EjnT3LWY2F3gA\nONvdH57STrNcREQiquksF3f/o7tvyfy8B+gHjiq1P5FqU0VBibuKjKGbWRI4Cbi/Ev2JVJqqLUor\nKPvGosxwSxr4f+5+c57XNeQidaVqi9KMShlymVXmBmcBNwD/ni+ZZ3V3d0/8nEqlSKVS5WxWJJJs\ntcWRkYOrLSqhS6NIp9Ok0+my+ijrCt3MrgOedvfPTtNGV+hSV7pCl2ZU0w9Fzex/AP8HeKeZbTaz\nB83sjFL7E6kWVVuUVqHiXNIyVFFQmomqLYqIxISqLYqItDAldBGRmFBCFxGJCSV0EZGYUEIXEYkJ\nJXQRkZhQQhcRiQkldBGRmFBCFxGJCSV0EZGYUEIXEYkJJXQRkZhQQhcRiQkldBGRmFBCFxGJCSV0\nEZGYUEIXEYkJJXQRkZgoK6Gb2Rlm9rCZ/dbMLq1UUCIiEl3JCd3MZgDfBd4NnAB82MyWVCqwRpFO\np+sdQlmaOf5mjh0Uf701e/ylKOcKfSXwqLsPuvsYsBY4uzJhNY5mf1M0c/zNHDso/npr9vhLUU5C\nPwrYlfP895llIiJSB/pQVEQkJszdS1vR7K1At7ufkXl+GeDufsWUdqVtQESkxbm7RWlfTkKfCTwC\nnAr8AegDPuzu/SV1KCIiZZlV6oruvt/MLgA2EAzd9CiZi4jUT8lX6CIi0liq+qGomXWa2fVm1m9m\nD5nZW6q5vUoxs+PMbLOZPZj59wUzu6jecUVhZpeY2W/MbJuZ/YeZtdc7pijM7GIz2555NPyxN7Me\nM9ttZttylh1uZhvM7BEzu8PMOusZ43QKxP8XmffQfjNbUc/4iikQ/5WZ3LPFzH5iZvPqGeN0CsT/\nZTPbmslBt5vZkcX6qfYsl+8At7n78cAbgaYYknH337r7cndfAbwJeBm4sc5hhWZmi4ALgRXuvoxg\naO2v6xtVeGZ2AvBJ4GTgJOB/mdnr6xtVUasJbrLLdRlwp7v/KfBfwOU1jyq8fPFvBz4A3FX7cCLL\nF/8G4AR3Pwl4lOY7/le6+xvdfTnwn8CqYp1ULaFnfhu+w91XA7j7Pnd/sVrbq6LTgMfcfVfRlo1l\nJjDHzGYBs4Gn6hxPFMcD97v7XnffD/wS+GCdY5qWu98NPDdl8dnAtZmfrwXeX9OgIsgXv7s/4u6P\nApFmWtTiCUn1AAACWklEQVRDgfjvdPfxzNP7gKNrHlhIBeLfk/N0DjBOEdW8Qj8GeNrMVmeGLv7N\nzBJV3F61fAhYU+8gonD3p4CrgCeAJ4Hn3f3O+kYVyW+Ad2SGLGYD7wUW1zmmUhzh7rsB3P2PwBF1\njqeVnQv8rN5BRGVmXzGzJ4BzgC8Va1/NhD4LWAF8LzN08QrBn6BNw8zagLOA6+sdSxRmdhjB1WEX\nsAiYa2bn1Deq8Nz9YeAK4OfAbcBmYH9dg6oMzUCoAzP7AjDm7r31jiUqd/+iu78O+A+CYdRpVTOh\n/x7Y5e6/zjy/gSDBN5P3AA+4+3C9A4noNGCnuz+bGbJYD7ytzjFF4u6r3f1kd08BzwO/rXNIpdht\nZgsBMh9oDdU5npZjZh8n+AuvaS5oCugF/nexRlVL6Jk/NXeZ2XGZRacCO6q1vSr5ME023JLxBPBW\nM+swMyM49k3xgXSWmS3I/Ps6gg/mmuHqypg83nwL8PHMzx8Dbq51QBFNjX/qa41uUvxmdgbwOeAs\nd99bt6jCmxr/G3Jeez8h/g9XdR66mb0RuAZoA3YCn3D3F6q2wQrKjN0OAq9395fqHU9UZraKYGbL\nGMGQxacyVTGbgpn9EphPEP8l7p6ub0TTM7NeIAW8BthNMCPhJoLhusUE76W/cvfn6xXjdArE/xzw\nT8BrCf5K2uLu76lXjNMpEP/ngXbgmUyz+9z9M3UJsIgC8Z8J/CnBcOMgcJ67/2HafnRjkYhIPKja\noohITCihi4jEhBK6iEhMKKGLiMSEErqISEwooYuIxIQSuohITCihi4jExP8HgN2fv0VnGAkAAAAA\nSUVORK5CYII=\n", | |
| "text/plain": [ | |
| "<matplotlib.figure.Figure at 0x1116ae048>" | |
| ] | |
| }, | |
| "metadata": {}, | |
| "output_type": "display_data" | |
| } | |
| ], | |
| "source": [ | |
| "# get params\n", | |
| "b1,b2,b3 = result.params\n", | |
| "# plot\n", | |
| "fig, ax = plt.subplots(1,1)\n", | |
| "raw_c = ax.scatter(df_c.x, df_c.y) # plot original data C\n", | |
| "raw_t = ax.scatter(df_t.x, df_t.y, c='r') # plot original data T\n", | |
| "fitted_c = ax.plot(df_c.x, df_c.x.map(lambda x:math.exp(b1+b2*x))) # plot the fitted line for C\n", | |
| "fitted_t = ax.plot(df_t.x, df_t.x.map(lambda x:math.exp(b1+b2*x+b3)),c='r') # plot the fitted line for T" | |
| ] | |
| } | |
| ], | |
| "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.5.2" | |
| } | |
| }, | |
| "nbformat": 4, | |
| "nbformat_minor": 0 | |
| } |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment