TIM Labs

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

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

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

dataSize <- 100
expSize <- 1000

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

createIndipendentData <- function() {
  createData(function(x, c) rpois(dataSize, lambda = 4))
}

createCorrelatedXData <- function() {
  createData(function(x, c) rpois(dataSize, lambda = exp(0.1 + 0.2 * x)))
}

createCorrelatedCData <- function() {
  createData(function(x, c) rpois(dataSize, lambda = exp(0.1 + 0.2 * c)))
}

createCorrelatedXCData <- function() {
  createData(function(x, c) rpois(dataSize, lambda = exp(0.1 + 0.2 * x + 0.2 * c)))
}

logLikelihood <- function(genData, b1, b2, b3) {
  genData <- match.fun(genData)
  data <- genData()
  prob <- vector(length = dataSize)
  for (j in 1:dataSize) {
    prob[j] <- log(dpois(data$y[j], lambda = exp(b1 + b2 * data$x[j] + b3 * data$c[j])))
  }
  sum(prob)
}

meanLogLikelihood <- function(genData, b1, b2, b3) {
  mean(replicate(dataSize, logLikelihood(genData, b1, b2, b3)))
}

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

  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
}

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

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

  cat("\t\t\tnull\t\tx\t\tc\t\txc\n")
  cat("log likelihood\t\t", r1["logLik"], "\t", r2["logLik"], "\t",
    r3["logLik"], "\t", r4["logLik"], "\n")
  cat("mean log likelihood\t", r1["meanLogLik"], "\t",
     r2["meanLogLik"], "\t", r3["meanLogLik"], "\t", r4["meanLogLik"], "\n")
  cat("bias\t\t\t", r1["bias"], "\t", r2["bias"], "\t", r3["bias"], "\t", r4["bias"], "\n")
  cat("aic\t\t\t", r1["aic"], "\t", r2["aic"], "\t", r3["aic"], "\t", r4["aic"], "\n\n")
}

cat("indipendent data\n")
showAic(createIndipendentData)

cat("correlated X data\n")
showAic(createCorrelatedXData)

cat("correlated C data\n")
showAic(createCorrelatedCData)

cat("correlated XC data\n")
showAic(createCorrelatedXCData)

少し長くなったので、解説は次回に回します。 ぜひ実行してみてください。

トラックバック(0)

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

コメントする

このブログ記事について

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

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

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

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