わかりやすいギブスサンプリング

ここでは、マルコフ連鎖モンテカルロ法（MCMC法）について説明していきます。

マルコフ連鎖モンテカルロ法とは

マルコフ連鎖モンテカルロ法とは、略称でMCMC法と呼ばれたりします。
MCMC法では「解析が難しい事後分布のパラメータをサンプリングして求める」ための手法です。

解析や計算処理が難しい事後分布からパラメータ（以下例えば\(\mu, \sigma^{2}\)))を求めたい場合などで使用したりします。
本記事ではMCMC法の1つであるギブスサンプリング法について説明しようと思います。

ベイズの更新式で、共役事前分布・共役事後分布という性質のものがあります。
これらでは、上記の右辺の尤度・事前分布を計算して事後分布を求めると、事後分布は事前分布と同じ分布になるという性質のものです。
これらのメリットは、同じ分布になるので、正規分布やガンマ分布などよくある確率密度関数で計算ができ、事後分布のパラメータをよくある確率密度関数からRやpythonなどのライブラリですでに用意されているので、簡単に事後分布のパラメータのサンプリングを行うことができます。
 

できる限り、この共役分布の性質のものを見つけることが大事ですが、
事後分布が全く知らないような独自分布になってしまった場合は、どうにかして事後分布のパラメータをサンプリングをしなければいけません。
そこで必要なのが、MCMC法です。

 
なので、最終的なMCMC法の目標は、
「事後分布からパラメータをサンプリングする。」
です。
そして、マルコフ連鎖で不変分布となることです。
 

ギブスサンプリング

通常正規分布のように、確率変数xとパラメータがあります。
ここに正規分布\(N(\mu, \sigma^{2})\)があるとします。
 

ギブスサンプリングは以下のような用途で使われます。

\(①\)はベイズ推定での事後分布サンプリングに使われることがあります。
\(②\)は通常サンプリングは1次元になります。なので多次元の確率変数などになった場合、サンプリングが難しくなるのでこのギブスサンプリングを用いたりします。
 

ギブスサンプリンでは、メトロポリス・ヘイスティングス法とは違い、提案分布を用いることはしません。
サンプリングしたい分布関数のみを使用します。
ただ直接分布関数を使用しても、そもそもサンプリングができないから別の方法でサンプリングをしようと思っているので直接使用はできません。
そのため、式変形などをして、その式変形した結果からサンプリングが直接できるような関数であれば、このギブスサンプリングを用いることができます。
 

確率変数でもパラメータでも良いのですが、
サンプリングしたいパラメータや確率変数を1つ選択して、それ以外を条件とした上で、条件付き確率を定義します。
この条件付き確率から1つサンプリングをして、
次にまた別のサンプリングしたいパラメータや確率変数を1つ選択して、それ以外を条件とした上で、条件付き確率を定義します。先程サンプリングしたパラメータは今回は条件のほうに入ってるので、サンプリングした最新の値を使います。
 

このようにして繰り返しサンプリングを行っていく手法がギブスサンプリングになります。
そのため、条件付き確率から直接サンプリングを行う必要があり、ここができないとギブスサンプリングを用いてのサンプリングはできないということになります。
 

パラメータを求める必要があるわけで、
例えば、パラメータが\(\mu, \sigma^{2}\)だとして、確率変数は\(x\)とした時、
それぞれ交互にサンプリングしていく方法がギブスサンプリングになります。
当然この2つのパラメータが同時分布\(P(\mu, \sigma^{2}|x)\)からはサンプリングはできません。
なので、それぞれ\(\mu\)のみの関数、\(\sigma^{2}\)のみの関数にしていきます。
要は周辺分布ということですね！
なので、
\(P(\mu|x, \sigma^{2})\)
\(P(\sigma^{2}|x, \mu)\)
を求めます。
それぞれ、
\(P(\mu|x, \sigma^{2}) = \displaystyle \frac{P(x, \mu, \sigma^{2})}{P(x, \sigma^{2})} \)
\(P(\sigma^{2}|x, \mu) = \displaystyle \frac{P(x, \mu, \sigma^{2})}{P(x, \mu)}\)
を計算して求めます。
 

