統計学 学習ノート

はじめに

「この統計量の分布、どう計算しても導けない……」

ここまでの記事では、分布を式で導いてきました。ところが実際の問題では、標本の中央値の分散のように、きれいな式にならないものがすぐに出てきます。導けないと、推定も検定もできないのでしょうか？

例えば、すごろくを思い浮かべてください。サイコロ2個の目の和なら、36通りを数えれば確率が分かります。でも、マスごとに「3つ戻る」「1回休み」が書かれたすごろくで、上がるまでの回数の分布を計算するのは大変です。一方で、1万回遊んで記録をとれば、分布の形はすぐに分かります。

実は、これは統計学そのものの考え方です。私たちはこれまで、標本から母集団の値を推定してきました。その同じ道具を計算機の中に持ち込み、「解けない計算」を「推定の問題」に変えることができるのです。しかも、その誤差も理論で評価できます。

そこで登場するのが、今回の主役 モンテカルロシミュレーション です。

前回の記事では、不完全なデータの尤度を書きました。今回の記事では、乱数 の作り方と証明、モンテカルロ法の誤差と分散減少法、MCMC の理論、そして ブートストラップ について解説します。

乱数を作る：一様乱数から好きな分布へ

はじめに、乱数の作り方について解説します。

乱数 とは、ある確率分布に従う独立な数の列です。計算機が作るのは、決まった式で次々に作られる 擬似乱数 です。その代表が線形合同法です。

$$x_{n+1} = (ax_n + c) \bmod M$$

具体的に見てみましょう。$a = 5$、$c = 3$、$M = 16$、$x_0 = 7$ とすると、数列は 7, 6, 1, 8, 11, 10, 5, 12, 15, 14, 9, 0, 3, 2, 13, 4 のあと、7に戻ります。$M$ で割ると、0以上1未満の一様乱数に見える数列になります。周期は最大でも $M$ なので、実際にはメルセンヌ・ツイスタのような周期の長い方法を使います。式で決まるので、同じ初期値からは同じ乱数が出て、結果を再現できるのが利点です。

一様乱数 $U$ が手に入ったら、それを別の分布に変換します。方法は主に3つです。1つずつ見ていきましょう。

1. 逆関数法

$F$ を連続な分布関数、$U$ を $(0, 1)$ の一様乱数とします。$X = F^{-1}(U)$ は分布関数 $F$ をもちます。

証明は1行です。$F$ は増加関数なので、次が成り立ちます。

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

イメージは、分布関数のグラフの縦軸にダーツを投げて、グラフに当たったところの横軸を読むことです。密度の高いところは $F$ の傾きが急なので、たくさん当たります。

例えば、指数分布なら $F(x) = 1 - e^{-\lambda x}$ で、$X = -\log(1 - U)/\lambda$ です。$\lambda = 2$、$U = 0.3$ なら、$X = 0.178$ です。

離散分布にも使えます。累積確率の区間のどこに $U$ が入るかで値を決めます。

2. ボックス・ミュラー法

正規分布の $F^{-1}$ は式で書けません。そこで、2次元で考えます。$U_1, U_2$ が独立な一様乱数のとき、次の $Z_1, Z_2$ は独立な標準正規乱数になります。

$$Z_1 = \sqrt{-2\log U_1}\cos(2\pi U_2), \qquad Z_2 = \sqrt{-2\log U_1}\sin(2\pi U_2)$$

なぜでしょうか？ 独立な標準正規変数の同時密度 $\frac{1}{2\pi}e^{-(z_1^2 + z_2^2)/2}$ を極座標に変換します。ヤコビアンは $r$ なので、密度は次のようになります。

$$\frac{1}{2\pi}\, re^{-r^2/2}$$

$R$ と $\Theta$ は独立で、$\Theta$ は $(0, 2\pi)$ の一様分布です。$R^2/2$ は平均1の指数分布なので、$R = \sqrt{-2\log U_1}$ と作れます。逆関数法の組み合わせでした。

3. 棄却法

密度 $f$ の乱数が直接作れないとします。乱数を作りやすい密度 $g$ と定数 $M$ で、$f(x) \le Mg(x)$ がすべての $x$ で成り立つものを用意します。

【棄却法の手順】

1. $g$ から候補 $Y$ を引き、一様乱数 $V$ を引く
2. $V \le f(Y)/\{Mg(Y)\}$ なら $Y$ を採択する
3. 棄却したら1に戻る

