統計学 学習ノート

はじめに

数えられなくても、何回も振って割合を見れば確率が分かる

「公式が分からないなら、実際に何回もやってみればいいのでは？」

サイコロを3個振って、目の合計が15以上になる確率はいくつでしょうか？ 数え上げれば、216通りのうち20通りなので、約0.093です。では、サイコロが30個で、条件がもっと複雑ならどうでしょう。数え上げは、すぐに手に負えなくなります。

ここで、コンピュータに1万回サイコロを振ってもらい、起きた割合を数える手があります。3個の場合なら、結果は0.09前後になります。理屈はとても素朴ですが、これが現代の統計学の計算を支えています。

実は、前回の記事のMCMCも、この発想の上に立っていました。分布の性質が式で分からなくても、その分布に従う数をたくさん作れば調べられるのです。

そこで登場するのが、今回の主役 モンテカルロ法 と ブートストラップ法 です。前回の記事では、事後分布から乱数を引くMCMCを学びました。今回の記事では、乱数の作り方、採択・棄却法、モンテカルロ法、そして推定量の精度を手元のデータだけで見積もるジャックナイフ法とブートストラップ法について解説します。

乱数：コンピュータでサイコロを振る

線形合同法の乱数は計算で作るので、必ず繰り返す

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

乱数 とは、でたらめに並んでいて、次の値を予測できない数の列のことです。ただし、コンピュータは決まった計算しかできません。そこで、計算で作るのに乱数のように見える 疑似乱数 を使います。

最も単純な作り方は、線形合同法 です。

$$X_{n+1} = (aX_n + c) \bmod m$$

$\bmod m$ は、$m$ で割った余りをとるという意味です。$X_n/m$ が0以上1未満の一様乱数になります。最初の値 $X_0$ を シード（種）と呼びます。シードが同じなら、同じ乱数列が再現できます。

具体的に見てみましょう。$a = 5$、$c = 3$、$m = 16$、シードを7とすると、次のようになります。

7 → 6 → 1 → 8 → 11 → 10 → 5 → 12 → 15 → 14 → 9 → 0 → 3 → 2 → 13 → 4 → 7 → ……

16個で1周して、7に戻りました。余りの種類は $m$ 個しかないので、必ず繰り返します。この繰り返しの長さを周期と呼びます。実用の乱数はメルセンヌ・ツイスタのように周期がとても長い方法を使いますが、「計算で作っている」という点は同じです。

さて、一様乱数 $U$ が作れれば、他の分布の乱数に変換できます。その基本が 逆関数法 です。分布関数を $F$ とすると、次のように作れます。

$$X = F^{-1}(U)$$

なぜでしょうか？ $P(X \le x) = P(F^{-1}(U) \le x) = P(U \le F(x)) = F(x)$ となり、$X$ の分布関数が $F$ に一致するからです。変数変換の記事で学んだ考え方の、裏返しです。

例えば、レート $\lambda$ の指数分布なら $F(x) = 1 - e^{-\lambda x}$ なので、$X = -\log(1 - U)/\lambda$ です。$\lambda = 2$ で $U = 0.5$ なら $X \approx 0.347$、$U = 0.8$ なら $X \approx 0.805$ です。

せっかくなので、別の例も見てみましょう。離散分布でも考え方は同じです。ガチャの結果が、確率0.2で金、0.5で銀、0.3で銅だとします。

【離散分布の乱数の作り方】

• $U < 0.2$ なら金
• $0.2 \le U < 0.7$ なら銀
• $0.7 \le U$ なら銅

要するに、0から1の物差しを確率の長さで区切り、$U$ が落ちた区間を結果とするだけです。

正規分布の $F^{-1}$ は式で書けません。そこで、独立な一様乱数 $U_1, U_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)$$

試験のポイント：逆関数法は、$F(x) = u$ を $x$ について解くだけです。指数分布は頻出です。

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