\begin{eqnarray}
\mu^{1} → (\sigma^{2})^{1} → \mu^{2} → (\sigma^{2})^{2} → ...
\end{eqnarray}
 

予測分布

もちろん、MCMC法の1つであるこのギブスサンプリングは事後分布はパラメータによる関数なので、このギブスサンプリングを用いることでパラメータを得ることができます。
確率変数\(x\)をサンプリングするのではなく、その\(x\)の確率分布であるパラメータをサンプリングします。
そして事前分布にて、その得たパラメータの出る確率を求め、そして尤度関数に適用します。
 

予測分布を定義します。
そもそも予測分布とは、今得られているデータから次に得られるデータを確率的に予測しようというものです。
 

今得られたデータ\(x_{n}\)から、次のデータ\(x_{n+1}\)を出したいとします。
この時、パラメータは\(\mu, \sigma^{2}\)とすると、
予測分布は以下のように設定できます。
\begin{eqnarray}
f(x_{n+1}|x_{n}) &=& \int\int f(x_{n+1}|\mu, \sigma^{2})f(\mu, \sigma^{2}|x_{n})d\mu d\sigma^{2}
\end{eqnarray}
 

この式の意味を考えていくと、
一番左はそれこそ次の\(x\)を求める予測分布で、現在のデータ\(x_{n}\)から次に取るデータ\(x_{n+1}\)を予測する条件付き分布になります。
一番右は今得られている\(x\)から求めた\(\mu\)と\(\sigma^{2}\)なので事後分布。
\(x\)→\(\mu\)と\(\sigma^{2}\)→次の\(x\)を求めるという形です。
 
事後分布である一番後ろは、今までのデータを条件にしてパラメータの推定を条件付き確率で構成してます。
今までのデータを用いて一番推定のよかったパラメータを使って、\(x_{n+1}\)の条件付き確率を作ってます。理にかなってます。
 

そして2つのパラメータを積分させることで、この2つのパラメータがとりうる全ての値で計算をして合算した値がこの予測分布になります。
実際この2つのパラメータで積分すれば右の項の真ん中は消えるので、
\(f(x_{n+1}|\mu, \sigma^{2})f(\mu, \sigma^{2}|x_{n})\)は、間の\(\mu, \sigma^{2}\)は消えて、\(f(x_{n+1}|x_{n})\)になります。
 

この計算をすることで、パラメータ\(\mu, \sigma^{2}\)は消えるので、あくまで\(x_{n}, x_{n+1}\)に依存した分布になります。
この定義式で良さそうですね！
 

予測分布の求め方\(①\)

実際これを計算して求めようとするとまずほぼほぼ無理です。
なのでどうするかというと、上記の積分の意味を考えると全ての\(\mu\)と\(\sigma^{2}\)の組み合わせでの和なので、
上のギブスサンプリングで得られた\(\mu\)と\(\sigma^{2}\)の組み合わせを以下のように定義します。
\begin{eqnarray}
(\mu^{(i)}, {\sigma^{2}}^{(i)})_{i=1,2,3,...}
\end{eqnarray}

iはサンプリングで得られた個数です。5個得たのであれば、\(i=1,2,3,4,5\)で、5組得られたことになります。
 
よって予測分布は以下のように表せます。
\begin{eqnarray}
f(x_{n+1}|x_{n}) &=& \sum_{i=1}^{5} f(x_{n+1}|\mu^{(i)}, {\sigma^{2}}^{(i)})f(\mu^{(i)}, {\sigma^{2}}^{(i)} |x_{n})
\end{eqnarray}
 

もちろん、5個ではダメで、burn_inや10000個などもっと大量のパラメータを生成して計算しないと精緻なものはできないです。
 

