マルコフ連鎖モンテカルロ法のお勉強

目次 (2023/4/27 ver.0.02)

1番目のタブが選択された！

1.はじめに

もくじ

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

1.はじめに
  • ベイズ推論の枠組みで、事後分布の推論などに用いるマルコフ連鎖モンテカルロ法（通称：MCMC）の概要について、簡単に説明してみたいと思います。
  • 対象は、最初の「ベイズ推論におけるMCMCの必要性」を読んで、理解いただける方になると思います。
2.ベイズ推論におけるMCMCの必要性
  • データをX、パラメータ（潜在変数）をθと書くこととします。 確率モデル（尤度モデル）を「$p(X|\theta)$」、事前分布を「$p(\theta)$」とすると、事後分布は以下のように計算できます。

    • $\Large \displaystyle p(\theta|X)=\frac{p(X|\theta)p(\theta)}{p(X)}$

  • ここで、重要となるのが「$p(X)$」の存在です。
  • 「$p(\theta|X)$」が確率分布であるためには以下の条件が必要となります。

    • $\Large \displaystyle \int p(\theta|X)d\theta=\int \frac{p(X|\theta)p(\theta)}{p(X)}d\theta=1$

  • 従って、

    • $\Large \displaystyle p(X)=\int p(X|\theta)p(\theta)d\theta$

  • となります。
  • 要するに1番最初の式で、事後分布を解析的に求める場合、上の積分計算によって「$p(X)$」を求めなければなりません。
  • しかし、一般にはこの積分計算は解析的に実行できません。そこで、登場するのが、MCMCです。
  • MCMCを端的に説明すると、「解析的に事後分布「$p(\theta|X)$」を計算せずとも、近似的に事後分布からサンプリングする方法」と言えます。
  • 事後分布からサンプリングできれば、以下のモンテカルロ法を用いて、平均・分散などの各種統計量を求めることができます。
  • 以降では、モンテカルロ法について説明した後に、「近似的な事後分布からサンプリング」を行う具体的な方法を説明していきます。

1.2 モンテカルロ法とは

1.モンテカルロ法とは
  • モンテカルロ法は、乱数生成による数値計算手法です。
  • ここでは、簡単な期待値計算の例について考えてみましょう。
  • 確率分布「$p(x)$」の期待値、すなわち確率変数「$x$」の厳密な平均値は以下の値で定義されます。

    • $\Large \displaystyle E[X]=\int xp(x)dx$

  • 定義から明らかなように、この値は積分計算によって求めますが、以下のように近似的に求めることもできます。

    • $\Large \displaystyle E[X] \approx \frac{1}{N}\sum_{i=1}^{N}x_i$

  • ここで、「$x_i$」 は「$p(x)$」からの独立なサンプルです。
  • 右辺の値は、確率変数「$X$」の標本平均であることを考えると、期待値と近い値であることが期待されます。
  • なお「$N\rightarrow\infty$」において、大数の法則から標本平均は期待値に（確率）収束します。
  • このように、厳密な積分によってではなく、確率分布からの乱数（サンプリング）によって期待値計算を行う手法がモンテカルロ法と呼ばれる手法です。（本来はもっと広義です）

1.3 マルコフ連鎖とは

1.マルコフ連鎖とは
  • マルコフ連鎖とは、確率変数列「$\{X_t\}_{t=0}^{\infty}$」 に対して、以下で定義される概念です。

    • $\Large \displaystyle p(x_t|x_{t-1},x_{t-2},\cdots)=p(x_t|x_{t-1})$

  • 要するに、「$t$」期における確率変数「$X_t$」の条件付き分布が、一期前の確率変数「$X_{t-1}$」の値にのみ依存するとき、 確率変数列「$\{X_t\}_{t=0^{\infty}}$」 はマルコフ連鎖であると言います。
  • マルコフ連鎖である確率変数列「$\{X_t\}_{t=0}^{\infty}$」の同時分布は、

    • $\Large \displaystyle p(x_0,\cdots,x_t)=p(x_0)p(x_1|x_0)\times\cdots\times p(x_t|x_{t-1},\times,x_0)$

    $\Large \displaystyle =p(x_0)\prod_{i=1}^{t}p(x_i|x_{i-1})$

  • と書き直すことができます。
  • ここで、「$p(x_i|x_{i-1})$」を遷移核と呼び、「$K(x_{i-1},x_i)$」と書くこととします。
  • 遷移核「$K(x_{i-1},x_i)$」は、「$(i-1)$期に$x_{i-1}$にいて、$i$期に$x_i$に移動する確率」と解釈できます。
  • すなわち、遷移核はマルコフ連鎖である確率変数列の移動ルールであると言えます。
  • スタート地点である「$x_0$」と移動ルール「$K(x_{i-1},x_i)$」を決めてしまえば、確率変数列「$\{X_t\}_{t=0}^{\infty}$」の挙動は完全に決定してしまいます。
  • このことから、遷移核「$K(x_{i-1},x_i)$」がマルコフ連鎖において極めて重要な要素であることがおわかりいただけると思います。
  • もう一度確認しておくと、上式に従う「確率変数列$\{X_t\}_{t=0}^{\infty}$」がマルコフ連鎖であり、その中でも移動ルールを記述している上式の「確率分布」が遷移核に対応しています。
  • 用語の定義を混同しないように気をつけてください。
