TIM Labs

sato.yuichiroによるエントリー一覧

このブログ記事は『データ解析のための統計モデリング入門』(久保拓弥 著、岩波書店)という、とても分かりやすい統計モデリングの入門書を、さらに分かりやすくするために読書メモをまとめたものです。 今回は第8章、マルコフ連鎖モンテカルロ(MCMC)法とベイズ統計モデルについてのまとめの四回目です。

この章ではMCMCやメトロポリス法について説明がされています。 その中で、メトロポリス法によるサンプリングと尤度が関係あることが解説されています。 なので、メトロポリス法によるサンプリングの分布と尤度を比べられるコードを用意しました。 コードはRで書きました。

dataSize <- 20
N <- 8

data <- rbinom(dataSize, N, 0.45)

likelihood <- function (q) {
  sum(log(dbinom(data, N, q)))
}

develop <- function(q) {
  if (runif(1, 0, 1) < 0.5) {
    q + 0.01
  } else {
    q - 0.01
  }
}

sampling <- function(n, init) {
  q <- init
  ql <- likelihood(q)
  samples <- vector(length = n)
  i <- 1
  while(i <= n) {
    p <- develop(q)   
    pl <- likelihood(p)
    r <- log(runif(1, 0, 1))
    if (ql < pl || r < pl - ql) {
      q <- p
      ql <- pl
    }
    samples[i] <- p
    i <- i + 1
  }
  samples
}

qs <- seq(0, 1, by = 0.01)
ls <- vector(length = length(qs))
for(i in 1:length(qs)) {
  ls[i] <- prod(dbinom(data, N, qs[i]))
}
plot(qs, ls, type = "l", xlab = "q", ylab = "Likelihood")

samples <- sampling(10000, 0.3)
hist(samples, breaks = seq(0, 1, by = 0.01), main = "10000 samples", xlab = "q")

サンプルデータとしては、本の中でやっているのと同じように、事象の発生確率が 0.45 で、事象の試行回数は 8 の二項分布から、20 個のサンプルを取って用意しました。 コードを実行すると次のような図がプロットされます。

metro_likelihood.png

この図は q を変化させた時の尤度のグラフです。 発生させたサンプルデータの偏り方によって多少ずれますが、だいたい 0.45 を中心とした山なりなグラフになります。

metro_hist.png

この図はサンプル数 10000 の場合のメトロポリス法によるサンプリングの結果のヒストグラムです。 確かに、さきほどプロットした尤度に比例した分布が得られていることが分かります。

このブログ記事は『データ解析のための統計モデリング入門』(久保拓弥 著、岩波書店)という、とても分かりやすい統計モデリングの入門書を、さらに分かりやすくするために読書メモをまとめたものです。 今回は第8章、マルコフ連鎖モンテカルロ(MCMC)法とベイズ統計モデルについてのまとめの三回目です。

この章ではMCMCやメトロポリス法について説明がされています。 その中で、メトロポリス法によるサンプリングが解説されています。 なので、同じようにメトロポリス法でサンプリングしてみるコードを用意しました。 コードはRで書きました。

dataSize <- 20
N <- 8

data <- rbinom(dataSize, N, 0.45)

likelihood <- function (q) {
  sum(log(dbinom(data, N, q)))
}

develop <- function(q) {
  if (runif(1, 0, 1) < 0.5) {
    q + 0.01
  } else {
    q - 0.01
  }
}

sampling <- function(n, init) {
  q <- init
  ql <- likelihood(q)
  samples <- vector(length = n)
  i <- 1
  while(i <= n) {
    p <- develop(q)   
    pl <- likelihood(p)
    r <- log(runif(1, 0, 1))
    if (ql < pl || r < pl - ql) {
      q <- p
      ql <- pl
    }
    samples[i] <- p
    i <- i + 1
  }
  samples
}

samples1 <- sampling(500, 0.3)
hist(samples1, breaks = seq(0, 1, by = 0.01), main = "500 samples", xlab = "q")

samples2 <- c(samples1, sampling(500, samples1[500]))
hist(samples2, breaks = seq(0, 1, by = 0.01), main = "1000 samples", xlab = "q")

