※現在、ブログ記事を移行中のため一部表示が崩れる場合がございます。
順次修正対応にあたっておりますので何卒ご了承いただけますよう、お願い致します。
2017年 05月 09日
linkFunction <- function(b1, b2) { function(x) { exp(b1 + b2 * log(x)) } } distribution <- function(x, s, linkFunction) { l <- match.fun(linkFunction) function(y) { dgamma(y, shape = s, rate = s / l(x)) } } xs <- c(0.1, 1.0, 2.0, 3.0, 4.0) ys <- seq(0.5, 10.0, by = 0.5) pchs <- 1:5 xl <- "y" yl <- "Probability" legends <- paste0("x = ", xs) for (b1 in seq(0.5, 2.5, by = 0.5)) { for (b2 in seq(0.5, 2.5, by = 0.5)) { for (s in 1:5) { l <- linkFunction(b1, b2) title <- paste0("shape = ", s, ", rate = ", s, " / exp(", b1, " + ", b2, "log(x))") plot(0, 0, type = "n", xlim = c(0, 10), ylim = c(0.0, 1.0), main = title, xlab = xl, ylab = yl) for (i in 1:5) { d <- distribution(xs[i], s, l) lines(ys, d(ys), type = "l") points(ys, d(ys), pch = pchs[i]) } legend("topright", legend = legends, pch = pchs) } } }リンク関数は `x` の単調増加関数です。 コードでは以下の部分です。
linkFunction <- function(b1, b2) { function(x) { exp(b1 + b2 * log(x)) } }これがガンマ分布の平均を与えます。 ガンマ分布にはshapeパラメータとrateパラメータがありますが、平均は `shape / rate` なのでリンク関数のパラメータだけを与えても、どんなガンマ分布になるか分かりません。 なので、リンク関数のパラメータ `b1` と `b2` の他にshapeパラメータも与えました。 コードを実行すると様々な `b1` と `b2` と `shape` の値に対してモデルをプロットします。 例えば下図のような図がプロットされます。