採択・棄却法：的に当たった分だけ採用する

山の下に刺さった矢だけを採用すると、f に従う

次に、逆関数が使えないときの方法を説明します。

ベータ分布のように、$F^{-1}$ が簡単に書けない分布はたくさんあります。そこで使うのが、採択・棄却法（acceptance-rejection method）です。作りやすい分布から候補を出し、一部を捨てることで、目標の分布に形を整えます。

ダーツに例えてみましょう。密度関数 $f$ のグラフを壁に描き、その山を囲む長方形に、ダーツをでたらめに投げます。山の下に刺さったものだけを採用すれば、採用した矢の横の位置は、山の高いところに多く集まります。つまり、$f$ に従います。

式で書くと、次のようになります。候補の分布の密度を $g$ とし、すべての $x$ で $f(x) \le M g(x)$ となる定数 $M$ を用意します。

【採択・棄却法の手順】

1. $g$ から候補 $X$ を、一様分布から $U$ を引く
2. $U \le \dfrac{f(X)}{M g(X)}$ なら $X$ を採択する
3. 棄却したら1に戻る

$M g$ が、$f$ の上にかぶせた屋根です。正しさを確かめましょう。候補が採択される確率は、次のとおりです。

$$P(\text{採択}) = \int g(t)\, \frac{f(t)}{M g(t)}\, dt = \frac{1}{M}$$

さらに、採択された $X$ の分布を計算します。

$$P(X \le x \mid \text{採択}) = \frac{\int_{-\infty}^{x} g(t)\, \dfrac{f(t)}{M g(t)}\, dt}{1/M} = \int_{-\infty}^{x} f(t)\, dt = F(x)$$

たしかに $f$ に従います。また、1個を得るのに必要な候補の数は、平均 $M$ 個です。屋根が低いほど無駄が減るので、$M$ はできるだけ小さくとります。

なお、$f$ は定数倍が分からなくても使えます。屋根を $f$ の定数倍に合わせればよいからです。

例題：ベータ分布の乱数を作る

$\text{Beta}(2, 2)$ の密度は $f(x) = 6x(1 - x)$（$0 < x < 1$）です。$g$ を $(0, 1)$ の一様分布として、採択・棄却法を使います。$M$ の最小値、候補 $x = 0.2$ の採択確率、全体の採択率を求めましょう。

【解き方】

1. $M$ の最小値は $f$ の最大値： $f(0.5) = 1.5$ より $M = 1.5$
2. 候補の採択確率： $f(0.2)/1.5 = 0.96 \div 1.5 = 0.64$
3. 全体の採択率： $1/M = 2/3$

答えは、$M = 1.5$、採択確率0.64、全体の採択率は約0.67です。つまり、3個の候補から平均2個が残ります。

以上が、採択・棄却法についての説明です。

モンテカルロ法：平均で積分を解く

味見の平均は、誤差1/√Nの速さで真の値に近づく

続いて、乱数で計算する方法を解説します。

モンテカルロ法 は、乱数を使って期待値や積分を近似する方法の総称です。名前は、カジノで有名な地名からとられました。

考え方は、大数の法則そのものです。$X_1, \dots, X_N$ が密度 $f$ に従う乱数なら、次のように近似できます。

$$E[h(X)] = \int h(x)\, f(x)\, dx \ \approx\ \frac{1}{N}\sum_{i=1}^{N} h(X_i)$$

鍋のスープを全部飲まなくても、よくかき混ぜて何口か味見をすれば、味が分かるのと同じです。

確率も期待値の一種です。$h$ を「条件を満たせば1、そうでなければ0」とすれば、平均は確率になります。冒頭のサイコロの例がこれです。

具体的に見てみましょう。積分 $\int_0^1 e^{-x^2} dx$ は、ふつうの方法では解けません。$U$ を $(0, 1)$ の一様分布とすると、これは $E[e^{-U^2}]$ なので、一様乱数を $N$ 個作って $e^{-U_i^2}$ を平均すれば求まります。真の値は約0.7468です。

