sage/データ解析のための統計モデリング入門勉強ノート第4章
をテンプレートにして作成
[
トップ
] [
新規
|
一覧
|
単語検索
|
最終更新
|
ヘルプ
]
開始行:
[[FrontPage]]
#contents
2014/03/30からのアクセス回数 &counter;
ここで紹介したSageワークシートは、以下のURLからダウンロー...
http://www15191ue.sakura.ne.jp:8000/home/pub/41/
また、Sageのサーバを公開しているサイト（http://www.sagenb...
アップロードし、実行したり、変更していろいろ動きを試すこ...
* Sageでグラフを再現してみよう：データ解析のための統計モ...
この企画は、雑誌や教科書にでているグラフをSageで再現し、
グラフの意味を理解すると共にSageの使い方をマスターするこ...
前回に続き、
[[データ解析のための統計モデリング入門>http://www.amazon....
（以下、久保本と書きます）
の第４章の例題をSageを使って再現してみます。
4章のメインは、図4.9のバイアス補正とその分布にあると思い...
&ref(fig-4.9.png);
数式処理システムSageのノートブックは、計算結果を表示する...
この機会にSageを分析に活用してみてはいかがでしょう。
** 前準備 [#o17d4d36]
最初に必要なライブラリーやパッケージをロードしておきます。
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')
}}
** あてはまりの良さとモデルの複雑さ [#y5237818]
「あてはまりの良さとモデルの複雑さ」を説明するために、３...
このような図を簡単に再現できれば、実際の分析の時にも役立...
3章のデータをも一度読み込み、statsmodelsのglmを使ってk=1...
sageへの入力:
#pre{{
# 3章のデータをもう一度読み込む
d = pd.read_csv('http://hosho.ees.hokudai.ac.jp/~kubo/sta...
}}
sageへの入力:
#pre{{
# statsmodelsを使ってglmを計算します
import statsmodels.formula.api as smf
import statsmodels.api as sm
from scipy.stats.stats import pearsonr
}}
sageへの入力:
#pre{{
# k=1でのglm回帰を実行
fit_1 = smf.glm('y ~ 1', data=d, family=sm.families.Poiss...
# λの予測値をλ_1にセット
d['lam_1'] = fit_1.predict()
}}
sageへの入力:
#pre{{
# λの予測値と一緒にプロット
ggplot(d, aes(x='x', y='y')) + geom_point() + geom_line(a...
}}
#pre{{
<ggplot: (39269997)>
}}
sageへの入力:
#pre{{
ggsave('fig-4.1-a.png', dpi=50)
}}
#pre{{
Saving 11.0 x 8.0 in image.
}}
&ref(fig-4.1-a.png);
*** 多次元（poly）のglmは、statsmodelsではできない [#qf23...
次にk=7（\(log \lambda = \beta_1 + \beta_2 x + ... + \bet...
残念ながら、statsmodelsは、線形項の多次元のglmをサポート...
Rのggplot2には、geom_smoothというlmやglmの計算結果を使っ...
これを使ってk=7の結果を表示してみましょう。線形項の多項式...
Sage（Pandas）とRのデータフレームを相互に変換したり、Rの...
ともて簡単にk=7の結果を図化できます。
k=1とk=7の結果をみるとk=7の結果は、あまりにも複雑な曲線を...
（Overfitting）の状態になっています。過学習では、学習用デ...
真のモデルから離れてしまいます。
sageへの入力:
#pre{{
# k=7のglm計算がpythonではできないため、Rで計算
PandaDf2RDf(d, 'd')
graph = preGraph("fig-4.1b-R.png")
r("p <- ggplot(d, aes(x=x, y=y)) + geom_point() + geom_sm...
r('plot(p)')
postGraph(graph)
}}
&ref(fig-4.1b-R.png);
** さまざまな逸脱度 [#w7564e95]
逸脱度は、「あてはまりの悪さ」を表現する指標と久保本では...
最大対数尤度を\(log L^*\)とすると、逸脱度(D)は以下の様に...
$$
D = -2 log L^*
$$
4章で登場する逸脱度を以下の表に示します。
最大の逸脱度は、一定モデルで、最小の逸脱度はフルモデルで...
sageへの入力:
#pre{{
# 様々な逸脱度
hdr = ['名前', '定義']
tbl = [ ['逸脱度(D)', '-2 log L*'],
['最小の逸脱度', 'フルモデルをあてはめたとき...
['残差逸脱度', 'D-最小のD'],
['最大の逸脱度', 'Nullモデルをあてはめたとき...
['Null逸脱度', '最大のD - 最小のD']]
html.table(tbl, hdr)
}}
&ref(tbl-1.png);
各モデルのglmのfit結果を求めて、その結果を表にまとめます。
ここで出てくるモデル選択の基準となるAIC(Akaike's informat...
$$
AIC = -2 { (最大対数尤度) - (最尤推定したパラメータの数)}...
$$
sageへの入力:
#pre{{
# 逸脱度（Deviance）は、-2 log L* と定義されており、glmの...
# 各モデルに対するfitを計算してDevianceの違いを整理してみる
fit_f = smf.glm('y ~ f', data=d, family=sm.families.Poiss...
fit_x = smf.glm('y ~ x', data=d, family=sm.families.Poiss...
fit_x_f = smf.glm('y ~ x + f', data=d, family=sm.families...
}}
データの個数と同じ数のパラメータを使ってあてはめたモデル...
これは、実際にそのようなモデルを作って計算するのでは無く...
その結果、残差は０であり、対数尤度は、以下のようになりま...
sageへの入力:
#pre{{
# フルモデルの対数尤度の計算
fullLogL = sum([r.dpois(y, y, log=True) for y in d.y]).n(...
}}
#pre{{
-192.889752524496
}}
結果をみやすくするために、小数点以下3桁を表示する関数prf3...
sageへの入力:
#pre{{
# 結果を小数点3桁で表示する
def prf3(f):
return '%.3f'%f
}}
sageへの入力:
#pre{{
# フルモデルは、k=Nで結果がすべて観測値と一致します。従っ...
# log L*: fit.llf, residual deviance: fit.deviance, AIC: ...
hdr = ['モデル', 'k', "log L*", 'deviance -2 log L*', 're...
tbl = [['一定', '1', prf3(fit_1.llf), prf3(-2*fit_1.llf)...
['f', '2', prf3(fit_f.llf), prf3(-2*fit_f.ll...
['x', '2', prf3(fit_x.llf), prf3(-2*fit_x.ll...
['x+f', '3', prf3(fit_x_f.llf), prf3(-2*fit_...
['フル', '100', prf3(fullLogL), prf3(-2*fullL...
html.table(tbl, hdr)
}}
&ref(tbl-2.png);
** 平均対数尤度を使ったモデルの予測の良さ [#l57acb92]
いよいよ図４．９の計算に入ります。久保本ではモデル選択の...
数値例で示しています。（これが久保本のすごいところ！）
平均の対数尤度は、真のモデルを使って生成されたサンプルか...
このようなことは不可能ですが、平均対数尤度と推定されたモ...
表現しています。bの定義は以下の通りです。
$$
b = log L^* - E(Log L)
$$
*** 観測データが一つの場合 [#v1566cfa]
図４．９の(A)観測データが一つの場合について、真のモデルλ=...
準備としてサンプル生成関数genSample, 対数尤度を計算するlo...
sageへの入力:
#pre{{
# サンプル生成λ=8ポアソン分布のサンプルをN=50個生成する
# sample = [ float(y.n()) for y in r.rpois(N, lambda_true)]
# sageobjを使って変換する方が速いのでgenSample関数を使用...
def genSample(lambda_target, sample_size):
return sageobj(r('rpois(%d, %f)'%(sample_size, lambda...
}}
sageへの入力:
#pre{{
# 真のλ=8から50個のポアソン分布のサンプルを生成する
N = 50
lambda_true = 8
sample = genSample(lambda_true, N)
sampleDf = pd.DataFrame({'y' : np.array(sample)})
}}
sageへの入力:
#pre{{
# 対数尤度の計算
def logL(sample, lambda_estimated):
return sum([r.dpois(y, lambda_estimated, log=True) fo...
}}
glmを使って一定モデルの結果を取得します。ここでは対数尤度...
sageへの入力:
#pre{{
# k=1でのglm回帰を実行
fit = smf.glm('y ~ 1', data=sampleDf, family=sm.families....
}}
sageへの入力:
#pre{{
# fit.summary2()
# 最大対数尤度
mxLogL = fit.llf
# 切片
beta_1 = fit.params[0]
print mxLogL, beta_1
}}
#pre{{
-119.642879665 2.00148000021
}}
この結果を図化し、logL_pltにセットします。
sageへの入力:
#pre{{
# β= [1.95, 2.20]の範囲の対数尤度を求める
logL_lst = [ (x,logL(sample, exp(x)) )for x in np.arange...
# log L*の曲線プロット
logL_plt = list_plot(logL_lst, figsize=5, plotjoined =True)
}}
つぎに、真のモデルから50個のサンプルを200セット生成し、平...
sageへの入力:
#pre{{
# 真のλ=8から生成されたサンプル200セットにβの推定値を使っ...
# logL200 = [ logL(genSample(lambda_true, N), exp(beta_1)...
# とても遅いので、
logL200 = [ r('sum(log(dpois(rpois(%d, %f), %f)))'%(N, la...
}}
求まった結果を一つにプロットします。logL*を赤○で、平均対...
各βの対数尤度曲線と一緒に表示したのが、以下の図です。
sageへの入力:
#pre{{
sample_plt = list_plot([ (beta_1, y) for y in logL200])
logL_pt = point((beta_1, mxLogL), size=40, color="red")
logL_mean_pt = point((beta_1, mean(logL200)), size=40, co...
(logL_plt + sample_plt + logL_pt + logL_mean_pt).show()
}}
&ref(sage0.png);
sageへの入力:
#pre{{
# 実際の解析では不可能だが、今回は真のモデルが分かってい...
}}
*** 繰り返しは簡単 [#c4679c67]
次に上記の処理を繰り返した図４．９(B)を試してみます。
Sageでは、一度試した処理をコピー＆ペーストしてループに入...
処理を進めることで効率良く目的のグラフを得ることができま...
最初にからのGraphics=gを生成し、上記の処理と同じことをし...
たったこれだけで、図４．９の(B)が再現できます。
久保本では、平均対数尤度が結構きれいな曲線で表示されてい...
sageへの入力:
#pre{{
# 上記の計算を25回繰り返しプロット図4.9(B)
g = Graphics()
for i in range(25):
sample = genSample(lambda_true, N)
sampleDf = pd.DataFrame({'y' : np.array(sample)})
# k=1でのglm回帰を実行
fit = smf.glm('y ~ 1', data=sampleDf, family=sm.famil...
mxLogL = fit.llf
# 切片
beta_1 = fit.params[0]
logL200 = [ r('sum(log(dpois(rpois(%d, %f), %f)))'%(N...
# プロット
g += list_plot([ (beta_1, y) for y in logL200])
g += point((beta_1, mxLogL), size=40, color="red")
g += point((beta_1, mean(logL200)), size=40, color="g...
# 結果を表示
g.show(figsize=5)
}}
&ref(sage1.png);
** バイアスbの計算 [#vfcdcf60]
時間がかかったので、２００セットのデータに対してバイアスb...
バイアスbは、0.7と久保本の1.01とはずれています。
sageへの入力:
#pre{{
# バイアスbの計算
def bias(lambda_true, N):
sample = genSample(lambda_true, N)
sampleDf = pd.DataFrame({'y' : np.array(sample)})
# k=1でのglm回帰を実行
fit = smf.glm('y ~ 1', data=sampleDf, family=sm.famil...
# 最大対数尤度
mxLogL = fit.llf
# 切片
beta_1 = fit.params[0]
logL200 = [ r('sum(log(dpois(rpois(%d, %f), %f)))'%(N...
return (mxLogL - mean(logL200).n() )
}}
sageへの入力:
#pre{{
# 非常に時間がかかります
bias_sample = [bias(lambda_true, N) for i in range(200)]
}}
sageへの入力:
#pre{{
mean(bias_sample)
}}
#pre{{
(0.69534825087697705+0j)
}}
sageへの入力:
#pre{{
# biasのヒストグラムを表示する
dd = pd.DataFrame({'y': np.array([ float(x) for x in bias...
ggplot(dd, aes(x='y')) + geom_histogram()
}}
#pre{{
binwidth defaulted to range/30. Use 'binwidth = x' to adj...
<ggplot: (39926417)>
}}
sageへの入力:
#pre{{
ggsave('fig-4.10.png', dpi=50)
}}
#pre{{
Saving 11.0 x 8.0 in image.
}}
&ref(fig-4.10.png);
*** Rで計算 [#nef39ffe]
久保本のサポートページにあるRの関数を使ってサンプル1000セ...
これでも 1.213と計算され、その分布をみると-20から20の区...
バイアスの定義を変形すると、平均対数尤度\(E(logL)\)は以下...
$$
E(log L) = log L^* - b
$$
ここで、一定モデルのk=1であることからバイアスbとkを入れ替...
AICが予測の「悪さ」を表す指標となっていることがうなずけま...
$$
AIC = -2 \times (log L^* -1)
$$
また、上記の数値実験の結果得られたバイアスb=１という値は...
sageへの入力:
#pre{{
# Rでbiasを計算する関数を定義すると非常に速く計算できる
r('source("%s") ' % (DATA + 'bias.txt'))
}}
#pre{{
$value
function (lambda.true, sample.size)
{
sample.rpois <- rpois(sample.size, lambda.true)
fit <- glm(sample.rpois ~ 1, family = poisson)
lambda.estimated <- exp(coef(fit))
likelihood.mean <- mean(sapply(1:200, function(i) {
sum(log(dpois(rpois(N, lambda.true), lambda.estim...
}))
logLik(fit) - likelihood.mean
}
$visible
[1] FALSE
}}
sageへの入力:
#pre{{
# サンプリングが1000くらいでもbiasの平均は、結構ぶれる
r('N <- 50')
r('lambda.true <- 8')
r('bias.sampled <- sapply(1:1000, function(i) bias(lambda...
r('mean(bias.sampled)')
}}
#pre{{
[1] 1.213076
}}
sageへの入力:
#pre{{
graph = preGraph("fig-4.10-R.png")
r('p <- qplot(bias.sampled, geom = "blank")
+ geom_histogram(aes(y = ..density..), colour = "bl...
+ geom_density(alpha = 0.2, fill = "#6666FF")
+ geom_vline(aes(xintercept = mean(bias.sampled)), c...
r('plot(p)')
postGraph(graph)
}}
&ref(fig-4.10-R.png);
** コメント [#zc1b7f5e]
#vote(おもしろかった[2],そうでもない[0],わかりずらい[0])
皆様のご意見、ご希望をお待ちしております。
#comment_kcaptcha
終了行:
[[FrontPage]]
#contents
2014/03/30からのアクセス回数 &counter;
ここで紹介したSageワークシートは、以下のURLからダウンロー...
http://www15191ue.sakura.ne.jp:8000/home/pub/41/
また、Sageのサーバを公開しているサイト（http://www.sagenb...
アップロードし、実行したり、変更していろいろ動きを試すこ...
* Sageでグラフを再現してみよう：データ解析のための統計モ...
この企画は、雑誌や教科書にでているグラフをSageで再現し、
グラフの意味を理解すると共にSageの使い方をマスターするこ...
前回に続き、
[[データ解析のための統計モデリング入門>http://www.amazon....
（以下、久保本と書きます）
の第４章の例題をSageを使って再現してみます。
4章のメインは、図4.9のバイアス補正とその分布にあると思い...
&ref(fig-4.9.png);
数式処理システムSageのノートブックは、計算結果を表示する...
この機会にSageを分析に活用してみてはいかがでしょう。
** 前準備 [#o17d4d36]
最初に必要なライブラリーやパッケージをロードしておきます。
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')
}}
** あてはまりの良さとモデルの複雑さ [#y5237818]
「あてはまりの良さとモデルの複雑さ」を説明するために、３...
このような図を簡単に再現できれば、実際の分析の時にも役立...
3章のデータをも一度読み込み、statsmodelsのglmを使ってk=1...
sageへの入力:
#pre{{
# 3章のデータをもう一度読み込む
d = pd.read_csv('http://hosho.ees.hokudai.ac.jp/~kubo/sta...
}}
sageへの入力:
#pre{{
# statsmodelsを使ってglmを計算します
import statsmodels.formula.api as smf
import statsmodels.api as sm
from scipy.stats.stats import pearsonr
}}
sageへの入力:
#pre{{
# k=1でのglm回帰を実行
fit_1 = smf.glm('y ~ 1', data=d, family=sm.families.Poiss...
# λの予測値をλ_1にセット
d['lam_1'] = fit_1.predict()
}}
sageへの入力:
#pre{{
# λの予測値と一緒にプロット
ggplot(d, aes(x='x', y='y')) + geom_point() + geom_line(a...
}}
#pre{{
<ggplot: (39269997)>
}}
sageへの入力:
#pre{{
ggsave('fig-4.1-a.png', dpi=50)
}}
#pre{{
Saving 11.0 x 8.0 in image.
}}
&ref(fig-4.1-a.png);
*** 多次元（poly）のglmは、statsmodelsではできない [#qf23...
次にk=7（\(log \lambda = \beta_1 + \beta_2 x + ... + \bet...
残念ながら、statsmodelsは、線形項の多次元のglmをサポート...
Rのggplot2には、geom_smoothというlmやglmの計算結果を使っ...
これを使ってk=7の結果を表示してみましょう。線形項の多項式...
Sage（Pandas）とRのデータフレームを相互に変換したり、Rの...
ともて簡単にk=7の結果を図化できます。
k=1とk=7の結果をみるとk=7の結果は、あまりにも複雑な曲線を...
（Overfitting）の状態になっています。過学習では、学習用デ...
真のモデルから離れてしまいます。
sageへの入力:
#pre{{
# k=7のglm計算がpythonではできないため、Rで計算
PandaDf2RDf(d, 'd')
graph = preGraph("fig-4.1b-R.png")
r("p <- ggplot(d, aes(x=x, y=y)) + geom_point() + geom_sm...
r('plot(p)')
postGraph(graph)
}}
&ref(fig-4.1b-R.png);
** さまざまな逸脱度 [#w7564e95]
逸脱度は、「あてはまりの悪さ」を表現する指標と久保本では...
最大対数尤度を\(log L^*\)とすると、逸脱度(D)は以下の様に...
$$
D = -2 log L^*
$$
4章で登場する逸脱度を以下の表に示します。
最大の逸脱度は、一定モデルで、最小の逸脱度はフルモデルで...
sageへの入力:
#pre{{
# 様々な逸脱度
hdr = ['名前', '定義']
tbl = [ ['逸脱度(D)', '-2 log L*'],
['最小の逸脱度', 'フルモデルをあてはめたとき...
['残差逸脱度', 'D-最小のD'],
['最大の逸脱度', 'Nullモデルをあてはめたとき...
['Null逸脱度', '最大のD - 最小のD']]
html.table(tbl, hdr)
}}
&ref(tbl-1.png);
各モデルのglmのfit結果を求めて、その結果を表にまとめます。
ここで出てくるモデル選択の基準となるAIC(Akaike's informat...
$$
AIC = -2 { (最大対数尤度) - (最尤推定したパラメータの数)}...
$$
sageへの入力:
#pre{{
# 逸脱度（Deviance）は、-2 log L* と定義されており、glmの...
# 各モデルに対するfitを計算してDevianceの違いを整理してみる
fit_f = smf.glm('y ~ f', data=d, family=sm.families.Poiss...
fit_x = smf.glm('y ~ x', data=d, family=sm.families.Poiss...
fit_x_f = smf.glm('y ~ x + f', data=d, family=sm.families...
}}
データの個数と同じ数のパラメータを使ってあてはめたモデル...
これは、実際にそのようなモデルを作って計算するのでは無く...
その結果、残差は０であり、対数尤度は、以下のようになりま...
sageへの入力:
#pre{{
# フルモデルの対数尤度の計算
fullLogL = sum([r.dpois(y, y, log=True) for y in d.y]).n(...
}}
#pre{{
-192.889752524496
}}
結果をみやすくするために、小数点以下3桁を表示する関数prf3...
sageへの入力:
#pre{{
# 結果を小数点3桁で表示する
def prf3(f):
return '%.3f'%f
}}
sageへの入力:
#pre{{
# フルモデルは、k=Nで結果がすべて観測値と一致します。従っ...
# log L*: fit.llf, residual deviance: fit.deviance, AIC: ...
hdr = ['モデル', 'k', "log L*", 'deviance -2 log L*', 're...
tbl = [['一定', '1', prf3(fit_1.llf), prf3(-2*fit_1.llf)...
['f', '2', prf3(fit_f.llf), prf3(-2*fit_f.ll...
['x', '2', prf3(fit_x.llf), prf3(-2*fit_x.ll...
['x+f', '3', prf3(fit_x_f.llf), prf3(-2*fit_...
['フル', '100', prf3(fullLogL), prf3(-2*fullL...
html.table(tbl, hdr)
}}
&ref(tbl-2.png);
** 平均対数尤度を使ったモデルの予測の良さ [#l57acb92]
いよいよ図４．９の計算に入ります。久保本ではモデル選択の...
数値例で示しています。（これが久保本のすごいところ！）
平均の対数尤度は、真のモデルを使って生成されたサンプルか...
このようなことは不可能ですが、平均対数尤度と推定されたモ...
表現しています。bの定義は以下の通りです。
$$
b = log L^* - E(Log L)
$$
*** 観測データが一つの場合 [#v1566cfa]
図４．９の(A)観測データが一つの場合について、真のモデルλ=...
準備としてサンプル生成関数genSample, 対数尤度を計算するlo...
sageへの入力:
#pre{{
# サンプル生成λ=8ポアソン分布のサンプルをN=50個生成する
# sample = [ float(y.n()) for y in r.rpois(N, lambda_true)]
# sageobjを使って変換する方が速いのでgenSample関数を使用...
def genSample(lambda_target, sample_size):
return sageobj(r('rpois(%d, %f)'%(sample_size, lambda...
}}
sageへの入力:
#pre{{
# 真のλ=8から50個のポアソン分布のサンプルを生成する
N = 50
lambda_true = 8
sample = genSample(lambda_true, N)
sampleDf = pd.DataFrame({'y' : np.array(sample)})
}}
sageへの入力:
#pre{{
# 対数尤度の計算
def logL(sample, lambda_estimated):
return sum([r.dpois(y, lambda_estimated, log=True) fo...
}}
glmを使って一定モデルの結果を取得します。ここでは対数尤度...
sageへの入力:
#pre{{
# k=1でのglm回帰を実行
fit = smf.glm('y ~ 1', data=sampleDf, family=sm.families....
}}
sageへの入力:
#pre{{
# fit.summary2()
# 最大対数尤度
mxLogL = fit.llf
# 切片
beta_1 = fit.params[0]
print mxLogL, beta_1
}}
#pre{{
-119.642879665 2.00148000021
}}
この結果を図化し、logL_pltにセットします。
sageへの入力:
#pre{{
# β= [1.95, 2.20]の範囲の対数尤度を求める
logL_lst = [ (x,logL(sample, exp(x)) )for x in np.arange...
# log L*の曲線プロット
logL_plt = list_plot(logL_lst, figsize=5, plotjoined =True)
}}
つぎに、真のモデルから50個のサンプルを200セット生成し、平...
sageへの入力:
#pre{{
# 真のλ=8から生成されたサンプル200セットにβの推定値を使っ...
# logL200 = [ logL(genSample(lambda_true, N), exp(beta_1)...
# とても遅いので、
logL200 = [ r('sum(log(dpois(rpois(%d, %f), %f)))'%(N, la...
}}
求まった結果を一つにプロットします。logL*を赤○で、平均対...
各βの対数尤度曲線と一緒に表示したのが、以下の図です。
sageへの入力:
#pre{{
sample_plt = list_plot([ (beta_1, y) for y in logL200])
logL_pt = point((beta_1, mxLogL), size=40, color="red")
logL_mean_pt = point((beta_1, mean(logL200)), size=40, co...
(logL_plt + sample_plt + logL_pt + logL_mean_pt).show()
}}
&ref(sage0.png);
sageへの入力:
#pre{{
# 実際の解析では不可能だが、今回は真のモデルが分かってい...
}}
*** 繰り返しは簡単 [#c4679c67]
次に上記の処理を繰り返した図４．９(B)を試してみます。
Sageでは、一度試した処理をコピー＆ペーストしてループに入...
処理を進めることで効率良く目的のグラフを得ることができま...
最初にからのGraphics=gを生成し、上記の処理と同じことをし...
たったこれだけで、図４．９の(B)が再現できます。
久保本では、平均対数尤度が結構きれいな曲線で表示されてい...
sageへの入力:
#pre{{
# 上記の計算を25回繰り返しプロット図4.9(B)
g = Graphics()
for i in range(25):
sample = genSample(lambda_true, N)
sampleDf = pd.DataFrame({'y' : np.array(sample)})
# k=1でのglm回帰を実行
fit = smf.glm('y ~ 1', data=sampleDf, family=sm.famil...
mxLogL = fit.llf
# 切片
beta_1 = fit.params[0]
logL200 = [ r('sum(log(dpois(rpois(%d, %f), %f)))'%(N...
# プロット
g += list_plot([ (beta_1, y) for y in logL200])
g += point((beta_1, mxLogL), size=40, color="red")
g += point((beta_1, mean(logL200)), size=40, color="g...
# 結果を表示
g.show(figsize=5)
}}
&ref(sage1.png);
** バイアスbの計算 [#vfcdcf60]
時間がかかったので、２００セットのデータに対してバイアスb...
バイアスbは、0.7と久保本の1.01とはずれています。
sageへの入力:
#pre{{
# バイアスbの計算
def bias(lambda_true, N):
sample = genSample(lambda_true, N)
sampleDf = pd.DataFrame({'y' : np.array(sample)})
# k=1でのglm回帰を実行
fit = smf.glm('y ~ 1', data=sampleDf, family=sm.famil...
# 最大対数尤度
mxLogL = fit.llf
# 切片
beta_1 = fit.params[0]
logL200 = [ r('sum(log(dpois(rpois(%d, %f), %f)))'%(N...
return (mxLogL - mean(logL200).n() )
}}
sageへの入力:
#pre{{
# 非常に時間がかかります
bias_sample = [bias(lambda_true, N) for i in range(200)]
}}
sageへの入力:
#pre{{
mean(bias_sample)
}}
#pre{{
(0.69534825087697705+0j)
}}
sageへの入力:
#pre{{
# biasのヒストグラムを表示する
dd = pd.DataFrame({'y': np.array([ float(x) for x in bias...
ggplot(dd, aes(x='y')) + geom_histogram()
}}
#pre{{
binwidth defaulted to range/30. Use 'binwidth = x' to adj...
<ggplot: (39926417)>
}}
sageへの入力:
#pre{{
ggsave('fig-4.10.png', dpi=50)
}}
#pre{{
Saving 11.0 x 8.0 in image.
}}
&ref(fig-4.10.png);
*** Rで計算 [#nef39ffe]
久保本のサポートページにあるRの関数を使ってサンプル1000セ...
これでも 1.213と計算され、その分布をみると-20から20の区...
バイアスの定義を変形すると、平均対数尤度\(E(logL)\)は以下...
$$
E(log L) = log L^* - b
$$
ここで、一定モデルのk=1であることからバイアスbとkを入れ替...
AICが予測の「悪さ」を表す指標となっていることがうなずけま...
$$
AIC = -2 \times (log L^* -1)
$$
また、上記の数値実験の結果得られたバイアスb=１という値は...
sageへの入力:
#pre{{
# Rでbiasを計算する関数を定義すると非常に速く計算できる
r('source("%s") ' % (DATA + 'bias.txt'))
}}
#pre{{
$value
function (lambda.true, sample.size)
{
sample.rpois <- rpois(sample.size, lambda.true)
fit <- glm(sample.rpois ~ 1, family = poisson)
lambda.estimated <- exp(coef(fit))
likelihood.mean <- mean(sapply(1:200, function(i) {
sum(log(dpois(rpois(N, lambda.true), lambda.estim...
}))
logLik(fit) - likelihood.mean
}
$visible
[1] FALSE
}}
sageへの入力:
#pre{{
# サンプリングが1000くらいでもbiasの平均は、結構ぶれる
r('N <- 50')
r('lambda.true <- 8')
r('bias.sampled <- sapply(1:1000, function(i) bias(lambda...
r('mean(bias.sampled)')
}}
#pre{{
[1] 1.213076
}}
sageへの入力:
#pre{{
graph = preGraph("fig-4.10-R.png")
r('p <- qplot(bias.sampled, geom = "blank")
+ geom_histogram(aes(y = ..density..), colour = "bl...
+ geom_density(alpha = 0.2, fill = "#6666FF")
+ geom_vline(aes(xintercept = mean(bias.sampled)), c...
r('plot(p)')
postGraph(graph)
}}
&ref(fig-4.10-R.png);
** コメント [#zc1b7f5e]
#vote(おもしろかった[2],そうでもない[0],わかりずらい[0])
皆様のご意見、ご希望をお待ちしております。
#comment_kcaptcha
ページ名:
SmartDoc
![[PukiWiki]](image/pukiwiki.png "[PukiWiki]")
