LoginSignup
4

More than 5 years have passed since last update.

Rでガウス分布の平均・精度パラメータを逐次ベイズ推定

Last updated at Posted at 2013-04-13

PRML 2.3.6に記載の通り、ガウス分布に従う観測値xが与えられたとき、そのガウス分布の平均・精度パラメータの共役事前分布を用いて、xの観測毎に事後分布のパラメータを逐次求めます。

以下のベイズ推定をそれぞれ行います。

  • xの精度が既知として、共役事前分布としてガウス分布p(u)=N(u|nu,nv)を用いて、xの平均パラメータの事後分布を求めます。
  • xの平均が既知として、共役事前分布としてガンマ分布p(λ)=Gam(λ|ga,gb)を用いて、xの精度パラメータの事後分布を求めます。
  • xの平均と精度が未知として、共役事前分布としてガウス-ガンマ分布(正規-ガンマ分布)p(u,λ)=N(u|ngu,(ngbeta*λ)^(-1))Gam(λ|nga,ngb)を用いて、xの平均・精度パラメータの同時事後分布を求めます。

xのプロットに重ねて、以下を示します。

  • xの精度が既知として、xの平均パラメータの事後分布の期待値と標準偏差
  • xの平均が既知として、xの精度パラメータの事後分布の期待値と標準偏差を、xにおける標準偏差に変換して既知のxの平均パラメータに足したもの
  • xの平均と精度が未知として、xの平均・精度パラメータの事後分布の平均周辺分布の期待値と標準偏差
  • xの平均と精度が未知として、xの平均・精度パラメータの事後分布の精度周辺分布の期待値と標準偏差を、xにおける標準偏差に変換して、事後分布の平均周辺分布の期待値に足したもの

また、一定の観測毎に、その時点の平均パラメータ事前分布、精度パラメータ事前分布、平均・精度パラメータ同時事前分布を、その事前分布の期待値と標準偏差、およびxを生成した真の値とともに示します。

なお、事前分布には、無情報事前分布に近くなるよう、大きな分散を持つ分布を設定しています。

frame()
set.seed(0)
par(mfrow=c(4, 6))
layout(matrix(c(rep(1, 6), t(matrix(2:(3 * 6 + 1), 3))), 4, byrow=T), heights=c(3, 1, 1, 1))
par(mar=c(2.3, 2.5, 1, 0.1))
par(mgp=c(1.3, .5, 0))
urange <- c(148, 152)
lrange <- c(0, 0.5)
x1 <- seq(urange[1], urange[2], 0.1)
x2 <- seq(lrange[1], lrange[2], 0.01)

# サンプルの生成
N <- 20
U <- 150
V <- 4
x <- rnorm(N, U, sqrt(V))
xrange <- c(min(x), max(x))
# 平均パラメータの事前分布
nu <- 0
nv <- 1.0E+6  # 平均パラメータの事前分布の分散nvを大きく取る
# 精度パラメータの事前分布
ga <- 1
gb <- 1.0E-6  # 精度パラメータの事前分布の分散ga/gbを大きく取る
# 平均・精度パラメータの事前分布
ngu <- 0
ngbeta <- 1.0E-6  # 平均パラメータの事前分布の分散(ngbeta*λ)^-1を大きく取る
nga <- 1
ngb <- 1.0E-6  # 精度パラメータの事前分布の分散nga/ngbを大きく取る

d <- data.frame(nu=nu, nv=nv, ga=ga, gb=gb, ngu=ngu, ngbeta=ngbeta, nga=nga, ngb=ngb)
for (n in 1:N) {
    Nml <- 1
    uml <- x[n] / Nml
    vml <- (x[n] - U) ^ 2 / Nml
    # 精度が既知の時の平均パラメータuの事後分布N(u|nu,nv)を更新
    nu <- (V * nu + Nml * nv * uml) / (Nml * nv + V)
    nv <- 1 / (1 / nv + Nml / V)
    # 平均が既知の時の精度パラメータλの事後分布Gam(λ|ga,gb)を更新
    ga <- ga + Nml / 2
    gb <- gb + Nml * vml / 2
    # 平均・精度パラメータ(u,λ)の事後分布N(u|ngu,(ngbeta*λ)^(-1))Gam(λ|nga,ngb)を更新
    nga <- nga + Nml / 2
    ngb <- ngb + (Nml * 0 + ngbeta * Nml * (uml - ngu) ^ 2 / (ngbeta + Nml)) / 2
    ngu <- (Nml * uml + ngbeta * ngu) / (Nml + ngbeta)
    ngbeta <- ngbeta + Nml
    d <- rbind(d, c(nu, nv, ga, gb, ngu, ngbeta, nga, ngb))
}

# 精度パラメータ事前分布の期待値
d$gexp <- d$ga / d$gb
# 精度パラメータ事前分布の分散
d$gvar <- d$ga / d$gb ^ 2
# 平均・精度パラメータ事前分布の平均周辺分布の分散
d$nguvar <- d$ngb / (d$ngbeta * (d$nga - 1))
# 平均・精度パラメータ事前分布の精度周辺分布の期待値
d$nglexp <- d$nga / d$ngb
# 平均・精度パラメータ事前分布の精度周辺分布の分散
d$nglvar <- d$nga / d$ngb ^ 2