予測分布の求め方\(②\)

上の「予測分布の求め方\(①\)」では、パラメータの組\((\mu^{(i)}, {\sigma^{2}}^{(i)})_{i=1,2,3,...}\)をサンプリングしていきました。

なので、であれば、この\(\mu\)と\(\sigma^{2}\)の組について、
直接、\(x_{n+1}^{(i)} \sim N(\mu^{(i)}, {\sigma^{2}}^{(i)})\)なので、
直接それぞれの\(i\)の組み合わせで、サンプリングをして、
ヒストグラムを書いて、一番出る最頻値(mode)を採用すれば良いかと思います。
 

なので、予測分布の\(x_{n}\)を条件としたことは、間接的に\(x_{n}\)を用いて\(\mu\)と\(\sigma^{2}\)は求まっているので、
サンプリングでヒストグラム書いて求めればそれでいい！
 

またこの予測分布の関数が正規分布などよく知られているような確率分布でないのであれば、サンプリングは少々難しいです。
有名な分布であれば、pythonやRで簡単にサンプリングができます。
知られていない場合は、それこそ上で説明した、モンテカルロ法を用いてサンプリングを行います。
「モンテカルロ法はさまざまな確率分布からの乱数を発生させる手法」なので。
 

ギブスサンプリングでは提案分布を使用せず、サンプリングしたい関数（この関数からサンプリングを直接するのが難しい）の
式変形を行って、そこから直接サンプリングを行うことができれば、ギブスサンプリングを用いることが可能です。

そしてメトロポリス・ヘイスティングス法は、提案分布を使用して、サンプリングしたい関数（この関数からサンプリングを直接するのが難しい）からサンプリングをする。
提案分布はこちらが用意した普通にサンプリングができる関数。例えば正規分布とか。
そこからサンプリングしたデータについて、そのデータがサンプリングしたい関数での確率を求める。（代入して確率計算できる）
→ こうすることで、サンプリングしたい関数からのサンプリングはマルコフ連鎖となり、不変分布となる。

マルコフ連鎖では以下のような性質がある。

【具体例】実際にギブスサンプリングを行ってみる

具体例を考えてみます！

 
まずベイズの定理から、事後分布 = 尤度 × 事前分布であることから、
\(\pi(\mu, \sigma^{2} | x) \sim f(x|\mu, \sigma^{2}) \cdot \pi(\mu, \sigma^{2}) \)
であり、事前分布はベイズの定理から\(\pi(\mu, \sigma^{2}) = \pi(\mu | \sigma^{2})\pi(\sigma^{2}) \)となるので、以下のように変換することができます。
\(\pi(\mu, \sigma^{2} | x) \sim f(x|\mu, \sigma^{2}) \cdot \pi(\mu | \sigma^{2}) \cdot \pi(\sigma^{2}) \)
 

そして、ここで仮定から尤度と事前分布は以下のように表せることができます。
\(f(x|\mu, \sigma^{2}) = \displaystyle \frac{1}{\sqrt{2\pi\sigma^{2}}} \exp{ - \displaystyle \frac{(x - \mu)^{2}}{2\sigma^{2}}} \)
\(\pi(\mu | \sigma^{2}) = \displaystyle \frac {\sqrt{m_{0}}} {\sqrt{2\pi} \sigma} \exp\Bigr\{ - \frac{m_{0}(\mu-\mu_{0})^{2}} {2\sigma^{2}} \Bigr\} \)
\(\pi(\sigma^{2}) = \displaystyle \frac{ (\frac{\beta_{0}}{2})^{\frac{\alpha_{0}}{2}}} {\Gamma(\frac{\alpha_{0}}{2})} (\sigma^{2})^{ -(\frac{\alpha_{0}}{2} + 1) } \exp \Bigr(-\frac{\beta_{0}}{2\sigma^{2}} \Bigr) \)
 

