kubo_book_chap4

3145 days ago by takepwave

Hiroshi TAKEMOTO (take@pwv.co.jp)

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

この企画は、雑誌や教科書にでているグラフをSageで再現し、 グラフの意味を理解すると共にSageの使い方をマスターすることを目的としています。

前回に続き、データ解析のための統計モデリング入門 (以下、久保本と書きます) の第4章の例題をSageを使って再現してみます。

4章のメインは、久保本図4.9(以下に引用)のバイアス補正とその分布にあると思います。

数式処理システムSageのノートブックは、計算結果を表示するだけではなく、実際に動かすことができるの大きな特徴です。 この機会にSageを分析に活用してみてはいかがでしょう。

前準備

最初に必要なライブラリーやパッケージをロードしておきます。

# 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') 
       

あてはまりの良さとモデルの複雑さ

「あてはまりの良さとモデルの複雑さ」を説明するために、3章のデータを使ってk=1とk=7の結果の比較が示されています。 このような図を簡単に再現できれば、実際の分析の時にも役立つと思います。

3章のデータをも一度読み込み、statsmodelsのglmを使ってk=1とその結果をプロットします。

# 3章のデータをもう一度読み込む d = pd.read_csv('http://hosho.ees.hokudai.ac.jp/~kubo/stat/iwanamibook/fig/poisson/data3a.csv') 
       
# statsmodelsを使ってglmを計算します import statsmodels.formula.api as smf import statsmodels.api as sm from scipy.stats.stats import pearsonr 
       
# k=1でのglm回帰を実行 fit_1 = smf.glm('y ~ 1', data=d, family=sm.families.Poisson(link=sm.families.links.log)).fit() # λの予測値をλ_1にセット d['lam_1'] = fit_1.predict() 
       
# λの予測値と一緒にプロット ggplot(d, aes(x='x', y='y')) + geom_point() + geom_line(aes(x='x', y='lam_1')) 
       
<ggplot: (8752412588061)>
<ggplot: (8752412588061)>
GgSave('fig-4.1-a.png') 
       
Saving 11.0 x 8.0 in image.
Saving 11.0 x 8.0 in image.

多次元(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)の状態になっています。過学習では、学習用データ(ここでいうサンプルデータ)に引きずられ、 真のモデルから離れてしまいます。

# k=7のglm計算がpythonではできないため、Rで計算 PandaDf2RDf(d, 'd') graph = preGraph("fig-4.1b-R.pdf") r("p <- ggplot(d, aes(x=x, y=y)) + geom_point() + geom_smooth(method=glm, family=poisson, formula =y ~ poly(x,7)) ") r('plot(p)') postGraph(graph) 
       

さまざまな逸脱度

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

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

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

# 様々な逸脱度 hdr = ['名前', '定義'] tbl = [ ['逸脱度(D)', '-2 log L*'], ['最小の逸脱度', 'フルモデルをあてはめたときのD'], ['残差逸脱度', 'D-最小のD'], ['最大の逸脱度', 'NullモデルをあてはめたときのD'], ['Null逸脱度', '最大のD - 最小のD']] html.table(tbl, hdr) 
       
名前 定義
逸脱度(D) -2 log L*
最小の逸脱度 フルモデルをあてはめたときのD
残差逸脱度 D-最小のD
最大の逸脱度 NullモデルをあてはめたときのD
Null逸脱度 最大のD - 最小のD
名前 定義
逸脱度(D) -2 log L*
最小の逸脱度 フルモデルをあてはめたときのD
残差逸脱度 D-最小のD
最大の逸脱度 NullモデルをあてはめたときのD
Null逸脱度 最大のD - 最小のD

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

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

# 逸脱度(Deviance)は、-2 log L* と定義されており、glmの結果の-2*llfから計算できる # 各モデルに対するfitを計算してDevianceの違いを整理してみる fit_f = smf.glm('y ~ f', data=d, family=sm.families.Poisson(link=sm.families.links.log)).fit() fit_x = smf.glm('y ~ x', data=d, family=sm.families.Poisson(link=sm.families.links.log)).fit() fit_x_f = smf.glm('y ~ x + f', data=d, family=sm.families.Poisson(link=sm.families.links.log)).fit() 
       

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

# フルモデルの対数尤度の計算 fullLogL = sageobj(sum([r.dpois(y, y, log=True) for y in d.y])); fullLogL 
       
-192.889752524496
-192.889752524496

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

# 結果を小数点3桁で表示する def prf3(f): return '%.3f'%f 
       
# フルモデルは、k=Nで結果がすべて観測値と一致します。従って残差は0となります。 # log L*: fit.llf, residual deviance: fit.deviance, AIC: fit.aicで取得できる hdr = ['モデル', 'k', "log L*", 'deviance -2 log L*', 'residual deviance', 'AIC'] tbl = [['一定', '1', prf3(fit_1.llf), prf3(-2*fit_1.llf), prf3(fit_1.deviance), prf3(fit_1.aic)], ['f', '2', prf3(fit_f.llf), prf3(-2*fit_f.llf), prf3(fit_f.deviance), prf3(fit_f.aic)], ['x', '2', prf3(fit_x.llf), prf3(-2*fit_x.llf), prf3(fit_x.deviance), prf3(fit_x.aic)], ['x+f', '3', prf3(fit_x_f.llf), prf3(-2*fit_x_f.llf), prf3(fit_x_f.deviance), prf3(fit_x_f.aic)], ['フル', '100', prf3(fullLogL), prf3(-2*fullLogL), '0', prf3(-2*fullLogL+2*100)]] html.table(tbl, hdr) 
       
