MCMC導入編 - Simulation Methodsの基本 - - Data Science by R and Python

こんにちは。
中間発表などで時間をとられたので、実に3ヶ月ぶりの更新となってしまいました。
でも、嬉しいことに、、、このブログ毎日300前後のアクセスを頂いていて、
書いている本人としてはとても嬉しいです。

この記事のテーマはマルコフ連鎖モンテカルロ法です。
通称、Markov Chain Monte Carloの頭の文字を取ってMCMCです。
長くなりそうなので、2回に分けるかもしれません。

追記* : この記事は、まだMCMCまで到達できていません。次回、本格的に書きます。

導入

MCMCが本格的に使われ始めたのは、1990年代以降という比較的新しい方法です。
ただ、最近のベイズ流を用いた解析においては、必ずと言っていいほどMCMCが登場します。
そんなMCMCを今日は、説明していこうと思います。

まず、MCMCは、何かの特定の手法を表すものではありません。
Robert, and Casella(2004)によると、マルコフ連鎖と、モンテカルロ近似に基づく
分布の近似を行う手法の総称としてMCMCという表現が使われているということです。

Monte Carlo Statistical Methods (Springer Texts in Statistics)

作者: Christian Robert,George Casella
出版社/メーカー: Springer
発売日: 2005/08/24
メディア: ハードカバー
この商品を含むブログ (2件) を見る

で、何の役に立つのということですが、、、まずはベイズの推論から確認します。
ベイズ流に基づくパラメータの予測というのは、データが与えられたもとでの、パラメータの分布を求め、その分布に基づいて推定量を構成するという方法です。推定量の例としては「期待値」「分散」などがありますね。

何を持ってベイズ流と解釈するのかというと、「パラメータ」が分布に従うという部分が、最尤法とは異なっています。そこで、私たちは、パラメータの事後分布なるものを考えることになるわけです。

普段よく用いている、線形回帰モデルは、こんな感じですよね。

${ \begin{equation} Y_{i} = X_{i}^{T}\beta + \varepsilon_{i} \\ \varepsilon_{i} \sim N(0, \sigma^{2}) \end{equation} }$

ベイズ流に、書き直せば次のようになります。

${ \displaystyle Y|\beta, \sigma \sim N(X\beta, \sigma^{2} I)\\ (\beta, \sigma) \sim 何かの分布 }$

通常、最尤法を行う際には「条件付き」の記述はありません。なぜなら、パラメータは定数だからです。一方でベイズ流では、パラメータは分布なので、線形回帰モデルにおけるYの分布にはbeta,sigmaの条件が必要になるというわけです。

ただ、この問題は「単純な線形回帰モデル」ですから、betaの推定は問題なく行えます。興味のある方は、次の本を参照してください。無料で、ダウンロードできます。

The Elements of Statistical Learning: Data Mining, Inference, and Prediction, Second Edition (Springer Series in Statistics)

作者: Trevor Hastie,Robert Tibshirani,Jerome Friedman
出版社/メーカー: Springer
発売日: 2009/03/01
メディア: ハードカバー
購入: 1人クリック: 222回
この商品を含むブログ (14件) を見る

マルコフ連鎖モンテカルロ法へのモチベーション

さて、ここまではいいのです。しかしながら、一般的に「複雑なモデル」を考える場合には、上記のように計算するわけには行きません。次のようなケースを考えてみましょう。

${ \begin{equation} Y_{i} |\beta,u \sim \mbox{Bernoulli}(p)\\ (g(\mu_{i})|\beta,u) = X_{i}^{T}\beta + Z_{i}^{T}u\\ (\beta, u) \sim \mbox{some distribution} \end{equation} }$

このようなモデルは一般化線形混合モデルと呼ばれます。私が扱っている環境データの解析でも、このように複雑化したモデルを用いた解析事例は、いくつか知っています。例えば、複数校の成績の分布をつくりたいときなどは、各学校毎にばらつきの異なりをモデル化したいという要望があったりするわけです。このようなケースでは、どうしてもモデルが複雑になりがちです。特に、このような階層を持つようなモデルを「Hierarchical Model」とよんだりもします。

さて、このようなモデルの事後分布を考えたいとなった場合には、問題が生じます。それは、事後分布が陽に計算できないということです。つまり、パラメータの分布が計算によって求まらないということになるわけです。もちろん、最尤法に基づく推定も難しくなります。このような場合に、事後分布を近似する手段としてMCMCがあるわけです。

事後分布：
${ \begin{equation} (\beta, u | Y) \sim \mbox{Some Posterior Distribution} \end{equation} }$

マルコフ連鎖モンテカルロ法への導入

さて、MCMCを詳しく議論するためにはMCMCを構成する2つの道具を説明する必要があります。1つ目は「マルコフ連鎖(Markov Chain)」、2つ目は「モンテカルロ法(Monte Carlo)」です。まずは、簡単な「モンテカルロ法」から説明しましょう。

1. モンテカルロ法

