※現在、ブログ記事を移行中のため一部表示が崩れる場合がございます。
順次修正対応にあたっておりますので何卒ご了承いただけますよう、お願い致します。
2017年 03月 10日
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つがあります。 `x` や `c` と関係のある無しで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"))後はサンプリングのための関数をこの関数に渡して、結果を表示しています。 得られた結果の解説は次回に回します。