長方形の的にダーツを投げて、曲線の下に刺さったものだけを残すイメージです。採択された $Y$ が密度 $f$ に従うことを、証明しましょう。

$$P(Y \le x \mid \text{採択}) = \frac{\int_{-\infty}^{x} g(y)\dfrac{f(y)}{Mg(y)}\, dy}{\int_{-\infty}^{\infty} g(y)\dfrac{f(y)}{Mg(y)}\, dy} = \frac{F(x)/M}{1/M} = F(x)$$

ここで、分母から採択確率が $1/M$ であると分かります。採択までの試行回数は幾何分布に従い、平均は $M$ です。したがって、$M$ はできるだけ小さくします。

例題：棄却法で乱数を作る

密度 $f(x) = 6x(1 - x)$（$0 < x < 1$）の乱数を、一様分布 $g(x) = 1$ を使う棄却法で作る。〔1〕最も小さい $M$ を求めよ。〔2〕採択の条件を書け。〔3〕採択確率と、1個の乱数を得るまでの平均試行回数を求めよ。

【解答】

〔1〕$f$ は $x = 1/2$ で最大値 $3/2$ をとるので、$M = 3/2$ である。

〔2〕$V \le f(Y)/M = 4Y(1 - Y)$ のとき採択する。

〔3〕採択確率は $1/M = 2/3$、平均試行回数は $3/2$ 回である。

試験のポイント：逆関数法と棄却法は、答案で正しさを証明できるようにします。棄却法の採択確率は $1/M$ です。

以上が、乱数の作り方についての解説です。

モンテカルロシミュレーション：誤差は √n で減る

次に、モンテカルロシミュレーションについて説明します。

モンテカルロシミュレーション は、乱数を使って期待値や確率を推定する方法の総称です。求めたい値が $\theta = E\{h(X)\}$ と書けるとき、$X_1, \dots, X_n$ を乱数で作り、次の平均で推定します。

$$\hat{\theta}_n = \frac{1}{n}\sum_{i=1}^{n} h(X_i)$$

これは、ただの標本平均です。したがって、これまでの理論がすべて使えます。不偏で、大数の法則から一致性があります。中心極限定理から、誤差の大きさも分かります。

$$\hat{\theta}_n \approx N\left(\theta, \frac{\sigma_h^2}{n}\right), \qquad \sigma_h^2 = V\{h(X)\}$$

要するに、誤差は $1/\sqrt{n}$ で減ります。精度を1桁上げるには、100倍の乱数が必要です。ただし、この誤差の式に次元が入っていないことがポイントです。格子点で数値積分をすると、10次元なら1軸10点でも $10^{10}$ 点が必要ですが、モンテカルロ法の誤差は次元で悪化しません。味噌汁の味見は、鍋の大きさによらず一口で足りるのと同じです。

具体的に見てみましょう。$\theta = \int_0^1 e^x\, dx = e - 1 = 1.718$ を推定します。$U$ を一様乱数として、$h(U) = e^{U}$ の平均をとります。分散は次のとおりです。

$$\sigma_h^2 = E(e^{2U}) - \theta^2 = \frac{e^2 - 1}{2} - (e - 1)^2 = 0.242$$

誤差を95%の確率で ±0.01 に収めたいなら、$n = (1.96 \times 0.492/0.01)^2 \approx 9300$ です。

せっかくなので、統計の問題の例も見てみましょう。モンテカルロシミュレーションは、モデルが崩れたときの手法の性能を調べるのに使えます。例えば、指数分布から10個の標本をとり、正規分布を仮定した両側5%の $t$ 検定を、正しい平均を $H_0$ として20万回くり返します。すると、棄却の割合は約0.10でした。名目の5%の約2倍です。分布の非対称性が効いているのです。

シミュレーションの結果そのものにも誤差があります。割合 $p$ を $B$ 回の試行で見積もるときの標準誤差は $\sqrt{p(1-p)/B}$ で、$p = 0.05$、$B = 10000$ なら約0.0022です。報告するときは、この誤差も添えます。

ここからは、同じ乱数の数で誤差を小さくする 分散減少法 を2つ見ましょう。

1. 対称変数法

$U$ を使ったら、$1 - U$ も使います。$(h(U) + h(1 - U))/2$ の分散は次の式です。

$$V\left\{\frac{h(U) + h(1 - U)}{2}\right\} = \frac{\sigma_h^2 + \operatorname{Cov}\{h(U), h(1 - U)\}}{2}$$

