モンテカルロ法の誤差と収束

推定はどれだけ外れうるか

モンテカルロ推定は n 個の乱数抽出の平均なので、それ自体がランダムです。典型的な誤差は標準誤差で、抽出値の標準偏差を sqrt(n) で割ったもの、t検定と同じ式です。2つの帰結が重要です。第一に、誤差は問題の見かけの難しさではなく、平均を取る量のばらつきだけで決まります。第二に、誤差はゆっくりしか縮みません。精度をもう1桁上げるには、抽出数が100倍必要です。

import math
import numpy as np

# Estimate E[exp(Z)] for a standard normal Z; the exact value is exp(0.5) = 1.6487
for n in (100, 10000, 1000000):
    values = np.exp(np.random.default_rng(5).standard_normal(n))
    print(n, round(values.mean(), 4), round(values.std(ddof=1) / math.sqrt(n), 5))
# 100 1.2167 0.14583
# 10000 1.6887 0.0217        the error falls roughly like 1 / sqrt(n)
# 1000000 1.6508 0.00216

誤差範囲と精度のコスト

シミュレーション値は信頼区間 平均 -/+ 1.96 * 標準誤差（真の値を約95%の確率で含む範囲）とともに報告します。許容誤差内の答えが必要なら、標準誤差の式を逆に使います：n = (std / 目標誤差)²（切り上げ）。ばらつきが約2.2、目標0.001なら約500万回、0.01なら4.7万回で済みます。まず小さな予備シミュレーションで std を知り、本番の抽出回数を決めましょう。

values = np.exp(np.random.default_rng(5).standard_normal(10000))
mean, std = values.mean(), values.std(ddof=1)
std_error = std / math.sqrt(len(values))

print(round(mean - 1.96 * std_error, 4), round(mean + 1.96 * std_error, 4))   # 1.6462 1.7312   contains the exact 1.6487

print(math.ceil((std / 0.01) ** 2))      # 47111      draws for a standard error of 0.01
print(math.ceil((std / 0.001) ** 2))     # 4710946    10 times the precision, 100 times the draws

安く精度を上げる：対称変量法

たくさん引く以外に、賢く引くこともできます。標準正規の Z に対して、鏡像の -Z も同じく有効な抽出です。f が単調なら f(Z) と f(-Z) は逆方向に動くので、その平均はどちら単独よりずっとノイズが小さくなります。対称変量推定量は抽出値ごとに (f(z) + f(-z)) / 2 の平均を取ります。例では、同じ数の乱数で標準誤差が0.0217から0.0119に下がり、3z + 1 のような奇関数ではノイズが完全に消えます。これは分散削減の技法の1つです。

z = np.random.default_rng(5).standard_normal(10000)
plain = np.exp(z)
antithetic = (np.exp(z) + np.exp(-z)) / 2           # each draw is paired with its mirror image

print(round(plain.mean(), 4), round(plain.std(ddof=1) / 100, 5))            # 1.6887 0.0217
print(round(antithetic.mean(), 4), round(antithetic.std(ddof=1) / 100, 5))  # 1.6565 0.01192   half the noise, same 10,000 numbers
print(round(math.exp(0.5), 4))                                              # 1.6487           the exact value

# The running mean of one long simulation settles down as n grows
values = np.exp(np.random.default_rng(5).standard_normal(100000))
running = np.cumsum(values) / np.arange(1, len(values) + 1)
print([round(float(running[k - 1]), 4) for k in (10, 100, 1000, 100000)])   # [1.5522, 1.2167, 1.634, 1.6502]

• 標準誤差 = std / sqrt(n)。抽出数を100倍にすると誤差は10分の1。mean -/+ 1.96 * se を報告する
• 目標誤差に必要な抽出数は ceil((std / 目標) ** 2)。まず小さな予備実行で std を知る
• 対称変量法：(f(z) + f(-z)) / 2 を平均すれば、乱数を増やさずにノイズを減らせる。別のシードで収束を必ず確認する

次は現在価値と割引です。債券への第一歩で、明日の1ユーロは今日の1ユーロより価値が低いことを学びます。

ラボ：採点される内容

• standard_error で推定値のノイズを測る25点
• mc_confidence_interval で答えに誤差範囲を付ける25点
• samples_needed でシミュレーションの規模を見積もる25点
• antithetic_mean で同じ精度をより少ない計算で得る25点