TIM Labs

データ解析のための統計モデリング入門 GLMのモデル選択 読書メモ5

このブログ記事は『データ解析のための統計モデリング入門』(久保拓弥 著、岩波書店)という、とても分かりやすい統計モデリングの入門書を、さらに分かりやすくするために読書メモをまとめたものです。 今回は第4章、GLMのモデル選択についてのまとめの五回目です。

この章ではバイアスやAICについて説明がされています。 データを説明するモデルとして様々なモデルが考えられる時に、データへのあてはまりの良さだけを基準にモデルを選ぶことはできません。 前回はAICを使ってモデル選択の様子を見るためのコードを紹介しました。 今回はそのコードの内容を解説します。

このコードを実行すると、4種類ある分布からサンプリングしたデータに対して、4種類のモデルを使ってポアソン回帰をし、AICやバイアスを計算して表示します。 計算されたAICを使って複数あるモデルから最適な物を選択する様子を見ることができます。

回帰に使うデータを作っているのは以下のコードです。 xc0 から 10 までの一様分布をします。 yxccreateY を通して何らかの関係を持つ数です。 xc が説明変数で y が応答変数です。

createData <- function(createY) {
  createY <- match.fun(createY)
  x <- runif(dataSize, 0, 10)
  c <- runif(dataSize, 0, 10)
  y <- createY(x, c)
  data.frame(x, c, y)
}

関係としては、以下の4つがあります。 xc と関係のある無しで4種類です。

function(x, c) rpois(dataSize, lambda = 4)
function(x, c) rpois(dataSize, lambda = exp(0.1 + 0.2 * x))
function(x, c) rpois(dataSize, lambda = exp(0.1 + 0.2 * c))
function(x, c) rpois(dataSize, lambda = exp(0.1 + 0.2 * x + 0.2 * c))

ポアソン回帰をやっているのは以下のコードです。 f でリンク関数を指定して、 genData によって生成されたデータ d に対してポアソン回帰をやっています。 推定されたパラメータを使いやすくするために、n という名前を付けています。

calculateAic <- function(d, genData, f, n) {
  fit <- glm(formula = f, family = poisson, data = d)
  co <- append(fit$coefficients, replicate(2, 0))
  names(co) <- n

このようにして得たパラメータを使って、最大対数尤度、平均対数尤度、バイアス、AICを以下のように計算しています。

  meanLogLik <- meanLogLikelihood(genData, co["b1"], co["b2"], co["b3"])
  bias <- logLik(fit) - meanLogLik
  aic <- -2 * (logLik(fit) - length(fit$coefficients))
  result <- c(logLik(fit), meanLogLik, bias, aic)
  names(result) <- c("logLik", "meanLogLik", "bias", "aic")
  result
}

これで1つのモデルに対して、1回分の実験が終わります。 以下のコードで実験を何度も繰り返して、結果を平均しています。

calculateMeanAic <- function(genData, f, n) {
  genData <- match.fun(genData)
  m <- replicate(expSize, calculateAic(genData(), genData, f, n))
  apply(m, 1, mean)
}

以下のコードで4種類のモデルに対して上記の計算をしています。 リンク関数は exp(b1)exp(b1 + b2 * x)exp(b1 + b3 * c)exp(b1 + b2 * x + b3 * c) の4種類です。

showAic <- function(genData) {
  r1 <- calculateMeanAic(genData, y ~ 1, c("b1","b2","b3"))
  r2 <- calculateMeanAic(genData, y ~ x, c("b1","b2","b3"))
  r3 <- calculateMeanAic(genData, y ~ c, c("b1","b3","b2"))
  r4 <- calculateMeanAic(genData, y ~ x + c, c("b1","b2","b3"))

後はサンプリングのための関数をこの関数に渡して、結果を表示しています。 得られた結果の解説は次回に回します。

トラックバック(0)

トラックバックURL: http://labs.timedia.co.jp/mt/mt-tb.cgi/582

コメントする

このブログ記事について

このページは、sato.yuichiroが2017年3月10日 12:00に書いたブログ記事です。

ひとつ前のブログ記事は「Chainer:超簡単なAutoEncoderを作ってみる(2)」です。

次のブログ記事は「Chainer:超簡単なAutoEncoderの脳の中を見てみよう」です。

最近のコンテンツはインデックスページで見られます。過去に書かれたものはアーカイブのページで見られます。