層化サンプリング｜シミュレーション

🎓 レベル：標準　|　重要度：A（必須）

📎 前提：制御変量法　|　関連：重点サンプリング

要点（BLUF）

層化サンプリング：定義域を**層（区間）**に分け、各層から決まった数だけサンプルします。全域でランダムに引くより配置が均等になります。
層への配分を設計で固定するので、層間のばらつきが消え、層内の変動だけが残る。これが分散減少の源です。
$\int_0^1 e^x dx$ を10層に均等配分すると、分散が約94倍縮小します。

1. アイデア：くじを均等に配る

プレーンなモンテカルロは $[0,1]$ 全域から一様に $n$ 点を引きますが、運次第で「左半分に偏る」ことが起きます。このサンプル配置のムラが余計な分散を生みます。

層化サンプリングは、定義域を $K$ 個の層 $S_1,\dots,S_K$ （ここでは等幅区間）に分け、各層から必ず $n_k$ 点ずつ引きます。配置のムラが設計で抑えられるので、推定が安定します。各層内で平均を取り、層の確率（幅）で重み付けて合算します。

\hat{I}_{\text{strat}} = \sum_{k=1}^{K} p_k\, \bar{g}_k,\qquad \bar{g}_k = \frac{1}{n_k}\sum_{i \in S_k} g(X_i),\quad p_k = P(X \in S_k)

2. なぜ分散が下がるか

全分散は「層間分散＋層内分散」に分解できます（分散分析と同じ）。

\text{Var}(g) = \underbrace{\sum_k p_k (\mu_k - \mu)^2}_{\text{層間}} + \underbrace{\sum_k p_k\,\sigma_k^2}_{\text{層内}}

層化サンプリングは各層の配分を固定するので、推定量の分散には層間の項が一切寄与せず、層内分散 $\sum_k p_k^2 \sigma_k^2 / n_k$ だけが残ります。関数が層をまたいで大きく動く（層間分散が大きい）ほど、消える量が大きく効果的。比例配分（ $n_k \propto p_k$ ）なら、分散は層内分散の加重平均まで縮みます。層を細かくすれば各層内で $g$ がほぼ一定になり、分散はさらに激減します。

3. 実測

$\int_0^1 e^x dx$ を等幅10層に分け、各層に $n/K$ 点ずつ均等配分します。

import numpy as np

# 乱数シードを固定
rng = np.random.default_rng(32)
n = 10_000
K = 10                                  # 層の数
reps = 2000

def plain(rng):
    u = rng.random(n)
    return np.exp(u).mean()

def stratified(rng):
    m = n // K                          # 各層のサンプル数
    total = 0.0
    for k in range(K):
        u = (k + rng.random(m)) / K     # 層 [k/K, (k+1)/K) 内の一様点
        total += np.exp(u).sum()
    return total / (m * K)              # 等幅なので単純平均でOK

p_ests = np.array([plain(rng) for _ in range(reps)])
s_ests = np.array([stratified(rng) for _ in range(reps)])

print(f"プレーン     平均={p_ests.mean():.5f}  分散={p_ests.var():.2e}")
print(f"層化         平均={s_ests.mean():.5f}  分散={s_ests.var():.2e}")
print(f"分散減少率   = {p_ests.var()/s_ests.var():.2f} 倍")

出力：

プレーン     平均=1.71817  分散=2.45e-05
層化         平均=1.71829  分散=2.62e-07
分散減少率   = 93.75 倍

出力の意味：同じ1万標本で、層化は分散を約94倍下げました。等幅10層に均等配分するだけで、これだけの改善。層内で $e^x$ がほぼ直線とみなせるほど各区間が狭いので、層内分散が小さく抑えられています。層を増やせばさらに下がりますが、各層に最低限のサンプルが要るので、層数とサンプル数のバランスを取ります。

4. 配分の最適化（ネイマン配分）

均等配分（ $n_k$ 一定）や比例配分（ $n_k \propto p_k$ ）が基本ですが、最適なのはネイマン配分 $n_k \propto p_k \sigma_k$ ——ばらつきの大きい層に多くサンプルを割く。層内分散 $\sigma_k$ が事前に分からなければ、予備実験で見積もるか、変動の大きそうな層を厚くします。層化は重点サンプリングと発想が近く、「重要な（変動の大きい）ところに資源を集中する」点で共通します。多次元ではラテン超方格法（各軸を層化）が実用的です。

数式の直観的意味

層化は「運任せの配置を、計画的な配置に置き換える」操作です。全域ランダムだと「たまたま左に偏る」リスクがあり、それが分散の一部（層間分散）。層ごとに枠を決めれば、その偏りは原理的に起きない——消えるのは「どの層から何点引くか」のばらつきで、残るのは「各層の中での $g$ の揺れ」だけ。だから層をまたぐ大きな構造を関数が持つほど効く。 $e^x$ は $[0,1]$ で約2.7倍変化する単調関数なので層間分散が大きく、それを消した効果が94倍として現れました。

⚠️ よくある誤解・落とし穴

「層の幅で重みを忘れる」：層が不等幅なら $p_k$ （幅）で重み付け必須。等幅だから単純平均で済んでいるだけ。
「層を増やせば無限に良くなる」ではない：層数を増やすと各層のサンプルが減り、層内平均の推定が不安定に。サンプル数とのバランスが要ります。
「効果は関数によらず一定」ではない：層間で大きく動く関数ほど効果大。層内でほぼ一定の関数（既にフラット）には効きません。
「高次元でそのまま使える」ではない：各軸を $K$ 分割すると層が $K^d$ に爆発。高次元はラテン超方格法や層化と他手法の併用で対処します。
「層化は不偏でない」ではない：各層を正しく重み付ければ不偏です。

対応シミュレーション参照

本文の10層・均等配分の層化（default_rng(32)、94倍減）。