ここで、事前分布\(\pi(\mu, \sigma^{2})\)について計算をすると、ベイズの定理から、
\begin{eqnarray}
\pi(\mu, \sigma^{2})
&=& \pi(\mu | \sigma^{2}) \cdot \pi(\sigma^{2}) \\
&=& \displaystyle \frac {\sqrt{m_{0}}} {\sqrt{2\pi} \sigma} \exp\Bigr\{ - \frac{m_{0}(\mu-\mu_{0})^{2}} {2\sigma^{2}} \Bigr\} \cdot \displaystyle \frac{ (\frac{\beta_{0}}{2})^{\frac{\alpha_{0}}{2}}} {\Gamma(\frac{\alpha_{0}}{2})} (\sigma^{2})^{ -(\frac{\alpha_{0}}{2} + 1) } \exp \Bigr(-\frac{\beta_{0}}{2\sigma^{2}} \Bigr) \\
&=&
\end{eqnarray}

事後分布を求める

従って、事後分布は以下のように計算できます。
\begin{eqnarray}
\pi(\mu, \sigma^{2} | x)
&\sim& f(x|\mu, \sigma^{2}) \cdot \pi(\mu, \sigma^{2}) \\
&=& \displaystyle \frac{1}{\sqrt{2\pi\sigma^{2}}} \exp{ - \displaystyle \frac{(x - \mu)^{2}}{2\sigma^{2}}} \cdot
\end{eqnarray}

 
そして事後分布に対して、パラメータごとの条件付き事後分布を求めます。
ギブスサンプリングを用いるために、パラメータを1つ選び、それ以外は全て条件に入れます。
今回の場合\(\mu, \sigma^{2}\)がパラメータになるので、
\(\pi(\mu | \sigma^{2}, x)\)と、\(\pi(\sigma^{2} | \mu, x)\)の2つを求めます。
 

\(\mu\)の条件付き事後分布

まず\(\sigma^{2}\)を条件にした上で、\(\mu\)の事後分布を求めてみます。
 
▼ \(\mu\)の条件付き事後分布 ▼
\begin{eqnarray}
\pi(\mu | \sigma^{2}, x)
&=& \displaystyle \frac{\pi(\mu, \sigma^{2} | x)}{\pi(\sigma^{2})} \\
&=& (\sigma^{2})^{- \frac{\alpha_{0}+1}{2}-1} \exp \Bigr\{ - \frac{ m_{1}(\mu - \mu_{1})^{2} + \beta_{1} } {2\sigma^{2}} \Bigr\} \\
&=& \frac{ (\frac{\beta_{0}}{2})^{\alpha_{0}}{2} }{\Gamma(\frac{\alpha_{0}}{2})} (\sigma^{2})^{-\frac{\alpha_{0}}{2} +1} \exp( - \frac{\beta}{2\sigma^{2}}) \\
&=& \frac{ (\sigma^{2})^{- \frac{\alpha_{0}+1}{2}-1} \exp \Bigr\{ - \frac{ m_{1}(\mu - \mu_{1})^{2} + \beta_{1} } {2\sigma^{2}} \Bigr\} }{ \frac{ (\frac{\beta_{0}}{2})^{\alpha_{0}}{2} }{\Gamma(\frac{\alpha_{0}}{2})} (\sigma^{2})^{-\frac{\alpha_{0}}{2} +1} \exp \Bigr( - \frac{\beta}{2\sigma^{2}} \Bigr) } \\
&=& (\sigma^{2})^{- \frac{\alpha_{0}+1}{2}-1} (\sigma^{2})^{-( -\frac{\alpha_{0}}{2} +1)} \cdot \exp \Bigr\{ - \frac{ m_{1}(\mu - \mu_{1})^{2} + \beta_{1} } {2\sigma^{2}} \Bigr\}
\exp \Bigr(\frac{\beta}{2\sigma^{2}} \Bigr) \\
&=& (\sigma^{2})^{-\frac{1}{2}} \exp \Bigr\{ - \frac{m_{1}(\mu - \mu_{1})^{2}}{2\sigma^{2}} - \frac{\beta}{2\sigma^{2}} + \frac{\beta}{2\sigma^{2}} \Bigr\} \\
&=& (\sigma^{2})^{-\frac{1}{2}} \exp \Bigr\{ - \frac{(\mu - \mu_{1})^{2}}{ \frac{2\sigma^{2}}{m_{1}}} \Bigr\} \\
&\sim& N \Bigr(\mu_{1}, \displaystyle \frac{\sigma^{2}}{m_{1}} \Bigr)
\end{eqnarray}

