多変量正規分布と事後分布 - ryamadaの遺伝学・遺伝統計学メモ

こちらでGaussian Sequence Modelをなぞっている。その一環
メモ

n.iter <- 10000
# n個の観察をするとする
n <- 20
# n個の真の値の種となる値を適当に決める。正規分布にしてもよい
theta_0 <- rexp(n)
# 正定値行列を作る
library(GPArotation)
R <- Random.Start(n)
S <- runif(n)*5
Tau <- R %*% diag(S) %*% t(R)
library(mvtnorm)
# 正定値対称行列を分散共分散行列とする、期待値theta_0の多変量正規乱数をn.iter個(組)作る
theta <- rmvnorm(n.iter,theta_0,Tau)
#plot(theta,asp=1)
# この「真値」を観測すると、theta[i,]を中心として正定値対称行列Sigmaなる分散共分散行列を持つ誤差の入った値(の組)が得られる
# 正定値行列を作る
S2 <- runif(n)*7
R2 <- Random.Start(n)
Sigma <- R2 %*% diag(S2) %*% t(R2)
# その作り方として、theta[i,]ごとに乱数を発生させてつくる
y <- matrix(0,n.iter,n)
for(i in 1:n.iter){
	y[i,] <- rmvnorm(1,theta[i,],Sigma)
}
# theta_0を中心に二つのばらつき行列の和をばらつき情報として与えて乱数を発生させる
y2 <- rmvnorm(n.iter,theta_0,Tau+Sigma)
# ２つの発生のさせ方でできる分布は同じ
plot(y2,asp=1,pch=20)
points(y,col=2,pch=20)
# 平均と分散が同じ
y.m <- apply(y,2,mean)
y2.m <- apply(y2,2,mean)
plot(y.m,y2.m,asp=1)
abline(0,1,col=2)

y.v <- apply(y,2,var)
y2.v <- apply(y2,2,var)
plot(y.v,y2.v,asp=1)
abline(0,1,col=2)

ここで見ているのは、ある値ベクトルから定めたルール(行列)で正規分布ずれしている真値をさらにルール(行列)で正規分布ずれしているような観測値の分布
この観測を元にして、元の値ベクトル $\theta_0$ の事後分布を知りたいとする
それは $\theta_y = (\Sigma^{-1}+T^{-1})^{-1} (\Sigma^{-1}y + T^{-1} \theta_0)$ を中心といて、分散共分散行列が $\Sigma_y = (\Sigma^{-1}+T^{-1})^{-1}$ なる正規分布になる
この中心は、 $\Sigma^{-1},T^{-1}$ という、元の２つのばらつきを定める行列の逆元に関する重み付き平均であり
その分散共分散行列は調和平均っぽいものになっている
行列の演算なので見にくいけれど $Q^{-1}$ という逆元を「割り算すること・分母にすること」と読んでしまえば、 $\theta_y = \frac{\frac{1}{\Sigma} y + \frac{1}{T} \theta_0}{\frac{1}{\Sigma}+\frac{1}{T}}$ であって、 $\Sigma_y = \frac{1}{\frac{1}{\Sigma} + \frac{1}{T}}$
この式もよいけれど、場合によっては別の式表現が関連領域で用いられることもあるので、その式に従った計算も実施し、２方法での計算が一致することを確認している

Tau.inv <- solve(Tau)
Sigma.inv <- solve(Sigma)
Sigma.y <- solve(Sigma.inv + Tau.inv)
Sigma.y2 <- Tau - Tau %*% solve(Tau + Sigma) %*% Tau
range(Sigma.y-Sigma.y2)
theta.y <- Sigma.y %*% (Sigma.inv %*% y[1,] + Tau.inv %*% theta_0)
theta.y2 <- Tau %*% solve(Tau+Sigma) %*% y[1,] + Sigma %*% solve(Tau+Sigma) %*% theta_0
range(theta.y-theta.y2)
plot(theta_0,theta.y)
theta.post <- rmvnorm(n.iter,theta.y,Sigma.y)
apply(theta.post,2,mean)
theta_0
apply(theta.post,2,var)