精度はどれくらいでしょうか？ 中心極限定理から、推定値の標準誤差は次のようになります。

$$\text{SE} = \frac{\sigma}{\sqrt{N}}, \qquad \sigma^2 = V[h(X)]$$

今の例では $\sigma \approx 0.201$ なので、$N = 10000$ で標準誤差は約0.002です。

ここがポイントです。誤差は $1/\sqrt{N}$ でしか減りません。精度を10倍にするには、乱数を100倍に増やす必要があります。その代わり、この速さは積分の次元によりません。格子点で数値積分すると、次元が増えるごとに点の数が爆発しますが、モンテカルロ法は高次元でも使えます。

せっかくなので、別の例も見てみましょう。正方形 $[0, 1]^2$ に点を $N$ 個落とし、半径1の四分円に入る割合 $\hat{p}$ を数えます。四分円の面積は $\pi/4$ なので、$4\hat{p}$ で $\pi$ を推定できます。$N = 10000$ なら標準誤差は $4\sqrt{(\pi/4)(1 - \pi/4)/N} \approx 0.016$ で、小数第1位がやっと安定する程度です。

補足（ここは読み飛ばしてOK）： 同じ $N$ で精度を上げる工夫を、分散減少法と呼びます。代表例は重点サンプリングで、$\int h f\, dx = E_g[h(X) f(X)/g(X)]$ と書き直し、重要な領域によく当たる $g$ から乱数を引きます。確率の小さい事象の計算で特に有効です。

試験のポイント：モンテカルロ法の標準誤差は $\sigma/\sqrt{N}$ です。比率を推定するなら $\sigma^2 = p(1 - p)$ です。

以上が、モンテカルロ法についての解説です。

ジャックナイフ法：1人ずつ抜いて確かめる

1人ずつ抜いて、結論がどれくらい変わるかを見る

さらに、推定量の偏りと精度を、データだけで調べる方法を説明します。

ここからは問題が変わります。手元の標本から推定量 $\hat{\theta}$ を計算したとき、その標準誤差を知りたいとします。$\bar{X}$ なら $s/\sqrt{n}$ という公式があります。でも、中央値や相関係数のように、公式が複雑な推定量もあります。母集団からもう一度標本をとれれば分かりますが、それは無理です。

そこで、手元の標本から標本を作り直します。これが 再標本化（リサンプリング）の考え方です。

ジャックナイフ法（jackknife）は、標本から1個ずつ取り除いて推定量を計算し直す方法です。会議で1人ずつ席を外してもらい、結論がどれくらい変わるかを見るイメージです。誰が抜けても結論が同じなら、その結論は安定しています。

$i$ 番目を除いた $n - 1$ 個で計算した推定量を $\hat{\theta}_{(i)}$、その平均を $\bar{\theta}$ とします。

$$\widehat{\text{bias}} = (n - 1)(\bar{\theta} - \hat{\theta}), \qquad \text{SE}_{\text{jack}} = \sqrt{\frac{n - 1}{n}\sum_{i=1}^{n} \left(\hat{\theta}_{(i)} - \bar{\theta}\right)^2}$$

偏りを修正した推定量は、$\hat{\theta} - \widehat{\text{bias}} = n\hat{\theta} - (n - 1)\bar{\theta}$ です。

なぜ $n - 1$ が掛かるのでしょうか？ 偏りが $E[\hat{\theta}] = \theta + a/n + \cdots$ の形だとします。すると、次のように $1/n$ の項が消えます。

$$E\left[n\hat{\theta} - (n - 1)\bar{\theta}\right] \approx n\left(\theta + \frac{a}{n}\right) - (n - 1)\left(\theta + \frac{a}{n - 1}\right) = \theta$$

