Sageでグラフを再現してみよう：データ解析のための統計モデリング入門第４章

多次元（poly）のglmは、statsmodelsではできない

次にk=7（$log \lambda = \beta_1 + \beta_2 x + ... + \beta_7 x^6$）の計算ですが、 残念ながら、statsmodelsは、線形項の多次元のglmをサポートしていないので、Rのglmを使って図化してみます。

Rのggplot2には、geom_smoothというlmやglmの計算結果を使ってスムーズ曲線を描画する機能があります。 これを使ってk=7の結果を表示してみましょう。線形項の多項式は、poly関数を使うと便利です。

Sage（Pandas）とRのデータフレームを相互に変換したり、Rの結果をSageで図化する関数をRUtilで定義してあるので、 ともて簡単にk=7の結果を図化できます。

k=1とk=7の結果をみるとk=7の結果は、あまりにも複雑な曲線を示しており、機械学習でいうところの「過学習」 （Overfitting）の状態になっています。過学習では、学習用データ（ここでいうサンプルデータ）に引きずられ、 真のモデルから離れてしまいます。

さまざまな逸脱度

逸脱度は、「あてはまりの悪さ」を表現する指標と久保本では説明されています。

最大対数尤度を$log L^*$とすると、逸脱度(D)は以下の様に定義されます。 $$ D = -2 log L^* $$

4章で登場する逸脱度を以下の表に示します。 最大の逸脱度は、一定モデルで、最小の逸脱度はフルモデルです。

各モデルのglmのfit結果を求めて、その結果を表にまとめます。

ここで出てくるモデル選択の基準となるAIC(Akaike's information criterion)は、以下の様に定義されています。 $$ AIC = -2 { (最大対数尤度) - (最尤推定したパラメータの数)} = D + 2k $$

データの個数と同じ数のパラメータを使ってあてはめたモデルがフルモデルです。 これは、実際にそのようなモデルを作って計算するのでは無く、フルモデルはすべての観測値が予測値と一致するモデルであり、 その結果、残差は０であり、対数尤度は、以下のようになります。

結果をみやすくするために、小数点以下3桁を表示する関数prf3を定義します。

平均対数尤度を使ったモデルの予測の良さ

いよいよ図４．９の計算に入ります。久保本ではモデル選択の指標としてAICが使えるのかを一定モデルを使って 数値例で示しています。（これが久保本のすごいところ！）

平均の対数尤度は、真のモデルを使って生成されたサンプルから求めた対数尤度の平均値であり、実際のモデルづくりでは このようなことは不可能ですが、平均対数尤度と推定されたモデルの関係とこの差（バイアス）bの分布からAICの意味するところを 表現しています。bの定義は以下の通りです。 $$ b = log L^* - E(Log L) $$

観測データが一つの場合

図４．９の(A)観測データが一つの場合について、真のモデルλ=8のポアソン分布から50個のサンプルを生成します。

準備としてサンプル生成関数genSample, 対数尤度を計算するlogLを定義します。

glmを使って一定モデルの結果を取得します。ここでは対数尤度が-119.6、 $\beta_1$の値は2.001と求まりました。

この結果を図化し、logL_pltにセットします。

つぎに、真のモデルから50個のサンプルを200セット生成し、平均の対数尤度を計算します。

求まった結果を一つにプロットします。logL*を赤○で、平均対数尤度を緑○で、各対数尤度を小さな青○で示し、 各βの対数尤度曲線と一緒に表示したのが、以下の図です。

繰り返しは簡単

次に上記の処理を繰り返した図４．９(B)を試してみます。 Sageでは、一度試した処理をコピー＆ペーストしてループに入れたり、関数にすることで結果を確認しながら 処理を進めることで効率良く目的のグラフを得ることができます。

最初に空のGraphics=gを生成し、上記の処理と同じことをして、プロット結果をgに加えていきます。 たったこれだけで、図４．９の(B)が再現できます。

久保本では、平均対数尤度が結構きれいな曲線で表示されていますが、私の計算ではぶれています。

バイアスbの計算

時間がかかったので、２００セットのデータに対してバイアスbを計算し、その分布を表示してみました。

バイアスbは、0.7と久保本の1.01とはずれています。

}}} {{{id=54| ggsave('fig-4.10.png', dpi=50) /// Saving 11.0 x 8.0 in image. }}}

Rで計算

久保本のサポートページにあるRの関数を使ってサンプル1000セットでバイアスbを計算してみました。 これでも　1.213と計算され、その分布をみると-20から20の区間でかなりぶれていることがわかります。

バイアスの定義を変形すると、平均対数尤度$E(logL)$は以下の様になります。 $$ E(log L) = log L^* - b $$ ここで、一定モデルのk=1であることからバイアスbとkを入れ替えるとAICの定義となり、 AICが予測の「悪さ」を表す指標となっていることがうなずけます。 $$ AIC = -2 \times (log L^* -1) $$

また、上記の数値実験の結果得られたバイアスb=１という値は、かなりぶれることも分かりました。

{{{id=56| # Rでbiasを計算する関数を定義すると非常に速く計算できる r('source("%s") ' % (DATA + 'bias.txt')) /// $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.estimated))) })) logLik(fit) - likelihood.mean } $visible [1] FALSE }}} {{{id=57| # サンプリングが1000くらいでもbiasの平均は、結構ぶれる r('N <- 50') r('lambda.true <- 8') r('bias.sampled <- sapply(1:1000, function(i) bias(lambda.true, N))') r('mean(bias.sampled)') /// [1] 0.9048361 }}} {{{id=58| graph = preGraph("fig-4.10-R.pdf") r('p <- qplot(bias.sampled, geom = "blank") + geom_histogram(aes(y = ..density..), colour = "black", fill = "white") + geom_density(alpha = 0.2, fill = "#6666FF") + geom_vline(aes(xintercept = mean(bias.sampled)), color = "red", linetype = "dashed", size = 2) ') r('plot(p)') postGraph(graph) /// }}} {{{id=59| /// }}}