\(\sigma^{2}\)の条件付き事後分布

まず\(\mu\)を条件にした上で、\(\sigma^{2}\)の事後分布を求めてみます。
 

\(\sigma^{2}\)の条件付き事後分布を求めるには、以下の逆ガンマ分布を用います。

 
▼ \(\sigma^{2}\)の条件付き事後分布 ▼
これを用いて以下を計算すると、
\begin{eqnarray}
\pi(\sigma^{2} | \mu, x)
&=& \displaystyle \frac{\pi(\sigma^{2}, \mu | x)}{\pi(\mu)} \\
&=& \displaystyle \frac{ (\sigma^{2})^{- \frac{\alpha_{0}+1}{2} - 1} \exp \Bigr\{ -\frac{m_{1}(\mu-\mu_{1})^{2}+\beta_{1}}{2\sigma^{2}} \Bigr\}} {定数} \\
&\sim& IG \Bigr(\frac{\alpha_{0}+1}{2}, \frac{m_{1}(\mu-\mu_{1})^{2}+\beta_{1}}{2} \Bigr)
\end{eqnarray}

 
以上のことから、\(\mu\)、\(\sigma^{2}\)のそれぞれの条件付き事後分布は
\(\pi(\mu | \sigma^{2}, x) = N \Bigr(\mu_{1}, \displaystyle \frac{\sigma^{2(0)}}{m_{1}} \Bigr) \)
\(\pi(\sigma^{2} | \mu, x) = IG \Bigr(\displaystyle \frac{\alpha_{0}+1}{2}, \displaystyle \frac{m_{1}(\mu-\mu_{1})^{2}+\beta_{1}}{2} \Bigr) \)
となります。

Rで実装してみる

#install.packages("MCMCpack")

library(MCMCpack)
if(!require('ECctmc')) {
  install.packages('ECctmc')
  library('ECctmc')
}
install.packages("tidyverse")
library(ggplot2) #綺麗なグラフが描けるライブラリ

# ギブスサンプリングの初期値（この値を更新かけていく）
mu0 <- 0
sigma <- 1
alpha0 <- 0.02
beta0 <- 0.02
m0 <- 1
m1 <- m0 + 100
n <- 100

mu_list <- c()  #list()
sigma_list <- c() #list()


mu1 <- (m0*mu0+n*100)/(m0+n)
alpha1 <- alpha0 + n
beta1 <- beta0 + 100 + n*m0*(100-mu0)*(100-mu0)/(n+m0) 

# 更新式
for (i in 1:n) {
  mu_list = c(mu_list, mu)
  sigma_list <- c(sigma_list, sigma)
  mu <- rnorm(1,m=mu1,sd=sigma/m1)
  sigma <- rinvgamma(1,shape=(alpha1+1)/2,scale=(m1*(mu-mu1)*(mu-mu1)+beta1)/2)
}
print(mu_list)
print(sigma_list)


print(typeof(mu_list))
datamu = as.data.frame(mu_list)
datasigma = as.data.frame(sigma_list)
print(typeof(datamu))
plot(unlist(datamu), main = "ギブスサンプリングによる平均パラメータの推定", ylim=c(90,110))
plot(unlist(datasigma), main = "ギブスサンプリングによる分散パラメータの推定", ylim=c(60,140))

t.test(datamu)
t.test(datasigma)