TIM Labs

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

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

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

createData <- function() {
  y <- rpois(10, lambda = 4)
  data.frame(y)
}

logLikelihood <- function(data, b1) {
  sum(log(dpois(data, lambda = exp(b1))))
}

printLogLikelihood <- function(data, i) {
  cat(paste0("Input data is data", i, "\n"))
  print(data[i]$y)

  fit <- glm(formula = y ~ 1, family = poisson, data = data[i])
  b1 <- fit$coefficients[1]

  cat("Estimated b1 =", b1, "\n")
  cat("Log likelihood for\n")
  cat("data1\t\tdata2\t\tdata3\t\n")
  cat(logLikelihood(data[1]$y, b1), "\t",
    logLikelihood(data[2]$y, b1), "\t",
    logLikelihood(data[3]$y, b1), "\n\n")
}

printData <- function(data, i) {
  cat(paste0("data", i))
  print(data[i]$y)
}

data <- c(createData(), createData(), createData())

for(i in 1:3) {
  printData(data, i)
}
cat("\n")
for(i in 1:3) {
  printLogLikelihood(data, i)
}

コードを実行すると、まずポアソン分布からサンプリングした 10 個のデータからなるデータセットを3個作ります。 それぞれのデータに対してポアソン回帰をして、推定されたパラメータを使って他の2つのデータセットに対して対数尤度を計算します。 対数尤度を計算しているコードは以下のコードです。

logLikelihood <- function(data, b1) {
  sum(log(dpois(data, lambda = exp(b1))))
}

b1 が最尤推定されたパラメータです。 data が対数尤度を計算する対象のデータです。 平均 exp(b1) を持つポアソン分布からデータが得られる確率の対数を計算しています。 実行結果は以下のようになります。

data1 [1] 5 1 2 4 9 3 3 2 3 4
data2 [1] 7 4 3 3 3 8 4 6 3 1
data3 [1] 9 1 5 6 1 3 3 8 3 2

Input data is data1
 [1] 5 1 2 4 9 3 3 2 3 4
Estimated b1 = 1.280934 
Log likelihood for
data1       data2       data3   
-20.59338    -21.43294   -24.32331 

Input data is data2
 [1] 7 4 3 3 3 8 4 6 3 1
Estimated b1 = 1.435085 
Log likelihood for
data1       data2       data3   
-21.04396    -20.95861   -24.00313 

Input data is data3
 [1] 9 1 5 6 1 3 3 8 3 2
Estimated b1 = 1.410987 
Log likelihood for
data1       data2       data3   
-20.91147    -20.97071   -23.99113 

傾向として、ポアソン回帰に使ったデータセットへの尤度が一番大きくて、他のデータセットへの尤度は小さくなります。 時々は、他のデータセットへの尤度が大きくなります。 最尤推定は与えられたデータに最もよく当てはまるパラメータを返すので、まだ見ぬデータにもよく当てはまるかどうかは分からないということです。

トラックバック(0)

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

コメントする

このブログ記事について

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

ひとつ前のブログ記事は「書評:"Learning Python 5th Edition"」です。

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

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