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

数理統計　第11章 §11.1　★ 目標分布を定常分布に持つマルコフ連鎖を作り、その標本平均で期待値を近似する。正しさの根拠は詳細つり合い

第6章から第10章までの議論では、最尤推定量は尤度方程式を手で解いて閉じた式で書け、事後分布の平均は共役事前分布のおかげで積分の公式から求まった。実際のデータ解析では、この2つはどちらも成り立たないことが多い。

尤度方程式が $\lambda/(1-\mathrm{e}^{-\lambda})=\bar x$ のような超越方程式になれば、解は数値的に求めるしかない。事後分布が共役でなければ、事後平均の分母に現れる規格化定数 $\int L(\theta)\pi(\theta)\,\mathrm{d}\theta$ が計算できない。この章では、こうした場合に使う計算の方法を3つ扱う。

1つめは、分布から乱数を発生させ、その平均で積分を近似するモンテカルロ法である。規格化定数が分からない分布からでも乱数を作れるようにする道具がマルコフ連鎖モンテカルロ法（MCMC）で、§11.1 で扱う。2つめは、推定量の分布そのものが手で求まらないとき、データから作った経験分布で母集団を置き換えて標準誤差や偏りを見積もるブートストラップ法で、§11.2 で扱う。3つめは、尤度方程式を反復で解く Newton–Raphson 法とEM アルゴリズムで、§11.3 で扱う。

1級の統計数理で計算機を使う問題は出ない。この章の内容が問われるとすれば、紙の上で導ける部分である。具体的には、棄却法で採択される乱数の密度が目標の密度になることの証明と採択確率、Metropolis–Hastings 法の推移核が詳細つり合いを満たすことの確認、ブートストラップ分布の平均と分散を経験分布から厳密に計算すること、Newton 法と EM アルゴリズムの更新式を尤度から導き、不動点が最尤推定量であることと収束の速さを示すこと、の4種類である。どれも第6章の尤度と第4章の条件付き分布の計算の応用なので、この章ではその計算を省かずに書く。

モンテカルロ積分

密度 $f$ を持つ確率変数 $X$ について、次の期待値を求めたいとする。

$$\mu_h=E[h(X)]=\int h(x)f(x)\,\mathrm{d} x$$

ここで $h$ は与えられた関数である。$f$ に従う独立な乱数 $X_1,\dots,X_N$ が作れるなら、

$$\hat\mu_h=\frac1N\sum_{t=1}^N h(X_t)\tag{11.1}$$

は $\mu_h$ の不偏推定量である。大数の法則により、$N\to\infty$ でこれは $\mu_h$ に確率収束する。分散 $\sigma_h^2=V[h(X)]$ が有限なら

$$V[\hat\mu_h]=\frac{\sigma_h^2}{N},\qquad \sqrt N(\hat\mu_h-\mu_h)\ \text{は近似的に}\ \mathrm N(0,\sigma_h^2)\ \text{に従う}$$

である（後半は中心極限定理による）。$\sigma_h^2$ 自体は $h(X_t)$ の標本分散で推定できる。そのため、近似の誤差の大きさが計算と同時に分かる。これをモンテカルロ積分と呼ぶ。

誤差が $1/\sqrt N$ でしか減らないことは不利に見えるが、この速さは積分の次元によらない。格子点の上で和をとる数値積分は、次元 $d$ が上がると同じ精度を得るための格子点の数が指数的に増える。事後分布の母数が数十次元になるベイズ推測でモンテカルロ法を使う理由はここにある。

確率 $p=P(X\in A)$ を求めるときは、$h=\bm 1_A$（指示関数）にとる。このとき $h(X_t)$ は成功確率 $p$ の Bernoulli 変数である。したがって $\sigma_h^2=p(1-p)$ で、相対誤差（標準誤差と $p$ の比）は

$$\frac{\sqrt{p(1-p)/N}}{p}=\sqrt{\frac{1-p}{Np}}$$

になる。$p$ が小さいほど、同じ相対誤差を得るために大きな $N$ が必要になる。

$f$ から直接乱数を作りにくいときは、別の密度 $g$ を使う手がある。$g$ は $f>0$ のところで正であるとする。$g$ から乱数 $Y_1,\dots,Y_N$ を作り、

$$\tilde\mu_h=\frac1N\sum_{t=1}^N h(Y_t)\frac{f(Y_t)}{g(Y_t)}\tag{11.2}$$

で近似してもよい。$E_g[h(Y)f(Y)/g(Y)]=\int h(y)f(y)\,\mathrm{d} y=\mu_h$ なので不偏である。これを重点サンプリングと呼び、比 $f/g$ を重みと呼ぶ。重みの分散が大きいと近似が悪くなるので、$g$ は $\lvert h\rvert f$ に形が近いものを選ぶ。

乱数の生成：逆関数法と棄却法

計算機は区間 $(0,1)$ 上の一様分布 $\mathrm U(0,1)$ に従う乱数を作れるとする。ほかの分布の乱数はこれを変換して作る。

証明