標準誤差の係数も同じです。$\hat{\theta}_{(i)}$ どうしはほとんど同じデータで計算しているので、差がとても小さく、それを $(n-1)$ の分だけ膨らませて補正します。例えば、$\hat{\theta} = \bar{x}$ なら $\hat{\theta}_{(i)} - \bar{\theta} = (\bar{x} - x_i)/(n - 1)$ なので、公式に入れるとちょうど $s/\sqrt{n}$ に一致します。

具体的に見てみましょう。データは 3, 5, 6, 8, 13 の5個です。

【標本平均の標準誤差】

1. $\bar{x} = 7$。1個を除いた平均は 8, 7.5, 7.25, 6.75, 5.5
2. これらの平均は7で、$\text{SE}_{\text{jack}} = \sqrt{0.8 \times 3.625} \approx 1.70$
3. 公式の $s/\sqrt{n} = \sqrt{14.5/5} \approx 1.70$ と一致する

せっかくなので、偏りのある推定量でも試してみましょう。$n$ で割る分散は $\hat{\theta} = 11.6$ です。1個を除いて同じ計算をすると 9.5, 13.25, 14.1875, 14.1875, 3.25 で、平均は10.875です。

$$\widehat{\text{bias}} = 4 \times (10.875 - 11.6) = -2.9, \qquad \text{修正後} = 11.6 + 2.9 = 14.5$$

修正後の値は、不偏分散の14.5に一致しました。

一方で、ジャックナイフ法には弱点があります。データの変化に対して滑らかでない推定量には向きません。例えば、中央値は1個を抜いてもごく少数の値しかとらないので、標準誤差の推定が不安定になります。

以上が、ジャックナイフ法についての説明です。

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

復元抽出を繰り返すと、中央値の揺れと信頼区間が見える

最後に、再標本化の主役を見ていきましょう。

ブートストラップ法（bootstrap）は、手元の標本を母集団の代わりと考え、そこから復元抽出で同じ大きさの標本を何度も作る方法です。エフロンが考案しました。名前は、「自分のブーツのつまみを引っ張って、自分を持ち上げる」という言い回しからきています。

くじ引きに例えてみましょう。$n$ 個のデータをそれぞれくじに書いて箱に入れます。1枚引いて記録し、箱に戻す。これを $n$ 回繰り返すと、1組のブートストラップ標本ができます。箱に戻すので、同じデータが何度も入ったり、1度も入らなかったりします。

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

1. 元の標本から、復元抽出で大きさ $n$ の標本を作る
2. その標本で推定量 $\hat{\theta}^{*}$ を計算する
3. 1と2を $B$ 回繰り返す（$B$ は1,000回以上が多い）
4. $\hat{\theta}^{*}_1, \dots, \hat{\theta}^{*}_B$ の標準偏差を標準誤差、平均と $\hat{\theta}$ の差を偏りの推定値とする

要するに、「母集団 → 標本」という関係を、「標本 → ブートストラップ標本」という関係でまねるのです。母集団の分布の代わりに、各データに確率 $1/n$ を置いた 経験分布 を使っています。

具体的に見てみましょう。10人の通学時間が 12, 15, 9, 20, 14, 11, 30, 13, 16, 10（分）だったとします。中央値は13.5分です。中央値の標準誤差には簡単な公式がありません。ブートストラップ法で20万回計算すると、標準誤差は約1.6分でした。

信頼区間も作れます。$\hat{\theta}^{*}$ を小さい順に並べ、下から2.5%と97.5%の値をとったものが パーセンタイル法 の95%信頼区間です。この例では、約 (11, 17.5) 分です。正規分布を仮定していないので、左右が非対称な区間になっています。

例題：ブートストラップ標本の確率を計算する

データ 3, 5, 6, 8, 13 からブートストラップ標本を作ります。次を求めましょう。

(1) 値13が1度も選ばれない確率　(2) ブートストラップ標本の中央値が3になる確率

【解き方】