モンテカルロ法は、大数の強法則に基づいて「積分」を和で近似する手法です。数学的に記述すれば次のようになります。以下では、Xが密度関数fに従う確率変数であるとします。最後の式は、大数の強法則によるものです。つまり、十分に大きなサンプルを、密度関数fから取ってくることができれば、任意のXの関数h(X)が近似できるというわけです。

$E_{f}[h(X)] = \int_{\chi}h(x)f(x)dx\\ \bar{h}_{m} = \sum_{i=1}^{m}h(x_{i})\\ \lim_{m\rightarrow\infty}\bar{h}_{m} \rightarrow E_{f}[h(X)]$

これをやってみるのは、非常に簡単なので、Rでコードを書いてみます。
今回は、正規分布の平均を近似したいと思います。
正規分布から、サンプルを発生させて、それに基づいて平均値を計算します。
その際、サンプル数の発生個数を1~10000まで変化させて収束する感じを図示してみましょう。

#====分散の収束性====
samples = c(1,5,10,30,50,100,1000,2000,5000,10000)
xx = matrix(0,length(samples),100) #期待値の保存用
for(i in 1:length(samples)){
	for(j in 1:100){
		x = rnorm(samples[i],0,5)
		xx[i,j] = mean(x)
	}
}
apply(xx,1,var)

#各サンプル数から計算された平均をBoxplotする。
boxplot(t(xx))

#平均の標本分散
21.812214409
3.649668231
2.558448360
0.606118897
0.369144894
0.209430104
0.027300966
0.013278143
0.006206799
0.001866782

f:id:tomoshige_n:20150315221604p:plain

結果を見れば、サンプル数を大きくすれば平均を近似できることがわかりました。このように、目標となる分布から、サンプルを発生させることができれば、その分布から様々な統計量を計算することができるわけです。これが、モンテカルロ法の基本的な考え方です。

しかしながら、次に問題になるのは「どうやって分布からサンプルを発生させるのか」という点です。そのための手法には、いくつかの方法があります。次は、サンプルを発生させる手法について議論します。この話のあとに、Markov Chainについて説明します。

サンプルを目標の分布から生成する

まず、目標の分布の関数F(x)の形が明確にわかっていて、かつ、逆変換がわかる場合です。この場合は、非常に簡単に発生させることができます。具体的には一様分布からの乱数を利用します。（※一様分布からの乱数発生には、いくつかの方法が用いられますが、これについては割愛することにしましょう）。この方法は逆変換法と呼ばれたりしているみたい。

逆変換法

まず、一様分布を用いた「逆変換」による乱数発生法を考えます。この方法は、言葉で説明すれば、分布関数Fは、実数空間R(定義域)から、[0,1](値域)区間への関数であるから、[0,1]区間の乱数を用いて、それをFの逆変換を用いれば、定義域の値を生成できるという考え方に基づきます。

具体的な例で確認する方が良いですね。では、指数分布からこの手法を用いて乱数を発生させます。まずは、指数分布の分布関数が、以下の形で表されることに注意しておきます。そして、このF(x)の定義域は[0, $\infty$ ), また、値域は[0,1]です。

$F(x) = 1 - exp(-\lambda x)$

これをXについて解いてみると、、、

$x = -\log(1-F(x))/\lambda$

となります。F(x)が[0,1]なので、1-F(x)も同様に[0,1]の値域を持ちます。さて、この部分を一様分布に置き直せば乱数を発生できるということです。*1

具体的に、rateが2の指数分布から、この手法を用いて乱数を発生させましょう。必要なものは、先ほどの議論から、次の2つです。
1．[0,1]区間の一様分布を発生させる関数 -> runif
2．発生させる目標分布の逆関数 -> さっきのxの式

u = runif(1000,0,1)
x = -log(u)/2
hist(x)

f:id:tomoshige_n:20150315232425p:plain

ばっちり、指数分布からデータを発生させることができました。でも、これって「Fの逆関数がわかっているということが前提」ですよね。では、次に目標となる分布Fがわかっているけど、先ほどのように逆関数を考えるのが難しい場合です。その場合、密度関数fを別の密度関数gでなんとかする、そんな方法がとられます。その代表格が、Accept Reject Methodが登場します。

accept-reject method

accept-reject methodを行うためには、目標分布fに対して、以下の性質を満たすようなgを用います。c(定数)として、
$\sup_{x}\frac{f(x)}{g(x)} = c$
つまり、目標分布fを、定数倍すれば、上から押さえ込むことができるような分布gを用いることになります。今回乱数を発生させたい標準正規分布は、次のような密度を持ちます。

$f(x) = \frac{1}{\sqrt{2\pi}}\exp{x^{2}}$

それを、定数倍したコーシー分布で押さえ込みます。実際、コーシー分布は次のような密度関数を持ちます。

$g(x) = \frac{1}{\pi(1+x^{2})}$

ここで、先ほどの式である

$\sup_{x}\frac{f(x)}{g(x)} = c$