2.マルコフ連鎖の具体例
  • マルコフ連鎖の具体例として、ARモデルが挙げられます。
  • 最もシンプルなARモデルであるAR(1)モデルは、

    • $\Large \displaystyle p(x_i|x_{i-1})=N(\rho x_{i-1},\sigma^2)$

  • と書くことができます。
  • 確率変数「$X_i$」の条件付き分布が一期前の値「$x_{i-1}$」にのみ依存しているので、この確率分布に従う確率変数列はマルコフ連鎖であり、 このマルコフ連鎖の遷移核（移動ルール）は、上式そのものです。
  • 下図は、「$(\rho,\sigma^2)=(0.5,1)$」とした場合の確率変数列「$\{X_t \}_{t=0}^{\infty}$」のサンプルをプロットしたものになります。
3.定常分布
  • MCMCを考える上で、重要な概念が定常分布です。
  • 定義から厳密に述べると少し難しいので、ここでは直感的な理解に留めたいと思います。
  • 下図は、先ほどの図にサンプルした確率変数列のヒストグラムを添えたものです。
  • 定常分布は、このヒストグラムに対応したような概念で、確率変数列を横から見たときに個々の確率変数「$X_i$」が従っている確率分布であると見なすことができます。
  • もう少し丁寧に述べると、「定常分布$P$を持つマルコフ連鎖に従う確率変数列$\{X_t\}_{t=0}^{\infty}$の各要素である確率変数「$X_t$」は、 定常分布「$P$」に従う」と言えます。
  • 要するに、

    • $\Large \displaystyle X_t\sim P$

  • のようなイメージで考えることができます。（厳密には少し異なるので、あくまでイメージとして認識してください）
    • 4.事後分布の推論
    • 最初に述べたように、MCMCを用いる目的は、「近似的に事後分布からサンプリングする」ことでした。
    • ここで仮に、定常分布として事後分布を持つマルコフ連鎖があったとすると、 このマルコフ連鎖の各要素「$X_t$」は定常分布である事後分布に従う確率変数であると見なすことができる、 すなわちマルコフ連鎖の各要素を近似的に事後分布からのサンプルと見ることができるわけです。
    • ここまでの話をまとめると、定常分布として事後分布を持つマルコフ連鎖を構築できれば、 マルコフ連鎖の各要素「$X_t$」は事後分布からのサンプルとみなせるという話でした。
    • また、マルコフ連鎖は遷移核によって定められるものであったことを思い出すと、 「定常分布として事後分布を持つようなマルコフ連鎖を構築する」というのは、 「定常分布として事後分布を持つマルコフ連鎖となるように、遷移核を設計する」ことに他なりません。
    • サンプリングしたい分布に対して適切な遷移核を設計するのが、MCMCの役割となります。
    5.メトロポリス-ヘイスティング・アルゴリズム(M-H)
    • MCMCには、様々な種類がありますが、ここではメトロポリス-ヘイスティング・アルゴリズム（以下：M-H）を具体例として示したいと思います。
    • Step 1.
      • 「$t=1$」として初期値「$\theta^{(0)}$」を設定する。
    • Step 2.
      • 「$\theta^{(t)}$」の候補「$\tilde{\theta}$」をマルコフ連鎖から生成する。

        • $\Large \displaystyle \tilde{\theta}\leftarrow q(\theta^{(t-1)},\theta)$

    • Step 3.
      • 「$\tilde{\theta}$」の採択確率を以下の式で計算する。
    • Step 4.
      • 「$\theta^{(t)}$」を以下のルールに従い更新する。（「$U(0,1)$」は一様分布）
    • Step 5.
      • tを1つ増やしてStep 2.に戻る

    以上

    2.PythonによるMCMC

    本節は、「 Pythonによるマルコフ連鎖モンテカルロ法(MCMC) 」を引用している。

    もくじ

    2.1 モンテカルロ法

    1.円の面積を求める
    • モンテカルロ法についてまずは見ていきます。
    • これは乱数を用いた数値計算となっています。
    • 例えば、円の面積を求めるといったものが代表的な例です。
    • また、モンテカルロ法では一般的に、ランダムに打つ点の個数を多くするほど精度が良くなっていくことが知られています。
    • 下記の実装では円の面積を求める計算をモンテカルロ法を用いて行っています。
    • この結果を見ると今回の場合の正解値(0.7854)に収束してい様子がわかります。 
      • # coding: UTF-8