samples3 <- c(samples2, sampling(1000, samples2[1000]))
hist(samples3, breaks = seq(0, 1, by = 0.01), main = "2000 samples", xlab = "q")

samples4 <- c(samples3, sampling(8000, samples3[2000]))
hist(samples4, breaks = seq(0, 1, by = 0.01), main = "10000 samples", xlab = "q")

サンプルデータとしては、本の中でやっているのと同じように、事象の発生確率が 0.45 で、事象の試行回数は 8 の二項分布から、20 個のサンプルを取って用意しました。 なので、このデータに対して q のサンプリングをすると、だいたい 0.45 を中心とした分布が得られるはずです。

コードではサンプル数が 5001000200010000 の時にヒストグラムをプロットしています。 コードを実行すると次のような図がプロットされます。

metro_500.png

この図はサンプル数 500 の場合です。 まだヒストグラムが凸凹しているのが分かります。

metro_10000.png

この図はサンプル数 10000 の場合です。 だいぶヒストグラムがなめらかになって、分布が収束してきているのが分かります。

このブログ記事は『データ解析のための統計モデリング入門』(久保拓弥 著、岩波書店)という、とても分かりやすい統計モデリングの入門書を、さらに分かりやすくするために読書メモをまとめたものです。 今回は第8章、マルコフ連鎖モンテカルロ(MCMC)法とベイズ統計モデルについてのまとめの二回目です。

この章ではMCMCやメトロポリス法について説明がされています。 なので、実際にメトロポリス法を動かしてみるコードを用意しました。 コードはRで書きました。

dataSize <- 20
N <- 8

data <- rbinom(dataSize, N, 0.45)
summary(data)
cat("Analytically, q =", sum(data) / (dataSize * N), "\n\n")

likelihood <- function (q) {
  sum(log(dbinom(data, N, q)))
}

develop <- function(q) {
  if (runif(1, 0, 1) < 0.5) {
    q + 0.01
  } else {
    q - 0.01
  }
}

q <- round(0.3 + runif(1, 0, 0.3), 2)
ql <- likelihood(q)
i <- 0
cat("start q =", q, "\n")
while(i < 30) {
  p <- develop(q)   
  pl <- likelihood(p)
  r <- log(runif(1, 0, 1))
  if (ql < pl || r < pl - ql) {
    q <- p
    ql <- pl
  }
  cat("q =", q, "likelihood", ql, "\n")
  i <- i + 1
}
cat("last q =", q, "\n")

サンプルデータとしては、本の中でやっているのと同じように、事象の発生確率が 0.45 で、事象の試行回数は 8 の二項分布から、20 個のサンプルを取って用意しました。 このデータに対して、事象の発生確率 q が分かっていないものとして推定していきます。

メトロポリス法ではランダムに q を変化させた後、尤度が低くなる場合でも、変化の前後で尤度の比を計算して、その確率で q を更新します。 今の場合、計算の途中で使っているのは対数尤度なので、比ではなくて差になります。 コードで言うと以下の部分です。

  r <- log(runif(1, 0, 1))
  if (ql < pl || r < pl - ql) {
    q <- p
    ql <- pl
  }

コードを実行すると次のような出力が得られます。

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
    2.0     3.0     4.0     3.7     4.0     6.0
Analytically, q = 0.4625

start q = 0.31
q = 0.31 likelihood -40.44572
q = 0.32 likelihood -39.35181
q = 0.32 likelihood -39.35181
q = 0.32 likelihood -39.35181
q = 0.33 likelihood -38.34881
q = 0.33 likelihood -38.34881
q = 0.34 likelihood -37.43294
q = 0.35 likelihood -36.60087
q = 0.36 likelihood -35.84958
q = 0.37 likelihood -35.17642
q = 0.38 likelihood -34.579
q = 0.39 likelihood -34.05522
q = 0.4 likelihood -33.60322
q = 0.39 likelihood -34.05522
q = 0.4 likelihood -33.60322
q = 0.41 likelihood -33.22138
q = 0.42 likelihood -32.90828
q = 0.42 likelihood -32.90828
q = 0.43 likelihood -32.66271
q = 0.43 likelihood -32.66271
q = 0.44 likelihood -32.48365
q = 0.43 likelihood -32.66271
q = 0.44 likelihood -32.48365
q = 0.43 likelihood -32.66271
q = 0.44 likelihood -32.48365
q = 0.45 likelihood -32.37025
q = 0.46 likelihood -32.32184
q = 0.45 likelihood -32.37025
q = 0.44 likelihood -32.48365
q = 0.43 likelihood -32.66271
last q = 0.43

