sage/データ解析のための統計モデリング入門勉強ノート第5章
をテンプレートにして作成
[
トップ
] [
新規
|
一覧
|
単語検索
|
最終更新
|
ヘルプ
]
開始行:
[[FrontPage]]
#contents
2014/04/06からのアクセス回数 &counter;
ここで紹介したSageワークシートは、以下のURLからダウンロー...
http://www15191ue.sakura.ne.jp:8000/home/pub/42/
また、Sageのサーバを公開しているサイト（http://www.sagenb...
アップロードし、実行したり、変更していろいろ動きを試すこ...
* Sageでグラフを再現してみよう：データ解析のための統計モ...
この企画は、雑誌や教科書にでているグラフをSageで再現し、
グラフの意味を理解すると共にSageの使い方をマスターするこ...
前回に続き、
[[データ解析のための統計モデリング入門>http://www.amazon....
（以下、久保本と書きます）
の第5章の例題をSageを使って再現してみます。
私が統計を嫌いになったのは、χ二乗検定が原因です。
どうしてそうなるのかも説明されないまま検定をすることに納...
久保本はそれをすぱっと解決してくれました。今回の目標はず...
&ref(Fig_5.4.png);
数式処理システムSageのノートブックは、計算結果を表示する...
この機会にSageを分析に活用してみてはいかがでしょう。
** 前準備 [#e6f6abc9]
最初に必要なライブラリーやパッケージをロードしておきます。
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 *
# statsmodelsを使ってglmを計算します
import statsmodels.formula.api as smf
import statsmodels.api as sm
from scipy.stats.stats import pearsonr
# RUtilにRとPandasのデータフレームを相互に変換する関数を...
load(DATA + 'RUtil.py')
}}
** 尤度比検定がすごい [#xe547c53]
どんな統計モデルにも使える尤度比検定が、5章のテーマです。...
：帰無仮説の尤度比を使って帰無仮説棄却の可否を判断します。
((久保本では「無に帰される」から帰無仮説と説明があり昔の...
今回は、3章の例題で一定モデルとｘモデルを比較します。対数...
$$
\frac{L^*_1}{L^*_2} = \frac{一定モデルの最大尤度}{ｘモデ...
$$
これの対数に-2を掛けた値（逸脱度の差）を使って検定をしま...
$$
\Delta D_{1,2} = -2 \times ( log L^*_1 - log L^*_2)
$$
逸脱度の差は、以下の計算で4.5とでました。
sageへの入力:
#pre{{
# 3章のデータをもう一度読み込む
d = pd.read_csv('http://hosho.ees.hokudai.ac.jp/~kubo/sta...
}}
sageへの入力:
#pre{{
# 一定モデル（k=1）でのglm回帰を実行
fit1 = smf.glm('y ~ 1', data=d, family=sm.families.Poisso...
}}
sageへの入力:
#pre{{
# xモデル
fit2 = smf.glm('y ~ x', data=d, family=sm.families.Poisso...
}}
sageへの入力:
#pre{{
# 逸脱度（-2logL）
-2*(fit1.llf - fit2.llf)
}}
#pre{{
4.5139410788522696
}}
sageへの入力:
#pre{{
# 切片beta_2の推定値2.058となり、exp(2.058)がd.yの平均と...
print fit1.summary2()
print exp(2.058), d.y.mean()
}}
#pre{{
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
======================================================
7.83029355218137 7.83
}}
** 帰無仮説を棄却させる手順 [#mebe844d]
第一種の過誤の検討に専念して、以下の手順で帰無仮説を棄却...
+ 帰無仮説（今回は一定モデル）が正しいものと仮定する
+ 観測データに帰無モデルをあてはめ、\(\hat{\beta_1}=2.058...
+ 上記の真のモデルからサンプルデータを生成し、その度に一...
+ 求まった分布から一定モデル（帰無モデル）とxモデルの逸脱...
*** サンプルの生成 [#s77e1737]
Sageには、ポアソン分布を生成する関数が無いので、Rのrpois...
sample_size個生成する関数genSampleを以下の様に定義します。
sageへの入力:
#pre{{
# sageobjを使って変換する方が速いのでgenSample関数を使用...
def genSample(lambda_target, sample_size):
return np.array(sageobj(r('rpois(%d, %f)'%(sample_siz...
}}
*** パラメトリックブートストラップ法 [#w46fc378]
パラメトリックブートストラップ法をSageで試してみます。4章...
一定モデルを真のモデルとしてサンプルデータを生成し、一定...
回帰分析を行い、対数尤度llfから逸脱度Dの差を求めます。
これをn_bootstrap回繰り返し、その分布を求めます。
sageへの入力:
#pre{{
# パラメトリックブートストラップ法
def pb(d, n_bootstrap):
N = len(d)
y_mean = d.y.mean()
def sampling(d, y_mean, N):
# 帰無仮説を真のモデルとしてサンプルデータを生成
d['rnd'] = genSample(y_mean, N)
fit1 = smf.glm('rnd ~ 1', data=d, family=sm.famil...
fit2 = smf.glm('rnd ~ x', data=d, family=sm.famil...
return -2*(fit1.llf - fit2.llf)
return [ sampling(d, y_mean, N) for i in range(n_boot...
}}
sageへの入力:
#pre{{
# 逸脱度の差を1000サンプル生成
dd12 = pb(d, 1000)
}}
*** サンプルの結果 [#b0f88167]
ブートストラップ法で求まった結果にdescribe（Rのsummaryに...
ggplotのヒストグラムでその分布を表示してみます。
sageへの入力:
#pre{{
dd12Df = pd.DataFrame({'x' : dd12})
dd12Df.describe()
}}
#pre{{
x
count 1.000000e+03
mean 1.016127e+00
std 1.438180e+00
min 2.531954e-08
25% 9.873545e-02
50% 4.598125e-01
75% 1.367469e+00
max 1.048328e+01
[8 rows x 1 columns]
}}
sageへの入力:
#pre{{
ggplot(dd12Df, aes(x="x")) + geom_histogram(color="grey")...
}}
#pre{{
binwidth defaulted to range/30. Use 'binwidth = x' to adj...
<ggplot: (13280017)>
}}
sageへの入力:
#pre{{
ggsave('fig-5.4.png', dpi=50)
}}
#pre{{
Saving 11.0 x 8.0 in image.
}}
&ref(fig-5.4.png);
*** 仮説の検証 [#y8954cad]
4.5以上のデータの個数は、37個でp値は、37/1000=0.037となり...
3.90と求まりました。
これで、有意水準5%の検定において帰無仮説（一定モデル）は...
sageへの入力:
#pre{{
# 4.5以上のデータの個数は
len(dd12Df[dd12Df.x >= 4.5])
}}
#pre{{
37
}}
sageへの入力:
#pre{{
# 95%となるxの値を求める
dd12Df.quantile(0.95)
}}
#pre{{
x 3.899619
dtype: float64
}}
** χ自乗分布を使かう方法 [#r56bd0ba]
残念ながら、statsmodelsにはglmに対するanovaがないため、R...
p値は、0.034となり、同様に帰無仮説は棄却されました。（Rを...
sageへの入力:
#pre{{
# χ自乗分布を使かう方法
fit1 = smf.glm('y ~ 1', data=d, family=sm.families.Poisso...
fit2 = smf.glm('y ~ x', data=d, family=sm.families.Poisso...
}}
sageへの入力:
#pre{{
# 残念ながらstatsmodelsにはglmに対するanovaがないみたい
from statsmodels.stats.anova import anova_lm
anova_lm(fit1, fit2)
}}
#pre{{
Traceback (most recent call last):
File "<stdin>", line 1, in <module>
File "_sage_input_80.py", line 10, in <module>
exec compile(u'open("___code___.py","w").write("# -*-...
途中省略
35, in __getattribute__
obj = getattr(results, attr)
AttributeError: 'GLMResults' object has no attribute 'ssr'
}}
sageへの入力:
#pre{{
# Rにｄを渡して計算
PandaDf2RDf(d, "d")
}}
sageへの入力:
#pre{{
# ANOVA(ANalysis Of VAriance)を使って逸脱度を調べる
r('fit1 <- glm(y ~ 1, data=d, family=poisson)')
r('fit2 <- glm(y ~ x, data=d, family=poisson)')
r('anova(fit1, fit2, test = "Chisq")')
}}
#pre{{
Analysis of Deviance Table
Model 1: y ~ 1
Model 2: y ~ x
Resid. Df Resid. Dev Df Deviance Pr(>Chi)
1 99 89.507
2 98 84.993 1 4.5139 0.03362 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ...
}}
** コメント [#n4fa83ed]
#vote(おもしろかった[1],そうでもない[0],わかりずらい[0])
皆様のご意見、ご希望をお待ちしております。
#comment_kcaptcha
終了行:
[[FrontPage]]
#contents
2014/04/06からのアクセス回数 &counter;
ここで紹介したSageワークシートは、以下のURLからダウンロー...
http://www15191ue.sakura.ne.jp:8000/home/pub/42/
また、Sageのサーバを公開しているサイト（http://www.sagenb...
アップロードし、実行したり、変更していろいろ動きを試すこ...
* Sageでグラフを再現してみよう：データ解析のための統計モ...
この企画は、雑誌や教科書にでているグラフをSageで再現し、
グラフの意味を理解すると共にSageの使い方をマスターするこ...
前回に続き、
[[データ解析のための統計モデリング入門>http://www.amazon....
（以下、久保本と書きます）
の第5章の例題をSageを使って再現してみます。
私が統計を嫌いになったのは、χ二乗検定が原因です。
どうしてそうなるのかも説明されないまま検定をすることに納...
久保本はそれをすぱっと解決してくれました。今回の目標はず...
&ref(Fig_5.4.png);
数式処理システムSageのノートブックは、計算結果を表示する...
この機会にSageを分析に活用してみてはいかがでしょう。
** 前準備 [#e6f6abc9]
最初に必要なライブラリーやパッケージをロードしておきます。
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 *
# statsmodelsを使ってglmを計算します
import statsmodels.formula.api as smf
import statsmodels.api as sm
from scipy.stats.stats import pearsonr
# RUtilにRとPandasのデータフレームを相互に変換する関数を...
load(DATA + 'RUtil.py')
}}
** 尤度比検定がすごい [#xe547c53]
どんな統計モデルにも使える尤度比検定が、5章のテーマです。...
：帰無仮説の尤度比を使って帰無仮説棄却の可否を判断します。
((久保本では「無に帰される」から帰無仮説と説明があり昔の...
今回は、3章の例題で一定モデルとｘモデルを比較します。対数...
$$
\frac{L^*_1}{L^*_2} = \frac{一定モデルの最大尤度}{ｘモデ...
$$
これの対数に-2を掛けた値（逸脱度の差）を使って検定をしま...
$$
\Delta D_{1,2} = -2 \times ( log L^*_1 - log L^*_2)
$$
逸脱度の差は、以下の計算で4.5とでました。
sageへの入力:
#pre{{
# 3章のデータをもう一度読み込む
d = pd.read_csv('http://hosho.ees.hokudai.ac.jp/~kubo/sta...
}}
sageへの入力:
#pre{{
# 一定モデル（k=1）でのglm回帰を実行
fit1 = smf.glm('y ~ 1', data=d, family=sm.families.Poisso...
}}
sageへの入力:
#pre{{
# xモデル
fit2 = smf.glm('y ~ x', data=d, family=sm.families.Poisso...
}}
sageへの入力:
#pre{{
# 逸脱度（-2logL）
-2*(fit1.llf - fit2.llf)
}}
#pre{{
4.5139410788522696
}}
sageへの入力:
#pre{{
# 切片beta_2の推定値2.058となり、exp(2.058)がd.yの平均と...
print fit1.summary2()
print exp(2.058), d.y.mean()
}}
#pre{{
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
======================================================
7.83029355218137 7.83
}}
** 帰無仮説を棄却させる手順 [#mebe844d]
第一種の過誤の検討に専念して、以下の手順で帰無仮説を棄却...
+ 帰無仮説（今回は一定モデル）が正しいものと仮定する
+ 観測データに帰無モデルをあてはめ、\(\hat{\beta_1}=2.058...
+ 上記の真のモデルからサンプルデータを生成し、その度に一...
+ 求まった分布から一定モデル（帰無モデル）とxモデルの逸脱...
*** サンプルの生成 [#s77e1737]
Sageには、ポアソン分布を生成する関数が無いので、Rのrpois...
sample_size個生成する関数genSampleを以下の様に定義します。
sageへの入力:
#pre{{
# sageobjを使って変換する方が速いのでgenSample関数を使用...
def genSample(lambda_target, sample_size):
return np.array(sageobj(r('rpois(%d, %f)'%(sample_siz...
}}
*** パラメトリックブートストラップ法 [#w46fc378]
パラメトリックブートストラップ法をSageで試してみます。4章...
一定モデルを真のモデルとしてサンプルデータを生成し、一定...
回帰分析を行い、対数尤度llfから逸脱度Dの差を求めます。
これをn_bootstrap回繰り返し、その分布を求めます。
sageへの入力:
#pre{{
# パラメトリックブートストラップ法
def pb(d, n_bootstrap):
N = len(d)
y_mean = d.y.mean()
def sampling(d, y_mean, N):
# 帰無仮説を真のモデルとしてサンプルデータを生成
d['rnd'] = genSample(y_mean, N)
fit1 = smf.glm('rnd ~ 1', data=d, family=sm.famil...
fit2 = smf.glm('rnd ~ x', data=d, family=sm.famil...
return -2*(fit1.llf - fit2.llf)
return [ sampling(d, y_mean, N) for i in range(n_boot...
}}
sageへの入力:
#pre{{
# 逸脱度の差を1000サンプル生成
dd12 = pb(d, 1000)
}}
*** サンプルの結果 [#b0f88167]
ブートストラップ法で求まった結果にdescribe（Rのsummaryに...
ggplotのヒストグラムでその分布を表示してみます。
sageへの入力:
#pre{{
dd12Df = pd.DataFrame({'x' : dd12})
dd12Df.describe()
}}
#pre{{
x
count 1.000000e+03
mean 1.016127e+00
std 1.438180e+00
min 2.531954e-08
25% 9.873545e-02
50% 4.598125e-01
75% 1.367469e+00
max 1.048328e+01
[8 rows x 1 columns]
}}
sageへの入力:
#pre{{
ggplot(dd12Df, aes(x="x")) + geom_histogram(color="grey")...
}}
#pre{{
binwidth defaulted to range/30. Use 'binwidth = x' to adj...
<ggplot: (13280017)>
}}
sageへの入力:
#pre{{
ggsave('fig-5.4.png', dpi=50)
}}
#pre{{
Saving 11.0 x 8.0 in image.
}}
&ref(fig-5.4.png);
*** 仮説の検証 [#y8954cad]
4.5以上のデータの個数は、37個でp値は、37/1000=0.037となり...
3.90と求まりました。
これで、有意水準5%の検定において帰無仮説（一定モデル）は...
sageへの入力:
#pre{{
# 4.5以上のデータの個数は
len(dd12Df[dd12Df.x >= 4.5])
}}
#pre{{
37
}}
sageへの入力:
#pre{{
# 95%となるxの値を求める
dd12Df.quantile(0.95)
}}
#pre{{
x 3.899619
dtype: float64
}}
** χ自乗分布を使かう方法 [#r56bd0ba]
残念ながら、statsmodelsにはglmに対するanovaがないため、R...
p値は、0.034となり、同様に帰無仮説は棄却されました。（Rを...
sageへの入力:
#pre{{
# χ自乗分布を使かう方法
fit1 = smf.glm('y ~ 1', data=d, family=sm.families.Poisso...
fit2 = smf.glm('y ~ x', data=d, family=sm.families.Poisso...
}}
sageへの入力:
#pre{{
# 残念ながらstatsmodelsにはglmに対するanovaがないみたい
from statsmodels.stats.anova import anova_lm
anova_lm(fit1, fit2)
}}
#pre{{
Traceback (most recent call last):
File "<stdin>", line 1, in <module>
File "_sage_input_80.py", line 10, in <module>
exec compile(u'open("___code___.py","w").write("# -*-...
途中省略
35, in __getattribute__
obj = getattr(results, attr)
AttributeError: 'GLMResults' object has no attribute 'ssr'
}}
sageへの入力:
#pre{{
# Rにｄを渡して計算
PandaDf2RDf(d, "d")
}}
sageへの入力:
#pre{{
# ANOVA(ANalysis Of VAriance)を使って逸脱度を調べる
r('fit1 <- glm(y ~ 1, data=d, family=poisson)')
r('fit2 <- glm(y ~ x, data=d, family=poisson)')
r('anova(fit1, fit2, test = "Chisq")')
}}
#pre{{
Analysis of Deviance Table
Model 1: y ~ 1
Model 2: y ~ x
Resid. Df Resid. Dev Df Deviance Pr(>Chi)
1 99 89.507
2 98 84.993 1 4.5139 0.03362 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ...
}}
** コメント [#n4fa83ed]
#vote(おもしろかった[1],そうでもない[0],わかりずらい[0])
皆様のご意見、ご希望をお待ちしております。
#comment_kcaptcha
ページ名:
SmartDoc
![[PukiWiki]](image/pukiwiki.png "[PukiWiki]")
