正規分布のパラメータ \(\mu, \Sigma\) の共役事前分布は Normal-inverse-wishart (NIW) 分布。データ数が増加した時に真のパラメータに収束していく様子を図示してみた (クリックで拡大する、かも)。
絵のせつめい。
- データ数 n を変化させた各グラフの中で、事後分布から20組のパラメータをサンプリング。
- 一組のパラメータを90%信頼楕円として表している
- 点線は真の分布。黒丸は事前分布のモード(最頻値)。
- データが少ないうちは事後分布がばらつく、つまりパラメータの不確実性が大きいが、データ数が256を超えたあたりからほぼ真の分布に収束する
下のコードのコメントにも書いたが、ハイパーパラメータを学習しない場合のレシピとして、NIW分布の平均パラメータ \(\mu_0\)、分散パラメータ \(V_0\) に自信がないときはそれぞれ \(k_0, \nu_0\) を小さく設定すればいい、はず。事前分布の影響を小さくできる。
参考その１：前回のエントリ → 逆Wishart分布を図示してみる
参考その２：事後分布のパラメータの求め方 → Wikipedia/Conjugate prior
お絵かきスクリプト in R。ellipse, MCMCpack, mvtnorm パッケージは入っていなければインストールする必要あり。
library(ellipse) # conffidence ellipse
library(MCMCpack) # wishart
library(mvtnorm)
#----------------------------------------------------------------------------
# Baysian estimation of 2 dimensional normal distribution
#----------------------------------------------------------------------------
d <- 2 # dimension
#----------------------------------------------------------------------------
# Hyperparameters for NIW (Normal inverse wishart)
k0 <- 0.1 # see below
mu0 <- rep(0,d) # hyper mean
v0 <- 3.5 # see below
V0 <- diag(rep(10,d)) # hyper variance
Phi0 <- V0*(v0-d-1)
prior.hyper.par <- list(k=k0,mu=mu0,v=v0,Phi=Phi0)
#------------------------------------------------------------------------------------
# [Recipe for determining hyperparameters v0 and k0]
# 1. Less confidence for the hyper mean(mu0), take smaller k0; typically, 0<k0<<1
# 2. Less confidence for the hyper variance(V0), take smaller v0; typically, v0~1+p
#------------------------------------------------------------------------------------
#-----------------------------------------
# Bayesian estimation ( bayesian update)
bayes.update <- function(X,hyper.par){
k0 <- hyper.par$k
mu0 <- hyper.par$mu
v0 <- hyper.par$v
Phi0 <- hyper.par$Phi
if(!is.matrix(X)){
# This is n=1 case that we have to deal with as special because of R issue
n <- 1
EX <- X
X <- t(X)
C <- matrix(rep(0,d*d),nc=d)
} else {
n <- dim(X)[1]
EX <- colMeans(X)
C <- (n-1)*cov(X)
}
mu <- (k0*mu0 + n*EX)/(k0+n)
k <- k0+n
v <- v0+n
Phi <- Phi0 + C + k0*n/(k0+n)*((EX-mu0) %*% t(EX-mu0))
list(k=k,mu=mu,v=v,Phi=Phi) # This new parameters list is a bayesian update result!!
}
#------------------------------------------------------------------
# Sampling from posterior
sample.posterior <- function(n,bayes.fit){
k <- bayes.fit$k
v <- bayes.fit$v
Phi <- bayes.fit$Phi
mu <- bayes.fit$mu
result <- list()
for(i in 1:n){
# Sampling mu and V from NIW
V <- riwish(v,Phi) # 1. sampling covariance matrix by inverse wishart
result[[i]] <- list()
result[[i]]$V <- V
result[[i]]$mu <- rmvnorm(1,mu,V/k) # 2. sampling mu by normal dist.
}
result
}
#--------------------------------------------------------------
# Test case
# --- true distribution
V.true <- matrix(c(3,-1.5,-1.5,1.8),nc=d)
mu.true <- rep(10,d)
# --- generating random variables
X <- rmvnorm(20000,mu.true,V.true)
par(mfrow=c(3,4)) # dividing graphic device
j <- 1
ids <- 2^seq(0,11)
hc <- rainbow(12)
for(i in ids){
Y <- X[1:i,]
bayes.fit <- bayes.update(Y,prior.hyper.par) # estimating posterior
post.samp <- sample.posterior(20,bayes.fit) # sampling from posterior
plot(Y,pch=".",col="gray",xlim=c(-10,20),ylim=c(-10,20),xlab="",ylab="") # just plotting raw data
# plotting posterior samples
for( theta in post.samp ){
mu <- theta$mu
V <- theta$V
elp <- t(apply(ellipse(V,level=0.9),1,function(x) x+mu))
lines(elp,col=hc[j],lw=2)
}
# plotting true distribution
true.dist <- t(apply(ellipse(V.true,level=0.9),1,function(x) x+mu.true))
lines(true.dist,type="l",lw=2,lt=4)
# plotting mode of the prior
prior <- ellipse(Phi0/(v0+3), level=0.9) + mu0
lines(prior,lw=2)
text(20,-7,paste("n =",i),pos=2,cex=2)
j <- j+1
}
dev2bitmap(file="normal_posterior.jpg",taa=4,gaa=4,width=15,height=12)
