諸注意
- 問題本文は公式サイトまたは公式問題集を参照してください
- 統計検定2級の資格を持つ方を前提に解説していきます
問題13-1
メトロポリス・ヘイスティング法を用いて、式1の目標分布（混合正規分布）から10000個の乱数を発生させることを考える。この時、メトロポリス・ヘイスティング法はStep1からStep6までの手順に従い、乱数を発生させるものとする。
- Step1: 初期値設定
⇨ 初期値\(x_0\)と定数a(\(>0\))を定めて、カウンタtを0に設定する。 - Step2: 乱数発生
⇨ 一様分布U(-a, a)から\(\epsilon\)を生成し、\(y=x_t+\epsilon\)とする。 - Step3: 採択確率
⇨ U(0, 1)からuを生成し、\(x_{t+1} = \begin{cases} y & u \leq \alpha(x_t, y) \\ x_t & \text{それ以外} \end{cases}\)とする。 - Step4: カウンタアップ
⇨ カウンタtを1増加させる。 - Step5: 出力制御
⇨ カウンタtが1000以下の場合は出力せず、1000より大きい場合は出力する。 - Step6: 終了条件
⇨ カウンタtが11000なら終了。そうでなければStep2に戻る。
メトロポリス・ヘイスティング法では、目標分布の確率密度関数が\(\pi(x)\)の時、採択確率\(\alpha(x_t, y)\)は式2の通り表される。\(\phi(\cdot)\)を標準正規分布の確率密度関数とした時、Step3の採択確率を\(\phi(\cdot)\)を用いて表せ。
式1（目標分布）
\begin{align}
&\frac{1}{4} N(0, 1) + \frac{3}{4} N(6, 1) \\
\end{align}
式2（採択確率）
\begin{align}
&\alpha(x_t, y) = min\left(1, \frac{\pi(y)}{\pi(x_t)}\right) \\
\end{align}
答え \(\alpha(x_t, y) = min\left(1, \frac{\frac{1}{4} \phi(y) + \frac{3}{4} \phi(y-6)}{\frac{1}{4} \phi(x_t) + \frac{3}{4} \phi(x_t-6)}\right)\)
解説
標準正規分布の確率密度関数
問題文より、標準正規分布の確率密度関数を\(\phi(\cdot)\)と表します。
この表現は統計の教科書などでよく見られるもので、正規分布の確率密度関数を\(\phi(x; \mu, \sigma^2)\)と表します。標準正規分布は\(\mu=0, \sigma^2=1\)になり、\(\phi(x; 0, 1)\)を\(\phi(x)\)と表記することが多いです。
目標分布
目標分布の確率密度関数(式1)は、\(\phi(\cdot)\)を用いて以下の通り表現できます。
式3
\begin{align}
\pi(x) = \frac{1}{4} \phi(x; 0, 1) + \frac{3}{4} \phi(x; 6, 1) \\
\end{align}
この時、第1項は標準正規分布に従いますが、第2項は平均が6, 分散が1の正規分布に従います。そのため、第2項のみ平均を0にする処理を追加します。
式4
\begin{align}
\pi(x) &= \frac{1}{4} \phi(x; 0, 1) + \frac{3}{4} \phi(x-6; 0, 1) = \frac{1}{4} \phi(x) + \frac{3}{4} \phi(x-6) \\
\end{align}
採択確率に代入
式4の結果を採択確率(式2)に代入します。
式5
\begin{align}
&\alpha(x_t, y) = min\left(1, \frac{\frac{1}{4} \phi(y) + \frac{3}{4} \phi(y-6)}{\frac{1}{4} \phi(x_t) + \frac{3}{4} \phi(x_t-6)}\right) \\
\end{align}
したがって、正解は式5の通りになります。
問題13-2
Step1で初期値を\(x_0=6\), aを0.1, 1.0, 6.0のいずれかに設定した時に得られた10000個の乱数のヒストグラムと時系列プロットの組み合わせとして適切なものを答えよ。
（グラフを用意できなかったため、特徴を解説します）
答え （解説参照）
解説
前提知識
メトロポリス・ヘイスティング法
マルコフ連鎖モンテカルロ（MCMC）手法の一種。直接サンプリングすることが難しい複雑な目標分布からランダムサンプルを生成するために使用します。
目的分布を山で例えると、現在の地点\(x_t\)よりも新しい地点\(x_{t+1}\)の方が高ければ移動します。そうでなければ、移動するかその場に留まるかを採択確率で決めます。この時、現在の地点よりも新しい地点が低いほど、採択確率も低くなります。
定数aの役割
本問では定数aが0.1, 1.0, 6.0のいずれかに設定されており、10000個の乱数がどのように生成されるかを問われています。
\(x_{t+1}\)が取り得る値は\(y(=x+t+\epsilon)\)か\(x_t\)自身であり、\(\epsilon\)はU(-a, a)に従います。これは\(x_{t+1}\)は(\(x_t\pm a\))の範囲内に必ず含まれることを意味します。
ここから定数aが小さいほど、時系列プロットの振れ幅が小さくなることが分かります。
したがって、aが(0.1, 1.0, 6.0)に設定した時のグラフは、時系列プロットの振れ幅が小さい(ウ), (ア), (イ)の順に対応します。
問題13-3
Step5にあるように、\(x_1, …, x_{1000}\)を出力に加えない理由を答えよ。
答え バーンイン期間を考慮したため
解説
前提知識
バーンイン期間
メトロポリス・ヘイスティングス法などの乱数生成サンプリングでは、初期値の影響を排除し、目標分布に収束させるためにある程度のサンプルを破棄することが一般的です。
バーンイン期間
初期値は適当に定めることも多く、初期値の影響を受けるサンプリングの初期をバーンイン期間と呼び、生成開始直後の乱数は破棄することが一般的です。
したがって、\(x_1, …, x_{1000}\)を出力に加えない理由はバーンイン期間を考慮するためです。