を計算すると、

$\sup_{x}\frac{f(x)}{g(x)} = \sqrt{\frac{2\pi}{e}}$

となることがわかりますから、確かにコーシー分布の定数倍で押さえ込めることがわかります。一応、図で確認しておけば、 f:id:tomoshige_n:20150317190835p:plain

x=seq(-10,10,0.001)
plot(x,dnorm(x),xlim=c(-10,10),ylim=c(0,0.5),col=1,lty=1,type="l")
par(new=T)
plot(x,cauchy(x),xlim=c(-10,10),ylim=c(0,0.5),col=2,lty=1,type="l")
par(new=T)
plot(x,ccauchy(x),xlim=c(-10,10),ylim=c(0,0.5),col=3,lty=1,type="l")

labels=c("Normal","Cauchy","c times Cauchy")
legend("topleft", legend = labels, col = c(1,2,3),lty=c(1,1,1))

目標分布に対して、このような分布g(x)を選ぶことができれば、次のアルゴリズムによってシュミレーションが可能となります。*2

====アルゴリズム====
1．X ~ gに従う乱数を発生させる
2．U ~ U[0,1]に従う乱数を発生させる
3．U < f(X)/ c * g(X)を計算し、成立すればXをAcceptする。違えば、Rejectする。
4．1に戻る

Rで書くと次のようになります。

set.seed(102)
#accept個数
accept = 0
#アルゴリズムの繰り返し回数
iter=1000
#結果の保存用
result=rep(0,iter)
all=rep(0,iter)

#accept reject
for(i in 1:iter){
	y=rcauchy(1,0,1)
	all[i] = y
	u=runif(1)
	c = sqrt(2*pi/exp(1))
	print(dnorm(y)/(c*dcauchy(y)))
	if(u <= (dnorm(y)/(c*dcauchy(y)))){
		accept = accept + 1
		result[accept] = y
	}
}
#サンプリングされたもの
result = result[1:accept]

サンプルのヒストグラムを書くと次のようになります。うまく、正規分布っぽくなっているのがわかります。めでたしですね。
f:id:tomoshige_n:20150317192541p:plain

ただ、気をつけないといけないのは、fを覆うようなgが必要であるという点で、その比の最大値が発散しないということが必要です。また、最大値が発散しない場合でも、その値が非常に大きい場合には、accept率が低くなります。つまり、1つのサンプルを生成するために、アルゴリズムを数十回回すなんていうはめになることもあります。実際、accept率の期待値は、cの逆数で表されます（Acceptされる確率は幾何分布に従うことから計算できます）。そのため、今回のaccept率は1/c = 65.7%です。で、acceptの値は、上の結果から確かに65.7%付近にあることがわかります

さて、問題はここからなのです。次には目標となる分布が陽に書くことができない場合です。先ほどの階層モデルのパラメータの事後分布は、その典型例と言えるでしょう。特に、分母の積分どうするねん！っていう問題があります。さらに、押さえ込む分布がわからない。このような場合にどうするのかというのが次のステップです。私たちは、ここで「マルコフ連鎖」というものを利用することになります。今回は、ここまでにして、次でMCMCを説明しましょう

次の記事へ

次の記事で、本格的にマルコフ連鎖からMCMCの設計まで説明していきます。
（ここまで書くのに、とりあえず1日かかった。。。ので笑）

まとめ

シュミレーション手法を僕なりに噛み砕いて書いたつもりなのですが、、、なんとなくてもわかっていただければ大変嬉しいです。
次回、MCMCへ！（来週頭ぐらいには、できていると思います。。。笑）。
WinBUGSやStanを使わないで、自力で設計してみるところがゴールです。

参考にしてる本や資料など

・日本語だと一番まとまってるかも。。。ただ、理論的な部分はあまり書いてない
http://www.omori.e.u-tokyo.ac.jp/MCMC/mcmc.pdf

・Monte Carlo Statistical Method．これが個人的には一番まとまってる気がします。ただ、途中で一気に難化してゲロ吐きそうになります・・・

Monte Carlo Statistical Methods (Springer Texts in Statistics)

作者: Christian Robert,George Casella
出版社/メーカー: Springer
発売日: 2005/08/24
メディア: ハードカバー
この商品を含むブログ (2件) を見る

・Nummelin 1984
これ、本の画像さえ出てこないレベルの参考書です。MarkovChainについて測度論で議論していきます。普通の人にはいらない一冊♡

・事後分布からのMCMCでのサンプリングに関する論文（個人的には一番好き）←次回これメインで書きます。
Tierney : Markov Chains for Exploring Posterior Distributions

それでは、読んでくださってありがとうございました。

*1:と、直感的に書いてしまいましたが、きちんと書くのであれば確率の計算を行う必要があります。興味のある人は、やってみてください。

*2:これがなぜ、f(x)を近似できるのかを知りたいというかたは、Robert, Casella(2004)を参照していただくとよいでしょう。