sage/データ解析のための統計モデリング入門勉強ノート第4章

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

2014/03/30からのアクセス回数 4454

ここで紹介したSageワークシートは、以下のURLからダウンロードできます。

http://www15191ue.sakura.ne.jp:8000/home/pub/41/

また、Sageのサーバを公開しているサイト（http://www.sagenb.org/, http://www15191ue.sakura.ne.jp:8000/）にユーザIDを作成することで、ダウンロードしたワークシートをアップロードし、実行したり、変更していろいろ動きを試すことができます。

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

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

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

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

数式処理システムSageのノートブックは、計算結果を表示するだけではなく、実際に動かすことができるの大きな特徴です。この機会に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')

↑

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

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

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

sageへの入力:

# 3章のデータをもう一度読み込む
d = pd.read_csv('http://hosho.ees.hokudai.ac.jp/~kubo/stat/iwanamibook/fig/poisson/data3a.csv')

sageへの入力:

# statsmodelsを使ってglmを計算します
import statsmodels.formula.api as smf
import statsmodels.api as sm
from scipy.stats.stats import pearsonr

sageへの入力:

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

sageへの入力:

# λの予測値と一緒にプロット
ggplot(d, aes(x='x', y='y')) + geom_point() + geom_line(aes(x='x', y='lam_1'))

<ggplot: (39269997)>

sageへの入力:

ggsave('fig-4.1-a.png', dpi=50)

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

sageへの入力:

# 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_smooth(method=glm, family=poisson, formula =y ~ poly(x,7)) ")
r('plot(p)')
postGraph(graph)

↑

さまざまな逸脱度 †

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

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

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

sageへの入力:

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

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

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

sageへの入力:

# 逸脱度（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()

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

sageへの入力:

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

-192.889752524496

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

sageへの入力:

# 結果を小数点3桁で表示する
def prf3(f):
    return '%.3f'%f

sageへの入力:

# フルモデルは、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)

↑

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

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

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

↑

観測データが一つの場合 †

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

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

sageへの入力:

# サンプル生成λ=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)))

sageへの入力:

# 真のλ=8から50個のポアソン分布のサンプルを生成する
N = 50
lambda_true = 8
sample = genSample(lambda_true, N)
sampleDf = pd.DataFrame({'y' : np.array(sample)})

sageへの入力:

# 対数尤度の計算
def logL(sample, lambda_estimated):
    return sum([r.dpois(y, lambda_estimated, log=True) for y in sample]).n()

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

sageへの入力:

# k=1でのglm回帰を実行
fit = smf.glm('y ~ 1', data=sampleDf, family=sm.families.Poisson(link=sm.families.links.log)).fit()

sageへの入力:

# fit.summary2()
# 最大対数尤度
mxLogL = fit.llf
# 切片
beta_1 = fit.params[0]
print mxLogL, beta_1

-119.642879665 2.00148000021

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

sageへの入力:

# β= [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セット生成し、平均の対数尤度を計算します。

sageへの入力:

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

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

sageへの入力:

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()

sageへの入力:

# 実際の解析では不可能だが、今回は真のモデルが分かっているので、評価用のデータを生成し、E(logL)を計算できる

↑

繰り返しは簡単 †

次に上記の処理を繰り返した図４．９(B)を試してみます。

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

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

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

sageへの入力:

#  上記の計算を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 = [ r('sum(log(dpois(rpois(%d, %f), %f)))'%(N, lambda_true, exp(beta_1))) .n() 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の計算 †

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

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

sageへの入力:

# バイアス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 = [ r('sum(log(dpois(rpois(%d, %f), %f)))'%(N, lambda_true, exp(beta_1))) .n() for i in range(200)]
    return    (mxLogL - mean(logL200).n() )

sageへの入力:

# 非常に時間がかかります
bias_sample =  [bias(lambda_true, N) for i in range(200)]

sageへの入力:

mean(bias_sample)

(0.69534825087697705+0j)

sageへの入力:

# biasのヒストグラムを表示する
dd = pd.DataFrame({'y': np.array([ float(x) for x in bias_sample])})
ggplot(dd, aes(x='y'))  + geom_histogram()

binwidth defaulted to range/30. Use 'binwidth = x' to adjust this.
<ggplot: (39926417)>

sageへの入力:

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=１という値は、かなりぶれることも分かりました。

sageへの入力:

# 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

sageへの入力:

# サンプリングが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] 1.213076

sageへの入力:

graph = preGraph("fig-4.10-R.png")
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)

↑

コメント †

皆様のご意見、ご希望をお待ちしております。

sage/データ解析のための統計モデリング入門勉強ノート第4章

最新の5件

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

前準備 †

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

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

さまざまな逸脱度 †

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

観測データが一つの場合 †

繰り返しは簡単 †

バイアスbの計算 †

Rで計算 †

コメント †