蒙特卡洛的误差与收敛

估计可能错多少？

蒙特卡洛估计是 n 个随机抽样的平均值，因此它本身是随机的。它的典型误差是标准误：抽样值的标准差除以 sqrt(n)，与 t 检验中的公式相同。有两个重要推论。第一，误差与问题看起来有多难无关，只取决于你所求平均的那个量的离散程度。第二，它缩小得很慢：要多得到一位精度，需要多 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 时需要近五百万次抽样；目标 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 这样的奇函数，噪声会完全消失。这是几种方差缩减技巧之一。

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 倍，误差只减小到十分之一；报告 mean -/+ 1.96 * se。
• 达到目标误差所需的抽样数：ceil((std / 目标) ** 2)；先做小规模试运行了解 std。
• 对偶变量：对 (f(z) + f(-z)) / 2 取平均，不增加随机数就能降低噪声；务必换种子检验收敛。

下一节：现值与折现，进入固定收益的第一步：明天的一欧元不如今天的一欧元值钱。

实验：评分内容

• 用 standard_error 衡量估计的噪声25 分
• 用 mc_confidence_interval 为答案加上误差范围25 分
• 用 samples_needed 为模拟做预算25 分
• 用 antithetic_mean 以更小代价获得同样精度25 分