import numpy as np
import matplotlib.pyplot as plt

#print(np.arange(3,10,2))
#forループで乱数の数を変えて計算(start=1000,stop=500000+1,step=1000)
N_monte_list=np.arange(1000,500000+1,1000)

#結果を格納するリストを作成
sampling_rate_list=[]

for N_monte in N_monte_list:
	x_data=np.random.rand(N_monte)
	y_data=np.random.rand(N_monte)
	print("N_monte",N_monte) #1000～5000000
	#print("x_data", x_data) #0～1の乱数(x用)
	#print("y_data", y_data) #0～1の乱数(y用)

	sampling_rate=np.sum((x_data**2+y_data**2)<1)/N_monte
	print("sampling_rate", sampling_rate)

	sampling_rate_list.append(sampling_rate)

plt.plot(N_monte_list,sampling_rate_list)
plt.xlabel('N_monte')
#plt.savefig('./figure/fig_montecalro.jpg')
plt.show()

    2.2 棄却サンプリング

    1.事後分布からモンテカルロ法にてサンプリングを行う
    • 上記の例では、円の面積を求めるということを行いましたが、ベイズ統計では、事後分布からモンテカルロ法にてサンプリングを行います。
    • そのための手法として棄却サンプリングがあります。
    • 例として事後分布がbeta分布に従うとして考えていきます。
    • Beta分布とは下記のような分布のことです。

      • $\Large \displaystyle \beta(x|\alpha,\beta)=\frac{x^{\alpha-1}(1-x)^{\beta-1}}{\int_0^1[t^{\alpha-1}(1-t)^{\beta-1}]dt}$ 

      #Beta分布
X=np.arange(0,1,0.02)
dist_1=beta.pdf(X,1.5,2)

plt.plot(X,dist_1,label='α=1.5,β=2')
plt.legend()
plt.xlabel('x')
plt.ylabel('Probabirity')
#plt.savefig('./figure/fig_beta.jpg')
plt.show()
    • 今回行うこととしてはモンテカルロ法を用いてBeta分布のサンプリングすることである。
    • 棄却分布ではBeta分布p(x)(目標分布)に従う乱数の発生は諦め、サンプルが簡単な一様分布q(x)(提案分布)を用意し、これが、下記を満たすように定数Mを調整するということを行います。

      • $\Large \displaystyle f(x)\leq Mg(x)$

    • つまりMg(x)はf(x)を覆うような一様分布です。
    • よってMの値はBeta分布の最大値の接するようにg(x)が与えられるように決めます。
    • 次に棄却サンプリングのアルゴリズムについてみていきます。
    • これは下記の5つのステップのようになります。
      1. g(x)の定義域から乱数xをサンプリング
      2. Mg(x)地域から乱数rをサンプリングし、これが下記を満たすかを判定

        3. $\Large \displaystyle r(x)\leq p(x)$

      3. 真の場合はそのxのカウントを1増やし、偽の場合はそのサンプルを棄却する。
      4. 2と3をN回繰り返す。
      5. 各xでカウント数のヒストグラムを作成
    • このようにすることで各xに対して乱数を生成した場合にBeta分布以下の範囲内に収まる確率を得ることができ、これが事後確率に比例するということを用いています。
    • これをpythonにて実装していきます。
    • まずはライブラリのインポート、グラフのスタイル設定、乱数のセード固定を行います。 
      • # coding: UTF-8
#ライブラリインポート
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import beta,uniform
from scipy import optimize

#グラフの見た目変更
plt.style.use("ggplot")

#乱数のSeedを固定
np.random.seed(999)
    • 次にベータ分布の最大値を取る点の探索を行います。
    • ここで最小化問題に帰着させるために、最小化を行うメソッド(fmin)にBeta分布に-1を掛けたものを渡しています。 
      • #Beta分布のパラメータ設定
a,b=1.5,2
#x軸を作成
x=np.linspace(beta.ppf(0.001,a,b),beta.ppf(0.999,a,b),1000)
#Beta関数の密度関数を定義
p=beta(a,b).pdf
#最大化された結果(x座標)を取得
res=optimize.fmin(lambda x:-p(x),0.3)
#最大化された結果(y座標)を取得
y_max=p(res)
    • これで最大値(今回の場合1.443)が求まったので、これでMg(x)までが求まったことになります。
    • よって、この値を用いて、棄却サンプリングを行います。
    • まずはサンプリングを行います。 
      • #サンプリング数