1. (1) 毎回 $4/5$ の確率で13以外が選ばれるので、$(4/5)^5 \approx 0.328$
2. (2) 中央値が3になるのは、5回のうち3回以上で3が選ばれるとき
3. 選ばれる回数は $B(5, 1/5)$ に従うので、$P(\text{3回以上}) = 0.0512 + 0.0064 + 0.0003 \approx 0.058$

答えは (1) 約0.328、(2) 約0.058です。

一般に、あるデータが選ばれない確率は $(1 - 1/n)^n$ で、$n$ が大きいと $e^{-1} \approx 0.368$ に近づきます。つまり、ブートストラップ標本には、元のデータの約63.2%が含まれます。

ブートストラップ法にも限界があります。

【ブートストラップ法の注意点】

• 標本が小さく母集団を代表していないと、結果も偏る
• 最大値のように分布の端を見る統計量はうまくいかない。元の最大値が約63%の確率で再び選ばれるため
• 時系列のような従属するデータでは、順番が壊れる。連続した塊ごとに抽出するブロック・ブートストラップを使う

また、母集団の分布の形を仮定できるなら、推定した分布から乱数でデータを作る パラメトリック・ブートストラップ も使えます。ここまでの方法は、これと区別してノンパラメトリック・ブートストラップとも呼ばれます。

2つの再標本化法を比べておきましょう。

	ジャックナイフ法	ブートストラップ法
作り方	1個ずつ除く	復元抽出
計算の回数	$n$ 回	$B$ 回（多い）
滑らかでない統計量	苦手（中央値など）	比較的使える
主な用途	偏りの修正、標準誤差	標準誤差、信頼区間

試験のポイント：ブートストラップ法は「復元抽出で大きさ $n$」です。非復元で $n$ 個を抽出すると、元の標本の並べ替えにしかなりません。

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

細かい式を忘れてしまっても問題ありません。ここでは、「分布は乱数で作れる」「積分は平均で近似できる」「標本から標本を作り直せば精度が分かる」という点だけ押さえておきましょう。

今回のまとめ

乱数で分布を作り、平均で解き、標本から標本を作り直す

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

今回は、計算機で統計の問題を解く方法について紹介しました。式で解けない問題も、乱数をたくさん作れば答えに近づけます。母集団に戻れないときでも、標本の中で再標本化すれば、推定量の揺れが見えてきます。

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

振り返り

• 疑似乱数は計算で作る。逆関数法は $F^{-1}(U)$ で一様乱数を目標の分布に変換する。
• 採択・棄却法は、屋根 $Mg$ の下でダーツを投げるイメージ。採択率は $1/M$ になる。
• モンテカルロ法は、積分を乱数の平均で近似する。誤差は $\sigma/\sqrt{N}$ で、高次元でも使える。
• ジャックナイフ法は1個ずつ抜いて、偏りと標準誤差を見積もる。
• ブートストラップ法は復元抽出で標本を作り直し、標本を母集団の代わりにする。

キーワードのおさらい

• 乱数: でたらめに並ぶ数の列。コンピュータでは計算で作る疑似乱数を使う
• 逆関数法: 一様乱数 $U$ を $F^{-1}(U)$ に変換して、分布 $F$ の乱数を作る方法
• 採択・棄却法: 作りやすい分布から候補を出し、確率 $f/(Mg)$ で採択して目標の分布の乱数を作る方法
• モンテカルロ法: 乱数を使って期待値や積分を近似する方法
• ジャックナイフ法: 1個ずつ除いた推定量から、偏りと標準誤差を見積もる方法
• ブートストラップ法: 標本からの復元抽出を繰り返して、推定量の分布を調べる方法

さて、準1級の旅も、次がいよいよ最後です。ここまでの手法は、ほとんどが「変数の間の関連」を調べるものでした。では、「原因と結果」は、データからどこまで言えるのでしょうか？

次回は、準1級の締めくくりとして「統計的因果推論」について紹介します。

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