sage/データ解析のための統計モデリング入門勉強ノート第3章
をテンプレートにして作成
[
トップ
] [
新規
|
一覧
|
単語検索
|
最終更新
|
ヘルプ
]
開始行:
[[FrontPage]]
#contents
2014/03/22からのアクセス回数 &counter;
ここで紹介したSageワークシートは、以下のURLからダウンロー...
http://www15191ue.sakura.ne.jp:8000/home/pub/40/
また、Sageのサーバを公開しているサイト（http://www.sagenb...
アップロードし、実行したり、変更していろいろ動きを試すこ...
* Sageでグラフを再現してみよう：データ解析のための統計モ...
この企画は、雑誌や教科書にでているグラフをSageで再現し、
グラフの意味を理解すると共にSageの使い方をマスターするこ...
前回に続き、
[[データ解析のための統計モデリング入門>http://www.amazon....
（以下、久保本と書きます）
の第３章の例題をSageを使って再現してみます。
数式処理システムSageのノートブックは、計算結果を表示する...
この機会にSageを分析に活用してみてはいかがでしょう。
** 前準備 [#s0b3015d]
最初に必要なライブラリーやパッケージをロードしておきます...
sageへの入力:
#pre{{
# Rの必要なライブラリ
#r("install.packages('ggplot2')")
r('library(ggplot2)')
r('library(jsonlite)')
# python用のパッケージ
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from ggplot import *
# RUtilにRとPandasのデータフレームを相互に変換する関数を...
load(DATA + 'RUtil.py')
}}
** 例題のデータ [#p9909584]
３章の例題に使われているデータは、架空の植物の第i番目の個...
にどう影響するかを調べます。
本の図３．１を以下に引用します。
&ref(Fig3.1.png);
３章のデータは、ネット公開されており、以下の例でもネット...
データをロードしたら、summary関数とvar関数を使ってデータ...
ここで、確認することは、以下の3項目です。
- データの個数と欠損値の有無
- データの平均
- データの分散
sageへの入力:
#pre{{
# 3章のデータをネットから取り込む
d = pd.read_csv('http://hosho.ees.hokudai.ac.jp/~kubo/sta...
}}
*** カテゴリ [#vcf0d954]
肥料の有無（施肥\(f_i\)）は、CとTという文字列で表されてお...
アンケートなどのカテゴリと言った方分かりやすいかも知れま...
PandasのデータフレームもRのData.Frameに劣らず、\(y_i\)は...
sageへの入力:
#pre{{
# どんなデータか調べる
d.head()
}}
#pre{{
y x f
0 6 8.31 C
1 6 9.44 C
2 6 9.50 C
3 12 9.07 C
4 10 10.16 C
[5 rows x 3 columns]
}}
sageへの入力:
#pre{{
# 個数と平均、最大、最小を、標準偏差をみる
d.describe()
}}
#pre{{
y x
count 100.000000 100.000000
mean 7.830000 10.089100
std 2.624881 1.008049
min 2.000000 7.190000
25% 6.000000 9.427500
50% 8.000000 10.155000
75% 10.000000 10.685000
max 15.000000 12.400000
[8 rows x 2 columns]
}}
** データの可視化 [#m8579711]
施肥の有無（C, T）別に種子種\(y_i\)とサイズ\(x_i\)の関係...
ggplotを使うと、このようなグラフも簡単に表示することがで...
sageへの入力:
#pre{{
# F別の分布をみる
ggplot(d, aes(x='x', y='y', color='f')) + geom_point()
}}
#pre{{
<ggplot: (17966489)>
}}
sageへの入力:
#pre{{
ggsave('fig-3.2.png', dpi=50)
}}
#pre{{
Saving 11.0 x 8.0 in image.
}}
&ref(fig-3.2.png);
私は箱ひげ図には馴染みがなく、最初にみたときには驚きまし...
慣れてくるとばらつきを知るには良いプロット方法だと感じる...
この分布の違いは、後で示す施肥別のヒストグラムをみると分...
sageへの入力:
#pre{{
# 箱ひげプロットはggplotにはないので、Rのggplot2でプロッ...
# junk = r('d <- read.csv("http://hosho.ees.hokudai.ac.jp...
# dをRに渡す
PandaDf2RDf(d, 'd')
}}
sageへの入力:
#pre{{
graph = preGraph("fig3.3.png")
r('p <- ggplot(d, aes(x=f, y=y)) + geom_boxplot()')
r('plot(p)')
postGraph(graph)
}}
&ref(fig3.3.png);
*** 可視化から見えてくるもの [#g0d0071e]
施肥別の散布図からは、久保本でも整理してあるとおり、
- サイズxが増加すると共に種子数yが増えているように思える
- 施肥fの効果はないように思われる
のように感じます。
このようにデータを整理することが大切です。Sageでは処理と...
バージョン管理もされているので、どのようにして結果を導い...
*** 施肥の違いをもう少しみてみる [#qd3cc819]
PandasのデータフレームもRのData.Frame同様にデータの絞り込...
施肥のあるもの\(f_T\)を、以下の様にfのフィールドが'T"のも...
#pre{{
d_T = d[d['f'] == 'T']
}}
また、Pandasでは、列（Series）をドット(.)で指定することも...
#pre{{
d_T = d[d.f == 'T']
}}
sageへの入力:
#pre{{
# 施肥の有無の違いをみる
d_T = d[d['f'] == 'T']
print d_T.describe(), d_T.var()
}}
#pre{{
y x
count 50.00000 50.000000
mean 7.88000 10.370600
std 2.65453 0.944307
min 3.00000 8.520000
25% 6.00000 9.757500
50% 7.50000 10.225000
75% 9.00000 10.947500
max 15.00000 12.400000
[8 rows x 2 columns] y 7.046531
x 0.891716
dtype: float64
}}
同様に施肥のないもの\(f_C\)を、抽出します。
平均と分散がほぼ同じであることから、どちらもポアソン分布...
sageへの入力:
#pre{{
d_C = d[d['f'] == 'C']
print d_C.describe(), d_C.var()
}}
#pre{{
y x
count 50.000000 50.000000
mean 7.780000 9.807600
std 2.620874 0.999813
min 2.000000 7.190000
25% 6.000000 9.115000
50% 8.000000 9.880000
75% 10.000000 10.412500
max 12.000000 11.930000
[8 rows x 2 columns] y 6.868980
x 0.999627
dtype: float64
}}
列の度数分布は、以下の様にgroupbyでy毎の頻度を求めること...
sageへの入力:
#pre{{
# 施肥なし（C）の度数を調べる
d_C.groupby('y').size()
}}
#pre{{
y
2 1
3 2
4 3
5 4
6 10
7 1
8 5
9 8
10 9
11 4
12 3
dtype: int64
}}
*** 施肥別のヒストグラム [#ydf741c0]
ggplotを使えば、施肥別のヒストグラムも以下の1行でプロット...
facet_wrapでfを指定することで施肥別のヒストグラムを図化で...
sageへの入力:
#pre{{
# F別のヒストグラムをみる
ggplot(d, aes(x="y")) + geom_histogram(color="grey") + fa...
}}
#pre{{
binwidth defaulted to range/30. Use 'binwidth = x' to adj...
binwidth defaulted to range/30. Use 'binwidth = x' to adj...
<ggplot: (17701093)>
}}
sageへの入力:
#pre{{
ggsave('fig-3.0.png', dpi=50)
}}
#pre{{
Saving 11.0 x 8.0 in image.
}}
&ref(fig-3.0.png);
** ポアソン回帰の統計モデル [#v916c054]
施肥別の平均、分散からも\(個体_i\)の種子数\(y_i\)の確率\(...
以下の式で表します。
$$
p(y_i | \lambda_i) = \frac{\lambda_i^{y_i}exp(-\lambda_i)...
$$
*** \(個体_iの平均種子数\lambda_i\)の式 [#x2fd3984]
久保本の3.4.1では、\(個体_iの平均種子数\lambda_i\)を以下...
$$
\lambda_i = exp( \beta_1 + \beta_2 x_i)
$$
ここで、\(\beta_1, \beta_2\)は、式を決定づけるパラメータ...
どうような曲線になるか示しています。これをSageを使って表...
直線よりも指数関数の表現力があるようにみえます。
sageへの入力:
#pre{{
# fig3.4
p1 = plot(lambda x: exp(-1+0.4*x), [x, -4, 5], rgbcolor='...
p2 = plot(lambda x: exp(-2-0.8*x), [x, -4, 5], rgbcolor='...
(p1 + p2).show(figsize=5)
}}
&ref(sage0.png);
*** 線形予測子とリンク関数 [#ued201eb]
glmで使われるリンク関数と線形予測子の関係を整理すると以下...
線形予測子は、パラメータが線形結合しているもので、ターゲ...
（ここでは\(\lambda_i\)）を線形予測子で表したときに、以下...
$$
\lambda_i = f(線形予測子)
$$
リンク関数は、fの逆関数を使って以下の関係にあります。
$$
f^{-1}(\lambda_i) = (線形予測子)
$$
ですから指数関数で表された\(\lambda_i\)のリンク関数は、対...
$$
log \lambda_i = \beta_1 + \beta_2 x_i
$$
*** ポアソン回帰の実行 [#zbbabd42]
久保本では、モデルから推定される予測の良さを重視しており...
今回のデータに対する対数尤度(Log-Likelihood)は、以下の様...
$$
log L(\beta_1, \beta_2) = \sum_i log \frac{\lambda_i^{y_i...
$$
ここでは、Pythonのstatsmodelsパッケージを使ってGLM（一般...
sageへの入力:
#pre{{
# statsmodelsを使ってglmを計算します
import statsmodels.formula.api as smf
import statsmodels.api as sm
from scipy.stats.stats import pearsonr
}}
*** glmの引数 [#k745ea39]
glmの引数について説明すると、
- fit: 回帰の結果をセットする変数
- y ~ x: モデル式を線形予測子で表現します
- data: dでデータフレームを指定します
- family: sm.families.Poissonで確率分布としてポアソン分布...
- link: sm.families.links.logで対数を指定します
- fit(): 最後にfit関数で回帰を計算します
回帰の結果は、summary2関数で表示します。結果は、
- 対数尤度(Log-Likelihood): -235.39
- 逸脱度(deviance): 84.993
- 切片\(\beta_1\)(Intercept): 1.2917、標準誤差は0.3637
- xの係数\(\beta_2\): 0.0757、標準誤差は0.0356
と求まります。
sageへの入力:
#pre{{
fit = smf.glm('y ~ x', data=d, family=sm.families.Poisson...
fit.summary2()
}}
#pre{{
<class 'statsmodels.iolib.summary2.Summary'>
"""
Results: Generalized linear model
=====================================================
Model: GLM AIC: 474.7725
Link Function: log BIC: -366.3137
Dependent Variable: y Log-Likelihood: -235.39
No. Observations: 100 Deviance: 84.993
Df Model: 1 Pearson chi2: 83.8
Df Residuals: 98 Scale: 1.0000
Method: IRLS
-----------------------------------------------------
Coef. Std.Err. t P>|t| [0.025 0.975]
-----------------------------------------------------
Intercept 1.2917 0.3637 3.5517 0.0004 0.5789 2.0045
x 0.0757 0.0356 2.1251 0.0336 0.0059 0.1454
=====================================================
"""
}}
sageへの入力:
#pre{{
# 逸脱度→deviance, 対数尤度logLik→llf, AIC→aicにセットさ...
print fit.deviance, fit.llf, fit.aic
}}
#pre{{
84.9929964907 -235.38625077 474.77250154
}}
*** 回帰結果のプロット [#wf4a31a4]
λの予測値は、predict関数で取得することができます（fitのnu...
施肥別の分布図に予測値を重ね書きします。
予測値をyhatにセットし、ggplotで以下の様にgeom_lineを追加...
sageへの入力:
#pre{{
# λの予測値nuをyhatにセット
d['yhat'] = fit.predict()
}}
sageへの入力:
#pre{{
# λの予測値と一緒にプロット
ggplot(d, aes(x='x', y='y', color='f')) + geom_point() + ...
}}
#pre{{
<ggplot: (15736953)>
}}
sageへの入力:
#pre{{
ggsave('fig-3.7.png', dpi=50)
}}
#pre{{
Saving 11.0 x 8.0 in image.
}}
&ref(fig-3.7.png);
*** 施肥処理Fをモデルに組み込む [#jabeb95d]
つぎに施肥処理fを使ってポアソン回帰をしてみましょう。xと...
因子型で線形予測子を作る場合、因子の数より1個少ない数の\(...
\(d_i\)が、\(f_i\)=Tの場合1，\(f_i!=T\)の場合0とします。
因子型ではいずれか一つの因子に当てはまる（昔のカーラジオ...
としているので、\(f_i\)がTでなければCとなるので、因子がT...
もし、肥料が複数の場合には、ダミー変数は、\(d_{i,A}, d_{i...
glmは、とてもお利口なので人がダミー変数を指定する代わりに...
モデル式をy ~ xからy ~ fに変えてポアソン回帰を実行してみ...
回帰の結果は、summary2関数で表示します。結果は、
- 対数尤度(Log-Likelihood): -237.63
- 逸脱度(deviance): 89.475
- 切片\(\beta_1\)(Intercept): 2.0516、標準誤差は0.0507
- fの係数\(\beta_3\): 0.0128、標準誤差は0.0715
先の例と比べると、対数尤度は-237.63と小さく求まっています。
久保本は、とても丁寧でこの結果をどう考えるかを解説してく...
$$
\lambda_i = exp(2.05 + 0.0128) = exp(2.05) exp(0.0128)
$$
施肥を行った場合には、施肥を行わない場合（0.0128の項が０...
$$
\lambda_i = exp(2.05 + 0.0128) = 1.0128 exp(2.05)
$$
sageへの入力:
#pre{{
# 施肥処理Fをモデルに組み込む
fit_f= smf.glm('y ~ f', data=d, family=sm.families.Poisso...
fit_f.summary2()
}}
#pre{{
<class 'statsmodels.iolib.summary2.Summary'>
"""
Results: Generalized linear model
=======================================================
Model: GLM AIC: 479.2545
Link Function: log BIC: -361.8317
Dependent Variable: y Log-Likelihood: -237.63
No. Observations: 100 Deviance: 89.475
Df Model: 1 Pearson chi2: 87.1
Df Residuals: 98 Scale: 1.0000
Method: IRLS
-------------------------------------------------------
Coef. Std.Err. t P>|t| [0.025 0.975]
-------------------------------------------------------
Intercept 2.0516 0.0507 40.4630 0.0000 1.9522 2.1509
f[T.T] 0.0128 0.0715 0.1787 0.8582 -0.1273 0.1529
=======================================================
"""
}}
** 個体のサイズと施肥を合わせたモデル [#w4a9a065]
最後に個体のサイズ\(x_i\)と施肥\(f_i\)を合わせたモデルを...
モデル式は、y ~ x + fとなります。
結果は以下の様になります。
$$
\lambda_i = 1.26 + 0.08 x_i - 0.032
$$
sageへの入力:
#pre{{
# 説明変数が数量型＋因子型のモデル
fit_all = smf.glm('y ~ x + f', data=d, family=sm.families...
fit_all.summary2()
}}
#pre{{
<class 'statsmodels.iolib.summary2.Summary'>
"""
Results: Generalized linear model
========================================================
Model: GLM AIC: 476.5874
Link Function: log BIC: -361.8936
Dependent Variable: y Log-Likelihood: -235.29
No. Observations: 100 Deviance: 84.808
Df Model: 2 Pearson chi2: 83.8
Df Residuals: 97 Scale: 1.0000
Method: IRLS
--------------------------------------------------------
Coef. Std.Err. t P>|t| [0.025 0.975]
--------------------------------------------------------
Intercept 1.2631 0.3696 3.4172 0.0006 0.5386 1.9876
f[T.T] -0.0320 0.0744 -0.4302 0.6670 -0.1778 0.1138
x 0.0801 0.0370 2.1620 0.0306 0.0075 0.1527
========================================================
"""
}}
*** Nullモデルとの比較 [#r47c8714]
ついでなので、λが定数の場合（Nullモデル）のポアソン回帰の...
ここで注目するのは、Log-Likelihoodの値-237.64です。
λが定数の場合の対数尤度と自分の考えたモデルの対数尤度を比...
自分のモデルがどの程度よいのか比較できることです。
sageへの入力:
#pre{{
# Nullモデルの逸脱度とAIC
fit_null = smf.glm('y ~ 1', data=d, family=sm.families.Po...
fit_null.summary2()
}}
#pre{{
<class 'statsmodels.iolib.summary2.Summary'>
"""
Results: Generalized linear model
======================================================
Model: GLM AIC: 477.2864
Link Function: log BIC: -366.4049
Dependent Variable: y Log-Likelihood: -237.64
No. Observations: 100 Deviance: 89.507
Df Model: 0 Pearson chi2: 87.1
Df Residuals: 99 Scale: 1.0000
Method: IRLS
------------------------------------------------------
Coef. Std.Err. t P>|t| [0.025 0.975]
------------------------------------------------------
Intercept 2.0580 0.0357 57.5862 0.0000 1.9879 2.1280
======================================================
"""
}}
sageへの入力:
#pre{{
# 逸脱度→deviance, 対数尤度logLik→llf, AIC→aicにセットさ...
print fit_null.deviance, fit_null.llf, fit_null.aic
}}
#pre{{
89.5069375696 -237.643221309 477.286442619
}}
** 3章注目のグラフ [#j416a0f2]
「Sageでグラフを再現してみよう」ではいろんなグラフをSage...
3章の注目のグラフは、図3.9（以下久保本から引用します）の...
種子数λがサイズxに対して指数関数で増加するとそのポアソン...
今回はこのグラフをSageを使って計算してみます。
&ref(Fig3.9.png);
GLMの結果が示されていないので、\(x\)が0.5, 1.1, 1.7のyの...
0.1, 0.5, 3.0としてポアソン分布をプロットしたのが、以下の...
とよく似た分布になっています。
sageへの入力:
#pre{{
# 3章のメインのグラフは、図3.9だと思う。
# 0.5, 1.1, 1.7のyを0.1, 0.5, 3.0としてポアソン分布を表示
# Fig. 2.6
# ポアソン関数p(y, λ)を定義
def p(y, lam):
return lam^y*e^(-lam)/factorial(y)
g = Graphics()
for lam, cls in [(0.1, 'red'), (0.5, 'blue'), (3.0, 'gree...
g = g + list_plot([p(y, lam) for y in range(8)], colo...
g.show(figsize=5)
}}
&ref(sage1.png);
** コメント [#l0dba599]
#vote(おもしろかった[1],そうでもない[0],わかりずらい[0])
皆様のご意見、ご希望をお待ちしております。
#comment_kcaptcha
終了行:
[[FrontPage]]
#contents
2014/03/22からのアクセス回数 &counter;
ここで紹介したSageワークシートは、以下のURLからダウンロー...
http://www15191ue.sakura.ne.jp:8000/home/pub/40/
また、Sageのサーバを公開しているサイト（http://www.sagenb...
アップロードし、実行したり、変更していろいろ動きを試すこ...
* Sageでグラフを再現してみよう：データ解析のための統計モ...
この企画は、雑誌や教科書にでているグラフをSageで再現し、
グラフの意味を理解すると共にSageの使い方をマスターするこ...
前回に続き、
[[データ解析のための統計モデリング入門>http://www.amazon....
（以下、久保本と書きます）
の第３章の例題をSageを使って再現してみます。
数式処理システムSageのノートブックは、計算結果を表示する...
この機会にSageを分析に活用してみてはいかがでしょう。
** 前準備 [#s0b3015d]
最初に必要なライブラリーやパッケージをロードしておきます...
sageへの入力:
#pre{{
# Rの必要なライブラリ
#r("install.packages('ggplot2')")
r('library(ggplot2)')
r('library(jsonlite)')
# python用のパッケージ
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from ggplot import *
# RUtilにRとPandasのデータフレームを相互に変換する関数を...
load(DATA + 'RUtil.py')
}}
** 例題のデータ [#p9909584]
３章の例題に使われているデータは、架空の植物の第i番目の個...
にどう影響するかを調べます。
本の図３．１を以下に引用します。
&ref(Fig3.1.png);
３章のデータは、ネット公開されており、以下の例でもネット...
データをロードしたら、summary関数とvar関数を使ってデータ...
ここで、確認することは、以下の3項目です。
- データの個数と欠損値の有無
- データの平均
- データの分散
sageへの入力:
#pre{{
# 3章のデータをネットから取り込む
d = pd.read_csv('http://hosho.ees.hokudai.ac.jp/~kubo/sta...
}}
*** カテゴリ [#vcf0d954]
肥料の有無（施肥\(f_i\)）は、CとTという文字列で表されてお...
アンケートなどのカテゴリと言った方分かりやすいかも知れま...
PandasのデータフレームもRのData.Frameに劣らず、\(y_i\)は...
sageへの入力:
#pre{{
# どんなデータか調べる
d.head()
}}
#pre{{
y x f
0 6 8.31 C
1 6 9.44 C
2 6 9.50 C
3 12 9.07 C
4 10 10.16 C
[5 rows x 3 columns]
}}
sageへの入力:
#pre{{
# 個数と平均、最大、最小を、標準偏差をみる
d.describe()
}}
#pre{{
y x
count 100.000000 100.000000
mean 7.830000 10.089100
std 2.624881 1.008049
min 2.000000 7.190000
25% 6.000000 9.427500
50% 8.000000 10.155000
75% 10.000000 10.685000
max 15.000000 12.400000
[8 rows x 2 columns]
}}
** データの可視化 [#m8579711]
施肥の有無（C, T）別に種子種\(y_i\)とサイズ\(x_i\)の関係...
ggplotを使うと、このようなグラフも簡単に表示することがで...
sageへの入力:
#pre{{
# F別の分布をみる
ggplot(d, aes(x='x', y='y', color='f')) + geom_point()
}}
#pre{{
<ggplot: (17966489)>
}}
sageへの入力:
#pre{{
ggsave('fig-3.2.png', dpi=50)
}}
#pre{{
Saving 11.0 x 8.0 in image.
}}
&ref(fig-3.2.png);
私は箱ひげ図には馴染みがなく、最初にみたときには驚きまし...
慣れてくるとばらつきを知るには良いプロット方法だと感じる...
この分布の違いは、後で示す施肥別のヒストグラムをみると分...
sageへの入力:
#pre{{
# 箱ひげプロットはggplotにはないので、Rのggplot2でプロッ...
# junk = r('d <- read.csv("http://hosho.ees.hokudai.ac.jp...
# dをRに渡す
PandaDf2RDf(d, 'd')
}}
sageへの入力:
#pre{{
graph = preGraph("fig3.3.png")
r('p <- ggplot(d, aes(x=f, y=y)) + geom_boxplot()')
r('plot(p)')
postGraph(graph)
}}
&ref(fig3.3.png);
*** 可視化から見えてくるもの [#g0d0071e]
施肥別の散布図からは、久保本でも整理してあるとおり、
- サイズxが増加すると共に種子数yが増えているように思える
- 施肥fの効果はないように思われる
のように感じます。
このようにデータを整理することが大切です。Sageでは処理と...
バージョン管理もされているので、どのようにして結果を導い...
*** 施肥の違いをもう少しみてみる [#qd3cc819]
PandasのデータフレームもRのData.Frame同様にデータの絞り込...
施肥のあるもの\(f_T\)を、以下の様にfのフィールドが'T"のも...
#pre{{
d_T = d[d['f'] == 'T']
}}
また、Pandasでは、列（Series）をドット(.)で指定することも...
#pre{{
d_T = d[d.f == 'T']
}}
sageへの入力:
#pre{{
# 施肥の有無の違いをみる
d_T = d[d['f'] == 'T']
print d_T.describe(), d_T.var()
}}
#pre{{
y x
count 50.00000 50.000000
mean 7.88000 10.370600
std 2.65453 0.944307
min 3.00000 8.520000
25% 6.00000 9.757500
50% 7.50000 10.225000
75% 9.00000 10.947500
max 15.00000 12.400000
[8 rows x 2 columns] y 7.046531
x 0.891716
dtype: float64
}}
同様に施肥のないもの\(f_C\)を、抽出します。
平均と分散がほぼ同じであることから、どちらもポアソン分布...
sageへの入力:
#pre{{
d_C = d[d['f'] == 'C']
print d_C.describe(), d_C.var()
}}
#pre{{
y x
count 50.000000 50.000000
mean 7.780000 9.807600
std 2.620874 0.999813
min 2.000000 7.190000
25% 6.000000 9.115000
50% 8.000000 9.880000
75% 10.000000 10.412500
max 12.000000 11.930000
[8 rows x 2 columns] y 6.868980
x 0.999627
dtype: float64
}}
列の度数分布は、以下の様にgroupbyでy毎の頻度を求めること...
sageへの入力:
#pre{{
# 施肥なし（C）の度数を調べる
d_C.groupby('y').size()
}}
#pre{{
y
2 1
3 2
4 3
5 4
6 10
7 1
8 5
9 8
10 9
11 4
12 3
dtype: int64
}}
*** 施肥別のヒストグラム [#ydf741c0]
ggplotを使えば、施肥別のヒストグラムも以下の1行でプロット...
facet_wrapでfを指定することで施肥別のヒストグラムを図化で...
sageへの入力:
#pre{{
# F別のヒストグラムをみる
ggplot(d, aes(x="y")) + geom_histogram(color="grey") + fa...
}}
#pre{{
binwidth defaulted to range/30. Use 'binwidth = x' to adj...
binwidth defaulted to range/30. Use 'binwidth = x' to adj...
<ggplot: (17701093)>
}}
sageへの入力:
#pre{{
ggsave('fig-3.0.png', dpi=50)
}}
#pre{{
Saving 11.0 x 8.0 in image.
}}
&ref(fig-3.0.png);
** ポアソン回帰の統計モデル [#v916c054]
施肥別の平均、分散からも\(個体_i\)の種子数\(y_i\)の確率\(...
以下の式で表します。
$$
p(y_i | \lambda_i) = \frac{\lambda_i^{y_i}exp(-\lambda_i)...
$$
*** \(個体_iの平均種子数\lambda_i\)の式 [#x2fd3984]
久保本の3.4.1では、\(個体_iの平均種子数\lambda_i\)を以下...
$$
\lambda_i = exp( \beta_1 + \beta_2 x_i)
$$
ここで、\(\beta_1, \beta_2\)は、式を決定づけるパラメータ...
どうような曲線になるか示しています。これをSageを使って表...
直線よりも指数関数の表現力があるようにみえます。
sageへの入力:
#pre{{
# fig3.4
p1 = plot(lambda x: exp(-1+0.4*x), [x, -4, 5], rgbcolor='...
p2 = plot(lambda x: exp(-2-0.8*x), [x, -4, 5], rgbcolor='...
(p1 + p2).show(figsize=5)
}}
&ref(sage0.png);
*** 線形予測子とリンク関数 [#ued201eb]
glmで使われるリンク関数と線形予測子の関係を整理すると以下...
線形予測子は、パラメータが線形結合しているもので、ターゲ...
（ここでは\(\lambda_i\)）を線形予測子で表したときに、以下...
$$
\lambda_i = f(線形予測子)
$$
リンク関数は、fの逆関数を使って以下の関係にあります。
$$
f^{-1}(\lambda_i) = (線形予測子)
$$
ですから指数関数で表された\(\lambda_i\)のリンク関数は、対...
$$
log \lambda_i = \beta_1 + \beta_2 x_i
$$
*** ポアソン回帰の実行 [#zbbabd42]
久保本では、モデルから推定される予測の良さを重視しており...
今回のデータに対する対数尤度(Log-Likelihood)は、以下の様...
$$
log L(\beta_1, \beta_2) = \sum_i log \frac{\lambda_i^{y_i...
$$
ここでは、Pythonのstatsmodelsパッケージを使ってGLM（一般...
sageへの入力:
#pre{{
# statsmodelsを使ってglmを計算します
import statsmodels.formula.api as smf
import statsmodels.api as sm
from scipy.stats.stats import pearsonr
}}
*** glmの引数 [#k745ea39]
glmの引数について説明すると、
- fit: 回帰の結果をセットする変数
- y ~ x: モデル式を線形予測子で表現します
- data: dでデータフレームを指定します
- family: sm.families.Poissonで確率分布としてポアソン分布...
- link: sm.families.links.logで対数を指定します
- fit(): 最後にfit関数で回帰を計算します
回帰の結果は、summary2関数で表示します。結果は、
- 対数尤度(Log-Likelihood): -235.39
- 逸脱度(deviance): 84.993
- 切片\(\beta_1\)(Intercept): 1.2917、標準誤差は0.3637
- xの係数\(\beta_2\): 0.0757、標準誤差は0.0356
と求まります。
sageへの入力:
#pre{{
fit = smf.glm('y ~ x', data=d, family=sm.families.Poisson...
fit.summary2()
}}
#pre{{
<class 'statsmodels.iolib.summary2.Summary'>
"""
Results: Generalized linear model
=====================================================
Model: GLM AIC: 474.7725
Link Function: log BIC: -366.3137
Dependent Variable: y Log-Likelihood: -235.39
No. Observations: 100 Deviance: 84.993
Df Model: 1 Pearson chi2: 83.8
Df Residuals: 98 Scale: 1.0000
Method: IRLS
-----------------------------------------------------
Coef. Std.Err. t P>|t| [0.025 0.975]
-----------------------------------------------------
Intercept 1.2917 0.3637 3.5517 0.0004 0.5789 2.0045
x 0.0757 0.0356 2.1251 0.0336 0.0059 0.1454
=====================================================
"""
}}
sageへの入力:
#pre{{
# 逸脱度→deviance, 対数尤度logLik→llf, AIC→aicにセットさ...
print fit.deviance, fit.llf, fit.aic
}}
#pre{{
84.9929964907 -235.38625077 474.77250154
}}
*** 回帰結果のプロット [#wf4a31a4]
λの予測値は、predict関数で取得することができます（fitのnu...
施肥別の分布図に予測値を重ね書きします。
予測値をyhatにセットし、ggplotで以下の様にgeom_lineを追加...
sageへの入力:
#pre{{
# λの予測値nuをyhatにセット
d['yhat'] = fit.predict()
}}
sageへの入力:
#pre{{
# λの予測値と一緒にプロット
ggplot(d, aes(x='x', y='y', color='f')) + geom_point() + ...
}}
#pre{{
<ggplot: (15736953)>
}}
sageへの入力:
#pre{{
ggsave('fig-3.7.png', dpi=50)
}}
#pre{{
Saving 11.0 x 8.0 in image.
}}
&ref(fig-3.7.png);
*** 施肥処理Fをモデルに組み込む [#jabeb95d]
つぎに施肥処理fを使ってポアソン回帰をしてみましょう。xと...
因子型で線形予測子を作る場合、因子の数より1個少ない数の\(...
\(d_i\)が、\(f_i\)=Tの場合1，\(f_i!=T\)の場合0とします。
因子型ではいずれか一つの因子に当てはまる（昔のカーラジオ...
としているので、\(f_i\)がTでなければCとなるので、因子がT...
もし、肥料が複数の場合には、ダミー変数は、\(d_{i,A}, d_{i...
glmは、とてもお利口なので人がダミー変数を指定する代わりに...
モデル式をy ~ xからy ~ fに変えてポアソン回帰を実行してみ...
回帰の結果は、summary2関数で表示します。結果は、
- 対数尤度(Log-Likelihood): -237.63
- 逸脱度(deviance): 89.475
- 切片\(\beta_1\)(Intercept): 2.0516、標準誤差は0.0507
- fの係数\(\beta_3\): 0.0128、標準誤差は0.0715
先の例と比べると、対数尤度は-237.63と小さく求まっています。
久保本は、とても丁寧でこの結果をどう考えるかを解説してく...
$$
\lambda_i = exp(2.05 + 0.0128) = exp(2.05) exp(0.0128)
$$
施肥を行った場合には、施肥を行わない場合（0.0128の項が０...
$$
\lambda_i = exp(2.05 + 0.0128) = 1.0128 exp(2.05)
$$
sageへの入力:
#pre{{
# 施肥処理Fをモデルに組み込む
fit_f= smf.glm('y ~ f', data=d, family=sm.families.Poisso...
fit_f.summary2()
}}
#pre{{
<class 'statsmodels.iolib.summary2.Summary'>
"""
Results: Generalized linear model
=======================================================
Model: GLM AIC: 479.2545
Link Function: log BIC: -361.8317
Dependent Variable: y Log-Likelihood: -237.63
No. Observations: 100 Deviance: 89.475
Df Model: 1 Pearson chi2: 87.1
Df Residuals: 98 Scale: 1.0000
Method: IRLS
-------------------------------------------------------
Coef. Std.Err. t P>|t| [0.025 0.975]
-------------------------------------------------------
Intercept 2.0516 0.0507 40.4630 0.0000 1.9522 2.1509
f[T.T] 0.0128 0.0715 0.1787 0.8582 -0.1273 0.1529
=======================================================
"""
}}
** 個体のサイズと施肥を合わせたモデル [#w4a9a065]
最後に個体のサイズ\(x_i\)と施肥\(f_i\)を合わせたモデルを...
モデル式は、y ~ x + fとなります。
結果は以下の様になります。
$$
\lambda_i = 1.26 + 0.08 x_i - 0.032
$$
sageへの入力:
#pre{{
# 説明変数が数量型＋因子型のモデル
fit_all = smf.glm('y ~ x + f', data=d, family=sm.families...
fit_all.summary2()
}}
#pre{{
<class 'statsmodels.iolib.summary2.Summary'>
"""
Results: Generalized linear model
========================================================
Model: GLM AIC: 476.5874
Link Function: log BIC: -361.8936
Dependent Variable: y Log-Likelihood: -235.29
No. Observations: 100 Deviance: 84.808
Df Model: 2 Pearson chi2: 83.8
Df Residuals: 97 Scale: 1.0000
Method: IRLS
--------------------------------------------------------
Coef. Std.Err. t P>|t| [0.025 0.975]
--------------------------------------------------------
Intercept 1.2631 0.3696 3.4172 0.0006 0.5386 1.9876
f[T.T] -0.0320 0.0744 -0.4302 0.6670 -0.1778 0.1138
x 0.0801 0.0370 2.1620 0.0306 0.0075 0.1527
========================================================
"""
}}
*** Nullモデルとの比較 [#r47c8714]
ついでなので、λが定数の場合（Nullモデル）のポアソン回帰の...
ここで注目するのは、Log-Likelihoodの値-237.64です。
λが定数の場合の対数尤度と自分の考えたモデルの対数尤度を比...
自分のモデルがどの程度よいのか比較できることです。
sageへの入力:
#pre{{
# Nullモデルの逸脱度とAIC
fit_null = smf.glm('y ~ 1', data=d, family=sm.families.Po...
fit_null.summary2()
}}
#pre{{
<class 'statsmodels.iolib.summary2.Summary'>
"""
Results: Generalized linear model
======================================================
Model: GLM AIC: 477.2864
Link Function: log BIC: -366.4049
Dependent Variable: y Log-Likelihood: -237.64
No. Observations: 100 Deviance: 89.507
Df Model: 0 Pearson chi2: 87.1
Df Residuals: 99 Scale: 1.0000
Method: IRLS
------------------------------------------------------
Coef. Std.Err. t P>|t| [0.025 0.975]
------------------------------------------------------
Intercept 2.0580 0.0357 57.5862 0.0000 1.9879 2.1280
======================================================
"""
}}
sageへの入力:
#pre{{
# 逸脱度→deviance, 対数尤度logLik→llf, AIC→aicにセットさ...
print fit_null.deviance, fit_null.llf, fit_null.aic
}}
#pre{{
89.5069375696 -237.643221309 477.286442619
}}
** 3章注目のグラフ [#j416a0f2]
「Sageでグラフを再現してみよう」ではいろんなグラフをSage...
3章の注目のグラフは、図3.9（以下久保本から引用します）の...
種子数λがサイズxに対して指数関数で増加するとそのポアソン...
今回はこのグラフをSageを使って計算してみます。
&ref(Fig3.9.png);
GLMの結果が示されていないので、\(x\)が0.5, 1.1, 1.7のyの...
0.1, 0.5, 3.0としてポアソン分布をプロットしたのが、以下の...
とよく似た分布になっています。
sageへの入力:
#pre{{
# 3章のメインのグラフは、図3.9だと思う。
# 0.5, 1.1, 1.7のyを0.1, 0.5, 3.0としてポアソン分布を表示
# Fig. 2.6
# ポアソン関数p(y, λ)を定義
def p(y, lam):
return lam^y*e^(-lam)/factorial(y)
g = Graphics()
for lam, cls in [(0.1, 'red'), (0.5, 'blue'), (3.0, 'gree...
g = g + list_plot([p(y, lam) for y in range(8)], colo...
g.show(figsize=5)
}}
&ref(sage1.png);
** コメント [#l0dba599]
#vote(おもしろかった[1],そうでもない[0],わかりずらい[0])
皆様のご意見、ご希望をお待ちしております。
#comment_kcaptcha
ページ名:
SmartDoc
![[PukiWiki]](image/pukiwiki.png "[PukiWiki]")