$h$ が単調なら、共分散は負です。大きい値と小さい値が組で現れて、打ち消し合うのです。上の例では $\operatorname{Cov} = e - (e - 1)^2 = -0.234$ なので、ペアの平均の分散は0.0039です。同じ2個の乱数を独立に使うと分散は0.121なので、約31倍の効率です。

2. 重点サンプリング

$X$ が密度 $f$ に従うとき、別の密度 $q$ から乱数を引いて、次のように重みをつけます。

$$E_f\{h(X)\} = \int h(x)\frac{f(x)}{q(x)}\, q(x)\, dx = E_q\left\{h(X)\frac{f(X)}{q(X)}\right\}$$

狙いは、まれな事象のあたりから重点的に乱数を引くことです。釣りで、魚のいる場所に集中して糸を垂らしてから、そこだけ釣ったことの偏りを重みで直すイメージです。分散を最も小さくする $q$ は、$|h(x)|f(x)$ に比例する密度です。

例題：まれな事象の確率を推定する

$p = P(Z > 3) = 0.00135$（$Z$ は標準正規変数）を推定する。〔1〕素朴なモンテカルロ法で、相対誤差（標準誤差 ÷ $p$）を10%にする $n$ を求めよ。〔2〕$N(3, 1)$ を提案分布とする重点サンプリングの重みを書け。〔3〕〔2〕の推定量の、1個あたりの分散は $6.17 \times 10^{-6}$ である。相対誤差を10%にする $n$ を求めよ。

【解答】

〔1〕推定量の分散は $p(1 - p)/n$ だから、$\sqrt{(1 - p)/(np)} = 0.1$ を解いて、$n = (1 - p)/(0.01p) \approx 74000$ である。

〔2〕重みは $\phi(x)/\phi(x - 3) = \exp(-3x + 9/2)$ である。推定量は、$x > 3$ のときこの重みをとり、それ以外で0をとる変数の平均である。

〔3〕$n = 6.17 \times 10^{-6}/(0.01 p^2) \approx 340$ である。

同じ精度を得るのに、必要な乱数が約220分の1になりました。

試験のポイント：モンテカルロ法の誤差は $\sigma_h/\sqrt{n}$ です。必要な $n$ を求める問題は、信頼区間の幅からサンプルサイズを決める問題と同じ手順で解けます。

以上が、モンテカルロシミュレーションについての説明です。

MCMC：独立に引けないなら連鎖で歩く

続いて、MCMCの理論について解説します。

目標の分布 $\pi$ が高次元で、正規化定数も分からないと、独立な乱数はまず作れません。MCMC（マルコフ連鎖モンテカルロ法）は、定常分布が $\pi$ になるマルコフ連鎖を作り、その軌跡を乱数の代わりに使う方法です。準1級で学んだMHアルゴリズムがその代表です。

正しさの証明は2段階です。推移確率を $P(x, y)$ とします。第1段は、詳細釣り合い条件 $\pi(x)P(x, y) = \pi(y)P(y, x)$ から定常性が出ることです。$x$ で足すと、次のようになります。

$$\sum_x \pi(x)P(x, y) = \sum_x \pi(y)P(y, x) = \pi(y)$$

第2段は、MHの推移確率がこの条件を満たすことです。$x \ne y$ のとき、$P(x, y) = q(y \mid x)\,\alpha(x, y)$ で、採択確率 $\alpha(x, y) = \min\{1, \pi(y)q(x \mid y)/(\pi(x)q(y \mid x))\}$ を代入します。

$$\pi(x)P(x, y) = \min\{\pi(x)q(y \mid x),\ \pi(y)q(x \mid y)\}$$

右辺は $x$ と $y$ の入れ替えで変わらないので、条件が成り立ちます。あとは、連鎖が既約で周期をもたなければ、定常分布への収束と、軌跡の平均が期待値に収束するエルゴード定理が保証されます。

ただし、連鎖の値は独立ではありません。平均の分散を計算しましょう。定常な連鎖の自己相関を $\rho_k$、$h(X)$ の分散を $\sigma^2$ とします。

$$V(\bar{h}_n) = \frac{\sigma^2}{n}\left\{1 + 2\sum_{k=1}^{n-1}\left(1 - \frac{k}{n}\right)\rho_k\right\} \approx \frac{\sigma^2}{n}\left(1 + 2\sum_{k=1}^{\infty}\rho_k\right)$$

右の括弧を 自己相関時間 と呼びます。独立な乱数に比べて分散が何倍になるかを表す数です。$n$ をこの数で割った値を 有効標本サイズ と呼びます。$\rho_k = \rho^k$ なら、括弧の中は $(1 + \rho)/(1 - \rho)$ です。