$F$ は右連続で非減少である。このことから、次の同値が成り立つ。

$$F^{-1}(u)\le x\iff u\le F(x)$$

• ($\Rightarrow$) $F^{-1}(u)\le x$ とする。下限の定義と右連続性から $F(F^{-1}(u))\ge u$ である。$F$ は非減少なので $F(x)\ge F(F^{-1}(u))\ge u$ となる。
• ($\Leftarrow$) $u\le F(x)$ とする。このとき $x$ は集合 $\{x:F(x)\ge u\}$ に属する。下限の定義から $F^{-1}(u)\le x$ である。

したがって

$$P(X\le x)=P(F^{-1}(U)\le x)=P(U\le F(x))=F(x)$$

である。

逆関数法は $F^{-1}$ が式で書けるとき（指数分布・Weibull 分布・Cauchy 分布など）に使う。正規分布やガンマ分布のように $F^{-1}$ が式で書けないときは、次の棄却法を使う。

証明

1回の試行で「$Y\le x$ かつ採択」となる確率を求める。$Y=y$ を与えたとき採択される条件付き確率は $f(y)/\{Mg(y)\}$ で、これは 0 以上 1 以下である。したがって

$$P\Bigl(Y\le x,\ U\le\frac{f(Y)}{Mg(Y)}\Bigr)=\int_{-\infty}^x\frac{f(y)}{Mg(y)}\,g(y)\,\mathrm{d} y=\frac1M\int_{-\infty}^xf(y)\,\mathrm{d} y .\tag{11.3}$$

$x\to\infty$ とすると採択確率 $P(\text{採択})=1/M$ を得る。採択された $Y$ の分布関数は、条件付き確率として

$$P(Y\le x\mid\text{採択})=\frac{P(Y\le x,\ \text{採択})}{P(\text{採択})}=\int_{-\infty}^xf(y)\,\mathrm{d} y$$

となる。これは $f$ の分布関数である。各試行は独立で毎回確率 $1/M$ で採択されるので、試行回数は幾何分布に従う。

$f\le Mg$ の両辺を積分すると $1\le M$ を得る。$M$ が1に近いほど棄却が少ない。最良の $M$ は $\sup_x f(x)/g(x)$ である。

$f$ の規格化定数が分からないときも、手続きはそのまま使える。$f=c\tilde f$ と書いて、$\tilde f$ だけが分かるとする。$\tilde f\le\tilde Mg$ となる $\tilde M$ を使うと、採択の判定は

$$U\le\frac{\tilde f(Y)}{\tilde Mg(Y)}$$

となり、$c$ が現れない。採択確率は $1/(c\tilde M)$ になる。

マルコフ連鎖と定常分布

ベイズ推測の事後分布 $\pi(\theta\mid\bm x)\propto L(\theta)\pi(\theta)$ のように、密度が定数倍を除いてしか分からず、しかも高次元の分布から乱数を作りたいとする。棄却法は次元が上がると $M$ が指数的に大きくなり、ほとんどの提案が棄却される。

そこで、独立な乱数を作ることはあきらめ、目標分布を定常分布に持つマルコフ連鎖を走らせて、その軌道の平均で期待値を近似する。これがマルコフ連鎖モンテカルロ法である。マルコフ連鎖の一般論は第12章 §12.5 で扱う。ここでは MCMC に要る定義と結果だけを書く。

状態空間 $\mathcal X$ 上の確率変数の列 $X_0,X_1,X_2,\dots$ を考える。どの時刻 $t$ についても

$$P(X_{t+1}\in A\mid X_t=x,X_{t-1},\dots,X_0)=P(X_{t+1}\in A\mid X_t=x)$$

が成り立ち、右辺が $t$ によらないとき、この列を（時間的に一様な）マルコフ連鎖と呼ぶ。右辺を $(x,A)$ の関数とみたものを推移核と呼び、$K(x,A)$ と書く。$K(x,\cdot)$ が密度 $k(x,y)$ を持つときは、その密度も推移核と呼ぶ。状態が有限個なら推移核は推移確率行列 $P=(p_{ij})$、$p_{ij}=P(X_{t+1}=j\mid X_t=i)$ である。

証明

(11.5) の両辺を $x$ について積分すると、次を得る。

$$\int\pi(x)k(x,y)\,\mathrm{d} x=\pi(y)\int k(y,x)\,\mathrm{d} x=\pi(y)$$

最後の等号は、$k(y,\cdot)$ が確率密度で積分が 1 になることによる。推移核が点 $y$ に確率の塊（とどまる確率）を持つ場合も、和と積分を合わせて同じ計算ができる。

定常条件 (11.4) は、$y$ ごとに $x$ についての積分を含むので確かめにくい。詳細つり合い (11.5) は2点 $x,y$ ごとの等式なので確かめやすい。MCMC の推移核は、詳細つり合いを満たすように設計する。

定常分布を持つだけでは、連鎖がそこへ近づく保証は無い。次の定理（証明は第12章 §12.5 に譲る）が、軌道の平均を使ってよい根拠である。