この例の場合は、最尤推定値を解析的に求めることができるので、まず最初にその値を出力しています。 メトロポリス法は確率的に尤度が低くなる変化もするので、q は一直線に最尤推定値に向かうのではなく、時々戻ったりするのが分かります。

このブログ記事は『データ解析のための統計モデリング入門』(久保拓弥 著、岩波書店)という、とても分かりやすい統計モデリングの入門書を、さらに分かりやすくするために読書メモをまとめたものです。 今回は第8章、マルコフ連鎖モンテカルロ(MCMC)法とベイズ統計モデルについてのまとめの一回目です。

この章ではMCMCやメトロポリス法について説明がされています。 本の中ではMCMCの前に、まずはランダムにパラメータを更新していく手法が紹介されています。 なので、本の中で紹介されているような推定法を実際に動かしてみるコードを用意しました。 コードはRで書きました。

dataSize <- 20
N <- 8

data <- rbinom(dataSize, N, 0.45)
summary(data)
cat("Analytically, q =", sum(data) / (dataSize * N), "\n\n")

likelihood <- function (q) {
  sum(log(dbinom(data, N, q)))
}

develop <- function(q) {
  if (runif(1, 0, 1) < 0.5) {
    q + 0.01
  } else {
    q - 0.01
  }
}

q <- round(0.3 + runif(1, 0, 0.3), 2)
ql <- likelihood(q)
i <- 0
cat("start q =", q, "\n")
while(i < 30) {
  p <- develop(q)   
  pl <- likelihood(p)
  if (ql < pl) {
    q <- p
    ql <- pl
  }
  cat("q =", q, "likelihood", ql, "\n")
  i <- i + 1
}
cat("last q =", q, "\n")

サンプルデータとしては、本の中でやっているのと同じように、事象の発生確率が 0.45 で、事象の試行回数は 8 の二項分布から、20 個のサンプルを取って用意しました。 このデータに対して、事象の発生確率 q が分かっていないものとして推定していきます。

推定方法は q をランダムに変化させて、対数尤度を計算し、尤度が大きくなっていたら変化させた先の q に更新する、という方法です。 develop 関数の中で q0.01 足すか引くかして q を変化させています。 コードを実行すると次のような出力が得られます。

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
      1       3       4       4       5       6
Analyticaly, q = 0.5

start q = 0.35
q = 0.36 likelihood -41.69397
q = 0.37 likelihood -40.76192
q = 0.38 likelihood -39.90848
q = 0.39 likelihood -39.13129
q = 0.39 likelihood -39.13129
q = 0.39 likelihood -39.13129
q = 0.39 likelihood -39.13129
q = 0.4 likelihood -38.42821
q = 0.4 likelihood -38.42821
q = 0.41 likelihood -37.79737
q = 0.42 likelihood -37.23712
q = 0.42 likelihood -37.23712
q = 0.43 likelihood -36.74602
q = 0.44 likelihood -36.32282
q = 0.45 likelihood -35.96647
q = 0.45 likelihood -35.96647
q = 0.45 likelihood -35.96647
q = 0.46 likelihood -35.67609
q = 0.47 likelihood -35.45097
q = 0.48 likelihood -35.29055
q = 0.49 likelihood -35.19445
q = 0.5 likelihood -35.16245
q = 0.5 likelihood -35.16245
q = 0.5 likelihood -35.16245
q = 0.5 likelihood -35.16245
q = 0.5 likelihood -35.16245
q = 0.5 likelihood -35.16245
q = 0.5 likelihood -35.16245
q = 0.5 likelihood -35.16245
q = 0.5 likelihood -35.16245
last q = 0.5