例えば、$\rho = 0.9$ なら19倍なので、1万回の連鎖は独立な乱数の約526個分の価値しかありません。毎日同じ人にだけ聞き込みをしても、新しい情報があまり増えないのと同じです。

例題：2状態のMHアルゴリズム

状態 $\{1, 2\}$ の上の分布 $\pi = (1/3, 2/3)$ を目標とする。提案は、確率 $1/2$ で今の状態を、確率 $1/2$ でもう片方の状態を選ぶ。〔1〕推移確率行列を求めよ。〔2〕詳細釣り合い条件と $\pi P = \pi$ を確かめよ。〔3〕第2固有値は $1/4$ である。状態の指示関数の $\rho_k$ を求め、独立な乱数に比べて平均の分散が何倍になるかを求めよ。

【解答】

〔1〕1から2への提案の採択確率は $\min(1, 2) = 1$、2から1へは $\min(1, 1/2) = 1/2$ である。よって、$P(1, 2) = 1/2$、$P(2, 1) = 1/4$ で、次のようになる。

$$P = \begin{pmatrix} 1/2 & 1/2 \\ 1/4 & 3/4 \end{pmatrix}$$

〔2〕$\pi(1)P(1, 2) = 1/6 = \pi(2)P(2, 1)$ である。また、$\pi P = (1/6 + 1/6,\ 1/6 + 1/2) = (1/3, 2/3)$ である。

〔3〕2状態の連鎖では $\rho_k = (1/4)^k$ である。倍率は $(1 + 1/4)/(1 - 1/4) = 5/3$ である。

実際に100万回動かすと、状態2の割合は約0.666、1次の自己相関は約0.25となり、理論と一致しました。

試験のポイント：棄却されたときは、現在の値をもう1回記録します。これが $P(x, x)$ の項で、MCMCの分散の計算にも、自己相関の形で効いてきます。

以上が、MCMCについての解説です。

ブートストラップ：標本を母集団の代わりにする

最後に、ブートストラップを見ていきましょう。

統計量 $\hat{\theta}$ の標準誤差を知りたいとします。理想は、母集団から標本をくり返し取ることです。しかし、手元には1組の標本しかありません。

そこで、ブートストラップ（bootstrap）の出番です。標本 $x_1, \dots, x_n$ の各値に確率 $1/n$ を置いた 経験分布 $\hat{F}_n$ を、母集団の分布 $F$ の代わりに使います。名前は「自分の靴ひもを引っ張って自分を持ち上げる」という無理な話に由来します。とはいえ、$\hat{F}_n$ は $F$ に近いので、この無理が通るのです。箱の中に標本のコピーを入れて、そこから復元抽出するイメージです。

【ブートストラップの手順】

1. 標本から $n$ 個を復元抽出して、$x_1^{*}, \dots, x_n^{*}$ を作る
2. その統計量 $\hat{\theta}^{*}$ を計算する
3. 1と2を $B$ 回くり返す
4. $\hat{\theta}^{*}$ の標準偏差を標準誤差の推定値に、$\bar{\hat{\theta}}^{*} - \hat{\theta}$ を偏りの推定値にする

信頼区間は、$\hat{\theta}^{*}$ の2.5%点と97.5%点を使う パーセンタイル法 が基本です。

ここで2種類の近似を区別しましょう。$F$ を $\hat{F}_n$ で置き換える近似は、$n$ で決まります。$B$ 回で打ち切る近似は、モンテカルロ法の誤差で、$B$ を増やせば消せます。全部で $n^n$ 通りの再標本があり、並びを無視しても $\binom{2n-1}{n}$ 通りあります。$n = 10$ なら、後者でも92,378通りです。

$\hat{\theta}$ が標本平均なら、全部を列挙しなくても理論で解けます。$\hat{F}_n$ のもとで $X^{*}$ の平均は $\bar{x}$、分散は $\hat{\sigma}^2 = \frac{1}{n}\sum(x_i - \bar{x})^2$ なので、次が成り立ちます。

$$V^{*}(\bar{X}^{*}) = \frac{\hat{\sigma}^2}{n} = \frac{n - 1}{n} \cdot \frac{s^2}{n}$$

いつもの $s^2/n$ の $(n-1)/n$ 倍です。ブートストラップは、平均に対しては、ほぼ教科書どおりの答えを返すわけです。

