※現在、ブログ記事を移行中のため一部表示が崩れる場合がございます。
順次修正対応にあたっておりますので何卒ご了承いただけますよう、お願い致します。
2017年 03月 03日
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傾向として、ポアソン回帰に使ったデータセットへの尤度が一番大きくて、他のデータセットへの尤度は小さくなります。 時々は、他のデータセットへの尤度が大きくなります。 最尤推定は与えられたデータに最もよく当てはまるパラメータを返すので、まだ見ぬデータにもよく当てはまるかどうかは分からないということです。