この例の場合は、最尤推定値を解析的に求めることができるので、まず最初にその値を出力しています。 q をランダムに変化させていくと、最終的にはその値に近づいていくのが分かります。

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

この章では複数の分布を混ぜて使う、一般化線形混合モデルついて説明がされています。 なので、サンプルデータを用意して、実際に一般化線形混合モデルを最尤推定してみるコードを用意しました。 コードはRで書きました。

library('glmmML')

dataSize <- 100
N <- 20

linkFunction <- function(b1, b2, x) {
  function(r) {
    1 / (1 + exp(-b1 - b2 * x - r))
  }
}

createData <- function() {
  x <- runif(dataSize, 0, 5)
  y <- vector(length = dataSize)
  for (i in 1:dataSize) {
    r <- rnorm(1, mean = 0, sd = 2)
    l <- linkFunction(-2, 1, x[i])
    y[i] <- rbinom(1, N, l(r))
  }
  id <- 1:dataSize
  data.frame(x, y, id)
}

data <- createData()
glmmML(cbind(y, N - y) ~ x, data = data, family = binomial, cluster = id)

コードを実行するには、glmmMLパッケージをインストールする必要があります。

リンク関数は logit(q) = b1 + b2 * x + r です。 r は平均 0 、標準偏差 s の正規分布から生成されます。

サンプルデータは、0 から 5 までの一様分布から x を生成し、標準偏差 s = 2 の正規分布から r を生成し、b1 = -2b2 = 1 のリンク関数を使って事象の発生確率 q を求め、試行回数 N = 20 の二項分布から値を生成しています。 なので、最尤推定すると、b1 = -2b2 = 1s = 2 に近い値を推定してくれるはずです。

コードを実行すると以下のような出力が出力されます。

Call:  glmmML(formula = cbind(y, N - y) ~ x, family = binomial, data = data,      cluster = id)


              coef se(coef)      z Pr(>|z|)
(Intercept) -2.163   0.4570 -4.732 2.22e-06
x            1.210   0.1604  7.539 4.73e-14

Scale parameter in mixing distribution:  2.144 gaussian
Std. Error:                              0.1654

        LR p-value for H_0: sigma = 0:  5.114e-131

Residual deviance: 311.6 on 97 degrees of freedom       AIC: 317.6

期待通りにパラメータが推定されているのが分かります。

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

この章では複数の分布を混ぜて使う、一般化線形混合モデルついて説明がされています。 その中で二項分布とと正規分布を混ぜて使う一般化線形混合モデルについて説明されています。 なので、サンプルデータを用意して、様々なパラメータにおける対数尤度の値をプロットしてくれるコードを用意しました。 コードはRで書きました。

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 = -2b2 = 1 のリンク関数を使って事象の発生確率 q を求め、試行回数 N = 20 の二項分布から値を生成しています。 なので、様々なパラメータについてモデルを試すと、b1 = -2b2 = 1s = 1 の時に最もデータへの当てはまりがよくなるはずです。 つまり尤度が最大になるはずです。

コードを実行すると5枚のグラフがプロットされます。 例えば次のようなグラフが得られます。

mixture_likelihood.png

図を見ると、期待通りに b1 = -2b2 = 1s = 1 の場合に尤度が最大になっているのが分かります。 実際の問題では、これらのパラメータを最尤推定すればいいわけです。

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

この章では複数の分布を混ぜて使う、一般化線形混合モデルついて説明がされています。 その中でポアソン分布と正規分布を混ぜた分布について説明されています。 なので、ポアソン分布と正規分布を混ぜ合わせた混合分布をグラフにしてくれるコードを用意しました。 コードはRで書きました。

linkFunction <- function(b1, b2, x) {
  function(r) {
    exp(b1 + b2 * x + r)
  }
}

distribution <- function(s, linkFunction) {
  l <- match.fun(linkFunction)
  function(y) {
    sum <- 0
    for (r in seq(-100, 100, by = 0.1)) {
      sum <- sum + dpois(y, l(r)) * dnorm(r, mean = 0, sd = s) * 0.1
    }
    sum
  }
}