NMCS=50000

#x軸方向のサンプリング（0～1の区間でサンプリング）
x_mcs=uniform.rvs(size=NMCS)

#y軸方向のサンプリング（0～y_maxの区間でサンプリング）
r=uniform.rvs(size=NMCS)*y_max
    • 次にこれらを用いて、棄却するかの判定を行います。
    • これを実行し結果を表示すると下記のようになります。 
      • #Beta分布以下かどうかの条件判定
accept=x_mcs[r<=p(x_mcs)]

#結果の表示
plt.hist(accept,density=True,bins=30,rwidth=0.8,label='rejection sampling')
plt.plot(x,beta.pdf(x,a,b),label='Target')
plt.legend()
#plt.savefig('./figure/fig_rejection_sampling.jpg')
plt.show()

      Optimization terminated successfully.
         Current function value: -1.443376
         Iterations: 10
         Function evaluations: 20
    • この結果、事後分布のサンプリングを行うことができました。

    以上

    3.自己相関

    もくじ

    3.1 t分布

    1.t分布とは

    以上

    4.メトロポリス法

    もくじ

    4.1 t分布とは

    1.t分布とは
    • 正規分布は、母平均μと母分散σ2のみで示されたわけですが、一般にこれらの値は知られていない場合が多いようです。
    • そこで、母平均の代わりに標本平均で代用します。
    • ここでは、この標本平均がどれほど母平均に近いかが問題となるわけです。
    • そこで、この標本平均x~はどのような分布を示すのでしょうか？
    • 「$N(\mu,\sigma^2)$」の正規分布を示す母集団から、n 個の標本を取り出しその標本平均を「$\tilde{x}$」とおくと、 「$\tilde{x}$」は「$N(\mu,\frac{\sigma^2}{n})$」の正規分布を示す。
    • つまり、母集団の分布が下図のような正規分布「$N(\mu,\sigma^2$」だとすると、
    • 標本平均の分布は下図のような正規分布「$N(\mu,\frac{\sigma^2}{n})$」となります。
    • つまり、「$\frac{(\tilde{x}-\mu)}{\sqrt{\frac{\sigma^2}{n}}}=(\tilde{x}-\mu)\frac{\sqrt{n}}{\sigma}$」は正規分布「$N(0,1)$」を示します。
    • このように、標本平均が正規分布を示すには標本数が非常に多いということが前提となります。

    4.2 t過程回帰モデル

    1.t過程回帰モデルとは
    • ガウス過程と同様に、確率過程「$f_{TP}$」が、入力の任意の有限集合「$X$」に対して、 関数値のベクトル「$f_{TP}$」の分布が多変量「$t$」分布であるとき、t過程という。
    • 「$f_{TP}$」の分布は、「$t$」分布であるので、自由度「$\nu\in R^{+}\\[0, 2] $」 ν∈R+[0, 2] と平均関数 µ および、共分散 (カーネル) 関数を要素とする共分散行列 K を指定するによって完全に特定できる。
    • これを、「$f_{TP} \sim T(\nu,\mu(x),K(x,x_0))$」と書く。

    以上

    5.Hamilton Monte Carlo(HMC)法

    「Scilab使用例(4) -グラフ描画の基本-」から引用。

    もくじ

    5.1 グラフ描画の基礎

    1.ウィンドウを1行2列に分割してグラフ描画
    • グラフを作成する為には描画するデータが必要ですが、要素は行方向, 列方向どちらに並べておいても(描画する際に転置すれば良いので)問題ありません。 
      • //ウィンドウを1行2列に分割した左側部分を選択
subplot(1,2,1);

x=[0:0.1:2*%pi]';
plot2d(sin(x));

//右側部分を選択
subplot(1,2,2);

x=[0:0.1:2*%pi]';
plot2d(x,sin(x));
    2.ウィンドウを1行4列に分割してグラフ描画
    • plot2d：指定した座標を直線で繋いでいく
    • plot2d2：次の点になるまで同じy座標となる階段状のカーブを描く
    • plot2d3：各点をx軸からの縦線で表現
    • plot2d4：plot2dと同じだが、点を矢印で繋いでいく
    • という風に異なるコマンドを用いることにより見せ方の異なるグラフの描画も可能です。
    • その他にも棒グラフ, 円グラフ, 散布図や等高線図など様々な種類のグラフが用意されていますので、順番に掘り下げた記事をあげていきたいと思います。 
      • clf()
x=[0:0.1:2*%pi]';
subplot(1,4,1); plot2d(x,sin(x));
subplot(1,4,2); plot2d2(x,sin(x));
subplot(1,4,3); plot2d3(x,sin(x));
subplot(1,4,4); plot2d4(x,sin(x));

    以上