※現在、ブログ記事を移行中のため一部表示が崩れる場合がございます。
順次修正対応にあたっておりますので何卒ご了承いただけますよう、お願い致します。
2017年 06月 13日
dataSize <- 50 N <- 20 linkFunction <- function(b1, b2, x) { function(r) { 1 / (1 + exp(-b1 - b2 * x - r)) } } dmixture <- function(b1, b2, s) { function(y, x) { l <- linkFunction(b1, b2, x) sum <- 0 for (r in seq(-100, 100, by = 0.1)) { sum <- sum + dbinom(y, N, l(r)) * dnorm(r, mean = 0, sd = s) * 0.1 } sum } } likelihood <- function(b1, b2, s) { d <- dmixture(b1, b2, s) function (y, x) { d(y, x) } } logLikelihood <- function(data, b1, b2) { function (s) { sl <- 0 ll <- likelihood(b1, b2, s) for (i in 1:dataSize) { sl <- sl + log(ll(data$y[i], data$x[i])) } sl } } createData <- function() { x <- runif(dataSize, 0, 5) y <- vector(length = dataSize) for (i in 1:dataSize) { r <- rnorm(1, mean = 0, sd = 1) l <- linkFunction(-2, 1, x[i]) y[i] <- rbinom(1, N, l(r)) } data.frame(x, y) } data <- createData() xl <- "s" yl <- "Lilelihood" ss <- seq(0.2, 2, by = 0.2) b2s <- seq(-2, 2, by = 1) legends <- paste0("b2 = ", b2s) for (b1 in seq(-2, 2, by = 1)) { title <- paste0("logit(q) = ", b1, " + b2 * x + r, r ~ N(0, s)") plot(0, 0, type = "n", xlim = c(0, 2), ylim = c(-1000, 0), main = title, xlab = xl, ylab = yl) for (i in 1:5) { ls <- logLikelihood(data, b1, b2s[i])(ss) lines(ss, ls, type = "l") points(ss, ls, pch = i) } legend("topleft", legend = legends, pch = 1:5) }リンク関数は `logit(q) = b1 + b2 * x + r` です。 `r` は平均 `0` 、標準偏差 `s` の正規分布から生成されます。 サンプルデータを生成するには、まず `0` から `5` までの一様分布から `x` を、標準偏差 `s = 1` の正規分布から `r` を生成します。 その後、`b1 = -2`、`b2 = 1` のリンク関数を使って事象の発生確率 `q` を求め、試行回数 `N = 20` の二項分布から値を生成しています。 なので、様々なパラメータについてモデルを試すと、`b1 = -2`、`b2 = 1`、`s = 1` の時に最もデータへの当てはまりがよくなるはずです。 つまり尤度が最大になるはずです。 コードを実行すると5枚のグラフがプロットされます。 例えば次のようなグラフが得られます。