xl <- "y"
yl <- "Probability"
xs <- 1:5
ys <- 0:10
b1 <- 0.1
legends = paste0("x = ", xs)
for (b2 in seq(0.1, 0.5, by = 0.1)) {
  for (s in c(0.5, 1.0, 1.5, 2.0)) {
    title <- paste0("lambda = exp(", b1, " + ", b2, " * x + r, r ~ N(0, ", s, ")")
    plot(0, 0, type = "n", xlim = c(0, 10), ylim = c(0.0, 0.5), main = title, xlab = xl, ylab = yl)
    for (i in 1:5) {
      l <- linkFunction(b1, b2, xs[i])
      d <- distribution(s, l)
      lines(ys, d(ys), type = "l")
      points(ys, d(ys), pch = i)
    }
    legend("topright", legend = legends, pch = 1:5)
  }
}

リンク関数は λ = exp(b1 + b2 * x + r) です。 r は平均 0 、標準偏差 s の正規分布から生成されます。 r を生成する正規分布の標準偏差が0に近い数の時は、これはポアソン回帰のリンク関数とほとんど同じです。 混合分布はポアソン分布と正規分布をかけて r で積分したものです。

コードを実行すると20枚のグラフがプロットされます。 例えば次のようなグラフが得られます。

mixture_poisson.png

この図は b1 = 0.1b2 = 0.5s = 1 の場合です。 ポアソン分布よりも広く分布した分布が得られているのが分かります。

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

この章では複数の分布を混ぜて使う、一般化線形混合モデルついて説明がされています。 その中でポアソン分布と正規分布を混ぜた分布について説明されています。 なので、ポアソン分布と正規分布を掛けた値をグラフにしてくれるコードを用意しました。 コードはRで書きました。

linkFunction <- function(b1, b2, x) {
  function(r) {
    exp(b1 + b2 * x + r)
  }
}

distribution <- function(y, s, linkFunction) {
  l <- match.fun(linkFunction)
  function(r) {
    dpois(y, l(r)) * dnorm(r, mean = 0, sd = s)
  }
}

xl <- "r"
yl <- "Probability"
ys <- seq(0, 8, by = 2)
rs <- seq(-5, 5, by = 0.5)
b1 <- 0.1
legends = paste0("y = ", ys)
for (b2 in seq(0.1, 0.5, by = 0.1)) {
  for (s in c(0.5, 1.0, 1.5, 2.0)) {
    l <- linkFunction(b1, b2, 2)
    title <- paste0("lambda  = exp(", b1, " + ", b2, " * 2 +  r), r ~ N(0, ", s, ")")
    plot(0, 0, type = "n", xlim = c(-5, 5), ylim = c(0.0, 0.2), main = title, xlab = xl, ylab = yl)
    for (i in 1:5) {
      d <- distribution(ys[i], s, l)
      lines(rs, d(rs), type = "l")
      points(rs, d(rs), pch = i)
    }
    legend("topright", legend = legends, pch = 1:5)
  }
}

リンク関数は λ = exp(b1 + b2 * x + r) です。 r は平均 0、標準偏差 s の正規分布から生成される変数です。 あとは、正規分布から r が生成される確率と、平均 λ のポアソン分布を掛け合わせれば、欲しい値になります。

コードを実行すると20枚のグラフがプロットされます。 例えば次のようなグラフが得られます。

mixed_poission.png

y によってグラフの山の大きさが変わってくるのが分かります。 混合分布はこれを r で積分したものです。 なので、グラフ上で面積の大きな y が確率の高い分布になります。

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

この章では複数の分布を混ぜて使う、一般化線形混合モデルついて説明がされています。 なので、二項分布と正規分布を混ぜ合わせた混合分布をグラフにしてくれるコードを用意しました。 コードはRで書きました。

N <- 20

linkFunction <- function(b1, b2, x) {
  function(r) {
    1 / (1 + exp(-b1 - b2 * x - r))
  }
}

distribution <- function(s, linkFunction) {
  l <- match.fun(linkFunction)
  function(y) {
    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
  }
}