# 観測値
plot(x, pch=19, cex=1, xlim=c(1, N), ylim=xrange, xlab="n")
abline(h=U)
abline(h=U + sqrt(V) * c(-1, 1), lty=2)
# xの平均パラメータの事後分布の期待値と標準偏差
par(new=T)
plot(d$nu[-1], type="o", col=2, xlim=c(1, N), ylim=xrange, xlab="", ylab="", axes=F)
par(new=T)
plot((d$nu + sqrt(d$nv))[-1], type="l", lty=2, col=2, xlim=c(1, N), ylim=xrange, xlab="", ylab="", axes=F)
par(new=T)
plot((d$nu - sqrt(d$nv))[-1], type="l", lty=2, col=2, xlim=c(1, N), ylim=xrange, xlab="", ylab="", axes=F)
# xの精度パラメータの事後分布の期待値と標準偏差
par(new=T)
plot((U + sqrt(1 / d$gexp))[-1], type="o", col=3, xlim=c(1, N), ylim=xrange, xlab="", ylab="", axes=F)
par(new=T)
plot((U + sqrt(1 / (d$gexp + sqrt(d$gvar))))[-1], type="l", lty=2, col=3, xlim=c(1, N), ylim=xrange, xlab="", ylab="", axes=F)
par(new=T)
plot((U + sqrt(1 / (d$gexp - sqrt(d$gvar))))[-1], type="l", lty=2, col=3, xlim=c(1, N), ylim=xrange, xlab="", ylab="", axes=F)
# xの平均・精度パラメータの事後分布の平均周辺分布の期待値と標準偏差
par(new=T)
plot(d$ngu[-1], type="o", col=4, xlim=c(1, N), ylim=xrange, xlab="", ylab="", axes=F)
par(new=T)
plot((d$ngu - sqrt(d$nguvar))[-1], type="l", lty=2, col=4, xlim=c(1, N), ylim=xrange, xlab="", ylab="", axes=F)
par(new=T)
plot((d$ngu + sqrt(d$nguvar))[-1], type="l", lty=2, col=4, xlim=c(1, N), ylim=xrange, xlab="", ylab="", axes=F)
# xの平均・精度パラメータの事後分布の精度周辺分布の期待値と標準偏差
par(new=T)
plot((d$ngu + sqrt(1 / d$nglexp))[-1], type="o", col=6, xlim=c(1, N), ylim=xrange, xlab="", ylab="", axes=F)
par(new=T)
plot((d$ngu + sqrt(1 / (d$nglexp + sqrt(d$nglvar))))[-1], type="l", lty=2, col=6, xlim=c(1, N), ylim=xrange, xlab="", ylab="", axes=F)
par(new=T)
plot((d$ngu + sqrt(1 / (d$nglexp - sqrt(d$nglvar))))[-1], type="l", lty=2, col=6, xlim=c(1, N), ylim=xrange, xlab="", ylab="", axes=F)
legend("bottomright", legend=c(
    "p(u):E(u|X)±√(Var(u|X))", 
    "p(λ):150+√(1/(E(λ|X)±√(Var(λ|X))))", 
    "p(u,λ):E(u|X)±√(Var(u|X))", 
    "p(u,λ):E(u|X)+√(1/(E(λ|X)±√(Var(λ|X))))"), 
    col=c(2, 3, 4, 6), pch=1)

for (n in c(1, (1:5) * N / 5)) {
    # 平均パラメータ事前分布とその期待値と標準偏差
    curve(dnorm(x, d$nu[n], sqrt(d$nv[n])), col=2, xlim=urange, xlab="u", ylab="p(u)")
    abline(v=d$nu[n], col=2)
    abline(v=d$nu[n] + sqrt(d$nv[n]) * c(-1, 1), col=2, lty=2)
    abline(v=U)
    title(paste0("p(u)#", n), cex.main=1.0)
    # 精度パラメータ事前分布とその期待値と標準偏差
    curve(dgamma(x, d$ga[n], rate=d$gb[n]), col=3, xlim=lrange, xlab="λ", ylab="p(λ)")
    abline(v=d$gexp[n], col=3)
    abline(v=d$gexp[n] + sqrt(d$gvar[n]) * c(-1, 1), col=3, lty=2)
    abline(v=1/V)
    title(paste0("p(λ)#", n), cex.main=1.0)
    # 平均・精度パラメータ事前分布とその各周辺分布の期待値と標準偏差
    p <- outer(x1, x2, Vectorize(function(x1, x2) 
        dnorm(x1, d$ngu[n], sqrt(1 / (d$ngbeta[n] * x2))) * dgamma(x2, d$nga[n], rate=d$ngb[n])
        ))
    image(x1, x2, p, xlab=expression(mu), ylab=expression(lambda), col=rainbow(450)[256:1])
    abline(v=d$ngu[n], col=4)
    abline(v=d$ngu[n] + sqrt(d$nguvar[n]) * c(-1, 1), col=4, lty=2)
    abline(h=d$nglexp[n], col=6)
    abline(h=d$nglexp[n] + sqrt(d$nglvar[n]) * c(-1, 1), col=6, lty=2)
    points(U, 1/V)
    title(paste0("p(u,λ)#", n), cex.main=1.0)
}

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
4