モデル k log L* deviance -2 log L* residual deviance AIC
一定 1 -237.643 475.286 89.507 477.286
f 2 -237.627 475.255 89.475 479.255
x 2 -235.386 470.773 84.993 474.773
x+f 3 -235.294 470.587 84.808 476.587
フル 100 -192.890 385.780 0 585.780
モデル k log L* deviance -2 log L* residual deviance AIC
一定 1 -237.643 475.286 89.507 477.286
f 2 -237.627 475.255 89.475 479.255
x 2 -235.386 470.773 84.993 474.773
x+f 3 -235.294 470.587 84.808 476.587
フル 100 -192.890 385.780 0 585.780

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

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

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

観測データが一つの場合

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

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

# サンプル生成λ=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_target))) 
       
# 真のλ=8から50個のポアソン分布のサンプルを生成する N = 50 lambda_true = 8 sample = genSample(lambda_true, N) sampleDf = pd.DataFrame({'y' : np.array(sample)}) 
       
# 対数尤度の計算 def logL(sample, lambda_estimated): return sageobj(sum([r.dpois(y, lambda_estimated, log=True) for y in sample])) 
       

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

# k=1でのglm回帰を実行 fit = smf.glm('y ~ 1', data=sampleDf, family=sm.families.Poisson(link=sm.families.links.log)).fit() 
       
# fit.summary2() # 最大対数尤度 mxLogL = fit.llf # 切片 beta_1 = fit.params[0] print mxLogL, beta_1 
       
-125.051229132 2.1377104498
-125.051229132 2.1377104498

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

# β= [1.95, 2.20]の範囲の対数尤度を求める logL_lst = [ (x,logL(sample, exp(x)) )for x in np.arange(1.95, 2.20, 0.01)] # log L*の曲線プロット logL_plt = list_plot(logL_lst, figsize=5, plotjoined =True) 
       

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

# 真のλ=8から生成されたサンプル200セットにβの推定値を使って対数尤度を計算 # logL200 = [ logL(genSample(lambda_true, N), exp(beta_1)) for i in range(50)] # とても遅いので、 logL200 = [ sageobj(r('sum(log(dpois(rpois(%d, %f), %f)))'%(N, lambda_true, exp(beta_1)))) for i in range(200)] 
       

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

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, color="green") (logL_plt + sample_plt + logL_pt + logL_mean_pt).show() 
       
# 実際の解析では不可能だが、今回は真のモデルが分かっているので、評価用のデータを生成し、E(logL)を計算できる 
       

繰り返しは簡単

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

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

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

# 上記の計算を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.families.Poisson(link=sm.families.links.log)).fit() # 最大対数尤度 mxLogL = fit.llf # 切片 beta_1 = fit.params[0] logL200 = [ sageobj(r('sum(log(dpois(rpois(%d, %f), %f)))'%(N, lambda_true, exp(beta_1)))) for i in range(200)] # プロット 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="green") # 結果を表示 g.show(figsize=5) 
       

バイアスbの計算

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

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

# バイアス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.families.Poisson(link=sm.families.links.log)).fit() # 最大対数尤度 mxLogL = fit.llf # 切片 beta_1 = fit.params[0] logL200 = [ sageobj(r('sum(log(dpois(rpois(%d, %f), %f)))'%(N, lambda_true, exp(beta_1)))) for i in range(200)] return (mxLogL - sageobj(mean(logL200)) ) 
       
# 非常に時間がかかります bias_sample = [bias(lambda_true, N) for i in range(200)] 
       
mean(bias_sample) 
       
1.071531596453642
1.071531596453642
# biasのヒストグラムを表示する dd = pd.DataFrame({'y': np.array([ float(x) for x in bias_sample])}) ggplot(dd, aes(x='y')) + geom_histogram() 
       
stat_bin: binwidth defaulted to range/30.
    Use 'binwidth = x' to adjust this.
/usr/local/sage-6.7/local/lib/python2.7/site-packages/pandas-0.15.2-py2.\
7-linux-x86_64.egg/pandas/util/decorators.py:81: FutureWarning: the
'rows' keyword is deprecated, use 'index' instead
  warnings.warn(msg, FutureWarning)
<ggplot: (8752413086621)>
stat_bin: binwidth defaulted to range/30.
    Use 'binwidth = x' to adjust this.
/usr/local/sage-6.7/local/lib/python2.7/site-packages/pandas-0.15.2-py2.7-linux-x86_64.egg/pandas/util/decorators.py:81: FutureWarning: the 'rows' keyword is deprecated, use 'index' instead
  warnings.warn(msg, FutureWarning)
<ggplot: (8752413086621)>
ggsave('fig-4.10.png', dpi=50) 
       
Saving 11.0 x 8.0 in image.
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=1という値は、かなりぶれることも分かりました。

# 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
$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
# サンプリングが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
[1] 0.9048361
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)