xl <- "y"
yl <- "Probability"
xs <- c(-5,-2,0,2,5) 
ys <- 0:20
legends = paste0("x = ", xs)
for (b1 in seq(-0.5, 0.5, by = 0.2)) {
  for (b2 in seq(-0.5, 0.5, by = 0.2)) {
    for (s in c(0.1, 0.5, 1, 1.5, 2)) {
      title <- paste0("logit(q) = ", b1, " + ", b2, " * x + r, r ~ N(0, ", s, ")")
      plot(0, 0, type = "n", xlim = c(0, 20), ylim = c(0.0, 0.5), main = title, xlab = xl, ylab = yl)
      for (i in 1:5) {
        l <- linkFunction(b1, b2, xs[i])
        d <- distribution(s, l)
        lines(ys, d(ys), type = "l")
        points(ys, d(ys), pch = i)
      }
      legend("topright", legend = legends, pch = 1:5)
    }
  }
}

リンク関数は logit(q) = b1 + b2 * x + r です。 r は平均 0 、標準偏差 s の正規分布から生成されます。 r を生成する正規分布の標準偏差が0に近い数の時は、これはロジスティック回帰のリンク関数とほとんど同じです。

混合分布は二項分布と正規分布をかけて r で積分したものです。 今は単純に r-100 から 100 まで 0.1 刻みで変化させて、面積の合計を求めておきました。 コードで言うと以下の部分です。

    for (r in seq(-100, 100, by = 0.1)) {
      sum <- sum + dbinom(y, N, l(r)) * dnorm(r, mean = 0, sd = s) * 0.1
    }

これで様々な b1b2s について混合分布のグラフをプロットできます。 コードを実行すると180枚のグラフがプロットされます。 例えば次のようなグラフが得られます。

mixture_distribution.png

この図は b1 = 0.5b2 = 0.5s = 1 の場合です。 二項分布とはたいぶ違った分布が得られているのが分かります。

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

この章では複数の分布を混ぜて使う、一般化線形混合モデルついて説明がされています。 なので、二項分布と正規分布を掛けた値をグラフにしてくれるコードを用意しました。 コードはRで書きました。

N <- 20

linkFunction <- function(b1, b2, x) {
  function(r) {
    1 / (1 + exp(-b1 - b2 * x - r))
  }
}

distribution <- function(y, linkFunction) {
  l <- match.fun(linkFunction)
  function(r) {
    dbinom(y, N, l(r)) * dnorm(r, mean = 0, sd = 2)
  }
}

xl <- "r"
yl <- "Probability"
ys <- seq(0, 20, by = 5)
rs <- seq(-5, 5, by = 0.5)
legends = paste0("y = ", ys)
for (b1 in seq(-0.5, 0.5, by = 0.2)) {
  for (b2 in seq(-0.5, 0.5, by = 0.2)) {
    for (x in 0:3) {
      l <- linkFunction(b1, b2, x)
      title <- paste0("logit(q) = ", b1, " + ", b2, " * ", x, " + r, r ~ N(0, 2)")
      plot(0, 0, type = "n", xlim = c(-5, 5), ylim = c(0.0, 0.1), main = title, xlab = xl, ylab = yl)
      for (i in 1:5) {
        d <- distribution(ys[i], l)
        lines(rs, d(rs), type = "l")
        points(rs, d(rs), pch = ys[i])
      }
      legend("topright", legend = legends, pch = ys)
    }
  }
}

リンク関数は logit(q) = b1 + b2 * x + r です。 r は平均 0、標準偏差 2 の正規分布から生成される変数です。 事象の試行回数は 20 にしました。 リンク関数を使って二項分布する事象の発生確率 q が求まれば、事象が y 回起こる確率を計算できます。 あとは、正規分布から r が生成される確率と、二項分布を掛け合わせれば、欲しい値になります。

様々な b1b2x を与えて、r に対して確率をプロットしてみました。 コードを実行すると144枚のグラフがプロットされます。 例えば次のようなグラフが得られます。

mixed_2.png

y によってグラフの山の大きさが変わってくるのが分かります。 混合分布はこれを r で積分したものです

このアーカイブについて

このページには、sato.yuichiroが最近書いたブログ記事が含まれています。

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