「データ解析のための統計モデリング入門」読書ノート 7章無限混合分布の生成過程をアニメーションにした

一言で言うと

「データ解析のための統計モデリング入門」P.157 図7.8とP.158 図7.9の分布を混ぜる過程をアニメーションにした．無限混合分布に関して手元の資料であまり情報が無かったことがあり，検算のためにと計算したときに作ったものだが，意外とうまく作れたんじゃないかと思う．

GLMMの最尤推定と積分

個体差 $r_i$ を考慮に入れた今回のGLMMでは，尤度 $L_i$ を以下の式のように積分をして $r_i$ を消す．

$L_i = \int_{-\infty}^{\infty}p(y_i|\beta_1,\beta_2,r_i)p(r_i|s)dr_i$

これはすなわち， $p(y_i|\beta_1,\beta_2,r_i)$ と $p(r_i|s)$ を混ぜあわせた無限混合分布となっている．

アニメーション作成過程

今回は，P.157 図7.8とP.158 図7.9の2つの図をアニメーションとして表現してみた．それぞれの実験のパラメータは以下の通り．

二項分布と正規分布
- $\beta = 0$
- $q = \frac{1}{1+\exp(-(\beta + r))}$
ポアソン分布と正規分布
- $\beta = 0.5$
- $\lambda = \exp(\beta + r)$

ただし，無限混合分布に関しては，きちんと積分をして値を求めているわけではなく，rを細かく分割してそれぞれ計算した後に足しあわせているだけなので，あくまでも模式図ということで．

結果

gifアニメーションは少し容量が大きいので，以下のリンク先で公開している．

以下のプロットは，その中から1枚だけ抜き出した例となっている．左の図は二項分布/ポアソン分布の変化，中央の図は正規分布と $r$ の値，右の図は混合分布の変化を示している．rを小さい値から大きい値にかけて変化させていったときの，二項分布/ポアソン分布の変化と混合分布の形の変化がわかる．

以下ソースコード．

library(animation)

# 図 7.8
infinitemix <- function(){
  par(mfrow=c(1,3))
  q <- numeric(1)
  for(r in seq(-5, 5, length=100)){
    b <- dbinom(0:8, 8, (1/(1+exp(-r))))
    n <- dnorm(r, 0, 3)
    q <- q + b * n

    plot(0:8, b, type="b", ylim=c(0,1), xlab="y", main="Binom")
    curve(dnorm(x), from=-5, to=5, xlab="r", main="Normal")
    points(r, dnorm(r), cex=2)
    segments(r, 0, r, dnorm(r), lty=2)
    plot(0:8, q, type="b", ylim=c(0,2), xlab="y",main="Infinite mixture")
  }
}
saveMovie(infinitemix(), interval=0.1, moviename="",
          movietype="gif", outdir=getwd())

# 図 7.9
infinitemixpois <- function(){
  par(mfrow=c(1,3))
  q <- numeric(1)
  for(r in seq(-3, 3, length=100)){
    b <- dpois(0:10, exp(0.5+r))
    n <- dnorm(r, 0, 1)
    q <- q + b * n

    plot(0:10, b, type="b", ylim=c(0,1), xlab="y", main="Poisson")
    curve(dnorm(x), from=-4, to=4, xlab="r", main="Normal")
    points(r, dnorm(r), cex=2)
    segments(r, 0, r, dnorm(r), lty=2)
    plot(0:10, q, type="b", ylim=c(0,5), xlab="y", main="Infinite mixture")
  }
}
saveMovie(infinitemixpois(), interval=0.1, moviename="infinitemixpois",
          movietype="gif", outdir=getwd())