再標本の中身の性質も確かめましょう。特定の値が1回も選ばれない確率は $(1 - 1/n)^n$ なので、$n$ が大きいと $e^{-1} = 0.368$ に近づきます。つまり、1回の再標本には、もとの値の約63.2%が現れます。

例題：3個のデータの完全ブートストラップ

標本は 2, 4, 9 である。〔1〕順序を区別した再標本の総数と、並びを無視した総数を求めよ。〔2〕標本平均のブートストラップ分散を求めよ。〔3〕標本中央値 $\hat{\theta}^{*}$ のブートストラップ分布を求め、標準誤差と偏りの推定値を求めよ。

【解答】

〔1〕$3^3 = 27$ 通りと、$\binom{5}{3} = 10$ 通りである。

〔2〕$\bar{x} = 5$、$\hat{\sigma}^2 = (9 + 1 + 16)/3 = 26/3$ だから、$V^{*} = 26/9 = 2.89$ である。

〔3〕中央値が2になるのは、3回のうち2回以上が2のときで、確率は $3 \times (1/3)^2(2/3) + (1/3)^3 = 7/27$ である。中央値が9になる確率も同じで、4になる確率は $13/27$ である。

$$E^{*}(\hat{\theta}^{*}) = \frac{2 \times 7 + 4 \times 13 + 9 \times 7}{27} = 4.78, \qquad V^{*}(\hat{\theta}^{*}) = \frac{803}{27} - 4.78^2 = 6.91$$

標準誤差の推定値は $\sqrt{6.91} = 2.63$ である。標本中央値は4なので、偏りの推定値は $4.78 - 4 = 0.78$ である。

中央値の分散は、理論で導くのが難しい量でした。それでも、手順そのものはいつも同じです。

一方で、ブートストラップは万能ではありません。代表的な失敗は、$U(0, \theta)$ の最大値です。本当の最大値の分布は連続です。しかし、再標本の最大値は、確率 $1 - (1 - 1/n)^n \to 0.632$ で標本の最大値にぴったり一致します。分布の端を決める統計量では、$\hat{F}_n$ が端の情報を持てないのです。

補足（ここは読み飛ばしてOK）： 経験分布の代わりに、推定した母数のモデル $F_{\hat{\theta}}$ から乱数を引く方法を、パラメトリック・ブートストラップと呼びます。モデルが正しければ精度は上がりますが、モデルの誤りをそのまま引き継ぎます。

試験のポイント：「$\hat{F}_n$ を $F$ の代わりにする」が本質です。標本平均のブートストラップ分散の分母が $n - 1$ でなく $n$ の分散になることも、答案で示せるようにします。

以上が、ブートストラップについての解説です。

細かい手順を忘れてしまっても問題ありません。ここでは、「解けない計算を、乱数を使った推定の問題に変える」という点だけ押さえておきましょう。

今回のまとめ

それでは、今回の内容を振り返ります。

今回は乱数を使って分布を作る方法を紹介しました。一様乱数を好きな分布に変換し、平均をとれば期待値が求まり、誤差は標本平均の理論で評価できます。独立に引けないならマルコフ連鎖で歩き、母集団がなければ標本を箱に入れ直します。どれも、解けない問題を「くじを引いて推定する問題」に変える工夫でした。

今回のポイントを振り返りましょう。

振り返り

• 逆関数法は $F^{-1}(U)$、棄却法は採択確率 $1/M$ で、どちらも数行で証明できる。
• モンテカルロ法の誤差は $\sigma_h/\sqrt{n}$ で、次元によらない。
• 対称変数法と重点サンプリングで、同じ乱数の数でも分散を減らせる。
• MCMCは詳細釣り合い条件で定常分布を保証し、自己相関の分だけ有効標本サイズが減る。
• ブートストラップは経験分布を母集団の代わりにする。分布の端の統計量では失敗する。

キーワードのおさらい

• 乱数: ある確率分布に従う独立な数の列。計算機では擬似乱数を使う
• モンテカルロシミュレーション: 乱数を使って期待値や確率を推定する方法
• MCMC: 目標の分布を定常分布にもつマルコフ連鎖で乱数を作る方法
• ブートストラップ: 標本からの復元抽出で統計量の分布を推定する方法

さて、今回のMCMCは、もともとある分野の計算を支えるために大きく広まりました。それが、母数そのものに確率分布を与えるベイズ法です。

次回は、統計数理の最後の記事として、事後分布の計算について紹介します。

それでは、次の記事へ進んでいきましょう。