Metropolis–Hastings 法

採択確率 (11.6) には $\pi$ が比 $\pi(y)/\pi(x)$ でしか現れないので、規格化定数は約分されて要らない。ベイズ推測で使える理由はこれである。

証明

$x$ にとどまる部分は $x=y$ のときだけ現れ、そこでは (11.5) の両辺は同じ式なので自明に等しい。$x\neq y$ を示せばよい。$\pi(y)q(y,x)\le\pi(x)q(x,y)$ の場合を考える。このとき

$$\alpha(x,y)=\frac{\pi(y)q(y,x)}{\pi(x)q(x,y)},\qquad \alpha(y,x)=1$$

なので

$$\pi(x)k(x,y)=\pi(x)q(x,y)\frac{\pi(y)q(y,x)}{\pi(x)q(x,y)}=\pi(y)q(y,x)=\pi(y)q(y,x)\alpha(y,x)=\pi(y)k(y,x).$$

逆向きの不等号の場合は $x$ と $y$ の役割を入れ替えれば同じ計算になる。2つの場合をまとめると$\pi(x)k(x,y)=\min\{\pi(x)q(x,y),\pi(y)q(y,x)\}$ で、右辺は $x,y$ について対称である。命題 11.4 から $\pi$ は定常分布になる。

提案密度の選び方で、よく使う形が2つある。

1. ランダムウォーク Metropolis 法：提案が対称で $q(x,y)=q(y,x)$ となるもの。例えば $Y=x+\varepsilon$、$\varepsilon\sim\mathrm N(0,s^2)$ とする。採択確率は

$$\alpha(x,y)=\min\{1,\pi(y)/\pi(x)\}$$

になる。密度が上がる候補は必ず採択し、下がる候補は比の確率で採択する。2. 独立型 MH 法：$q(x,y)=g(y)$ のように、現在の点によらない提案を使うもの。$w=\pi/g$ とおくと、採択確率は $\min\{1,w(y)/w(x)\}$ になる。

Gibbs サンプラー

多次元の分布 $\pi(x_1,\dots,x_d)$ を考える。$x_j$ 以外の成分をまとめて $x_{-j}$ と書く。各成分の完全条件付き分布 $\pi(x_j\mid x_{-j})$ からは乱数を作れるとする。Gibbs サンプラーは、成分を1つずつ順に完全条件付き分布から引き直す方法である。2次元なら、現在の $(x_1^{(t)},x_2^{(t)})$ から

$$x_1^{(t+1)}\sim\pi(x_1\mid x_2^{(t)}),\qquad x_2^{(t+1)}\sim\pi(x_2\mid x_1^{(t+1)})$$

の順に更新する。

証明

$(X_1',X_2)$ の同時密度は、$X_2$ の周辺密度と条件付き密度の積で

$$\pi_2(x_2)\,\pi(x_1'\mid x_2)=\pi(x_1',x_2)$$

である。

2つの成分の更新はどちらも $\pi$ を保つので、それを続けて行う1回の Gibbs の推移も $\pi$ を保つ。また、1つの成分の更新は、提案を完全条件付き分布にとったMH 法で採択確率が常に 1 になるもの、とみることもできる（確認テスト 11.1-4）。

確認テスト

オリジナル

1. 11.1-1 $U\sim\mathrm U(0,1)$ とする。$X=-\lambda^{-1}\log U$ が $\mathrm{Ex}(\lambda)$ に従うことを示せ。逆関数法の $F^{-1}(U)$ と比べて何が違うか説明せよ。
2. 11.1-2 密度 $f(x)=6x(1-x)\ (0<x<1)$ の $\mathrm{Beta}(2,2)$ 乱数を、棄却法で作る。提案は $\mathrm U(0,1)$ とする。最小の $M$、採択の条件、採択確率、1つの乱数を得るまでの試行回数の平均と分散を求めよ。
3. 11.1-3 状態 1, 2 の2状態マルコフ連鎖の推移確率行列を

$$P=\begin{pmatrix}1-a&a\\ b&1-b\end{pmatrix}\qquad(0<a,b<1)$$

とする。定常分布 $\pi$ を求め、詳細つり合いを満たすことを確かめよ。さらに $P^n$ を求めよ。各行が $\pi$ に近づく速さが $\lvert 1-a-b\rvert ^n$ であることを示せ。4. 11.1-4 2次元の目標分布 $\pi(x_1,x_2)$ に対する MH 法で、第2成分を固定したまま第1成分だけを動かすものを考える。提案密度は

$$q((x_1,x_2),(y_1,x_2))=\pi(y_1\mid x_2)$$

とする。採択確率 (11.6) が常に 1 になることを示せ。5. 11.1-5 $Z\sim\mathrm N(0,1)$ とし、$p=P(Z>2)$ を推定する。推定量は、独立な $N$ 個の正規乱数の指示関数の平均とする。推定量の相対誤差（標準誤差と $p$ の比）を0.1 以下にするのに必要な最小の $N$ を求めよ。$p=0.02275$ を使ってよい。