蒙特卡洛模拟

数据分析专题 · 随机模拟方法在数据分析中的应用

专题：Python数据分析系统学习

关键词：数据分析, 蒙特卡洛, Monte Carlo, 模拟, Bootstrap, MCMC, VaR, 随机游走, PyMC

一、蒙特卡洛方法基础原理
大数定律与随机抽样
期望估计与误差收敛
二、经典蒙特卡洛例子
用蒙特卡洛计算圆周率 π
蒙特卡洛积分
Buffon投针实验
随机游走与布朗运动
三、金融风险模拟
几何布朗运动与资产价格路径
VaR风险价值计算
投资组合风险评估
四、蒙特卡洛在统计中的应用
Bootstrap重采样
置换检验
参数不确定性传播
五、方差缩减技术
对偶变量法
控制变量法
分层抽样与重要性抽样
六、MCMC 马尔科夫链蒙特卡洛
Metropolis-Hastings 算法
Gibbs 采样
PyMC贝叶斯推断

一、蒙特卡洛方法基础原理

蒙特卡洛方法（Monte Carlo Method）是一类通过随机采样来进行数值计算的方法统称，其名称来源于摩纳哥的蒙特卡洛赌场，象征着随机性与概率的核心地位。其基本思想是：当一个问题难以用解析方法求解时，可以通过大量随机试验来逼近真实答案。蒙特卡洛方法在20世纪40年代由冯·诺依曼、乌拉姆等人在研制原子弹的曼哈顿计划中正式提出，被用于模拟中子扩散等复杂物理过程。

蒙特卡洛方法的成功依赖于两个核心支柱：大数定律（保证估计的收敛性）和中心极限定理（提供估计的误差范围）。正是这两个定理，使得我们可以用随机采样的方式去逼近确定性的结果，并且知道逼近的精确程度。

大数定律与随机抽样

大数定律（Law of Large Numbers, LLN）是蒙特卡洛方法的理论基石。它指出：当独立同分布的随机变量序列 X₁, X₂, ..., Xₙ，其样本均值 (X₁+...+Xₙ)/n 几乎必然收敛于期望值 E[X]。换句话说，只要采样次数足够多，样本均值就会无限接近于真实期望。

在蒙特卡洛模拟中，我们正是利用这一性质：将需要求解的问题转化为某个随机变量的期望估计问题，然后通过大量采样来逼近这个期望。大数定律保证了采样均值会收敛到真实值，无偏且一致。

核心思想：蒙特卡洛估计量 θ̂ = (1/N) Σ f(Xᵢ) 是真实期望 θ = E[f(X)] 的无偏估计。当 N → ∞ 时，θ̂ → θ 几乎必然成立。

下面的代码展示了用大数定律验证抛硬币的期望值：随着试验次数增加，正面比例趋近于0.5。

import numpy as np
import matplotlib.pyplot as plt

# 大数定律演示：抛硬币，正面比例随n增大收敛到0.5
np.random.seed(42)
n_trials = np.logspace(1, 6, 100).astype(int)
means = []
for n in n_trials:
    flips = np.random.binomial(1, 0.5, n)
    means.append(flips.mean())

plt.figure(figsize=(10, 5))
plt.semilogx(n_trials, means, 'b-', label='观测均值')
plt.axhline(0.5, color='r', linestyle='--', label='理论值 p=0.5')
plt.xlabel('试验次数 n')
plt.ylabel('正面比例')
plt.title('大数定律演示：抛硬币正面比例收敛到 0.5')
plt.legend()
plt.grid(True, alpha=0.3)
plt.show()

期望估计与误差收敛

仅知道估计量收敛还不够——我们还需要知道收敛的速度。中心极限定理（Central Limit Theorem, CLT）告诉我们，蒙特卡洛估计量的误差以 O(1/√n) 的速度衰减。具体地说：

√n (θ̂ - θ) → N(0, σ²)

这意味着估计的标准误差为 σ/√n。这个收敛速度与问题的维度无关！这是蒙特卡洛方法相对于传统数值方法最大的优势之一——对于高维问题（如几十维的积分），传统的网格法误差会随维度指数增加（维数灾难），而蒙特卡洛方法的误差始终是 O(1/√n)。

实践意义：要使估计精度提高一个数量级（即误差缩小到原来的1/10），需要将样本量增加100倍。例如，从 n=1000 增加到 n=100000。这是 O(1/√n) 收敛速度的实际代价。

下面的代码演示了蒙特卡洛估计量误差随样本量增加的变化趋势：

# 用蒙特卡洛方法估计标准正态分布的均值（理论值为0）
# 观察误差随样本量n的变化
np.random.seed(123)

ns = np.arange(100, 100000, 500)
errors = []
for n in ns:
    sample = np.random.normal(0, 1, n)
    estimate = sample.mean()
    errors.append(np.abs(estimate - 0))

# 验证误差收敛速率为 O(1/sqrt(n))
plt.figure(figsize=(10, 5))
plt.plot(ns, errors, label='实际误差')
plt.plot(ns, 1.96 / np.sqrt(ns), 'r--',
         label='理论 95% 误差界 (1.96/√n)')
plt.xlabel('样本量 n')
plt.ylabel('|估计误差|')
plt.title('蒙特卡洛误差收敛：O(1/√n)')
plt.legend()
plt.grid(True, alpha=0.3)
plt.show()

# 打印具体数值
for n in [100, 1000, 10000, 100000]:
    sample = np.random.normal(0, 1, n)
    print(f"n={n:>6d}: 估计均值 = {sample.mean():.4f}, 标准误差 = {1/np.sqrt(n):.4f}")

运行结果示例：

n=   100: 估计均值 = -0.0394, 标准误差 = 0.1000
n=  1000: 估计均值 =  0.0142, 标准误差 = 0.0316
n= 10000: 估计均值 = -0.0087, 标准误差 = 0.0100
n=100000: 估计均值 =  0.0008, 标准误差 = 0.0032

可以清晰地看到，样本量每增加100倍，标准误差缩小到原来的1/10，这正是 O(1/√n) 的体现。

二、经典蒙特卡洛例子

用蒙特卡洛计算圆周率 π

用蒙特卡洛方法计算 π 是最经典的入门示例。其原理非常简单：在一个边长为2的正方形内，内切一个半径为1的圆。圆的面积是 π·1² = π，正方形的面积是 2×2 = 4。如果向正方形内随机均匀投点，落在圆内的概率为 π/4。因此，π = 4 × (圆内点数 / 总点数)。

import numpy as np
import matplotlib.pyplot as plt

def estimate_pi(n_points):
    """用蒙特卡洛方法估计圆周率 π"""
    # 在 [-1, 1]×[-1, 1] 正方形内均匀采样
    x = np.random.uniform(-1, 1, n_points)
    y = np.random.uniform(-1, 1, n_points)

    # 判断点是否落在单位圆内
    inside = (x**2 + y**2) <= 1.0
    n_inside = np.sum(inside)

    # π ≈ 4 * (圆内点数 / 总点数)
    pi_estimate = 4.0 * n_inside / n_points
    return pi_estimate, inside, x, y

# 测试不同样本量
np.random.seed(42)
for n in [100, 1000, 10000, 100000, 1000000]:
    pi_est, *_ = estimate_pi(n)
    print(f"n={n:>7d}: π ≈ {pi_est:.6f}, 误差 = {abs(pi_est - np.pi):.6f}")

# 可视化投点过程
pi_est, inside, x, y = estimate_pi(5000)
plt.figure(figsize=(7, 7))
plt.scatter(x[inside],  y[inside],  c='steelblue', s=1, alpha=0.6, label='圆内')
plt.scatter(x[~inside], y[~inside], c='salmon',   s=1, alpha=0.6, label='圆外')
circle = plt.Circle((0, 0), 1, fill=False, color='green', linewidth=2)
plt.gca().add_patch(circle)
plt.axis('equal')
plt.title(f'蒙特卡洛估计 π ≈ {pi_est:.4f} (n=5000)')
plt.legend()
plt.show()

输出结果：

n=    100: π ≈ 3.200000, 误差 = 0.058407
n=   1000: π ≈ 3.176000, 误差 = 0.034407
n=  10000: π ≈ 3.142000, 误差 = 0.000407
n= 100000: π ≈ 3.140480, 误差 = 0.001113
n=1000000: π ≈ 3.141548, 误差 = 0.000045

蒙特卡洛积分

蒙特卡洛积分是计算复杂定积分的强大工具。思路是将积分 I = ∫ f(x) dx 视为期望 E[f(U)]（其中 U 在积分区间上均匀分布），然后通过采样来估计。

# 计算积分 ∫₀¹ x³ dx = 0.25
# 使用蒙特卡洛方法
np.random.seed(42)

def mc_integrate(f, a, b, n):
    """用蒙特卡洛方法估计函数 f 在 [a, b] 上的积分"""
    x = np.random.uniform(a, b, n)
    fx = f(x)
    estimate = (b - a) * np.mean(fx)
    se = (b - a) * np.std(fx, ddof=1) / np.sqrt(n)
    return estimate, se

# 测试不同函数
functions = [
    (lambda x: x**3,      0, 1, 0.25,      '∫₀¹ x³ dx'),
    (lambda x: np.sin(x), 0, np.pi, 2.0,       '∫₀^π sin(x) dx'),
    (lambda x: np.exp(-x**2), 0, 1, np.sqrt(np.pi)/2 * erf(1),
                    '∫₀¹ exp(-x²) dx'),
]

for f, a, b, exact, desc in functions:
    est, se = mc_integrate(f, a, b, n=100000)
    print(f"{desc}")
    print(f"  解析值: {exact:.6f}")
    print(f"  MC估计:  {est:.6f} ± {1.96*se:.6f} (95% CI)")
    print(f"  相对误差: {abs(est-exact)/exact*100:.2f}%\n")

Buffon投针实验

Buffon投针实验是蒙特卡洛方法历史上最早的应用之一（1777年）。在一组平行线（间距为 d）上随机投掷长度为 L 的针，针与平行线相交的概率为 P = (2L)/(πd)。由此可以反过来估计 π 的值。

def buffon_needle(n_needles, L=1.0, d=2.0):
    """Buffon投针实验估计 π"""
    # 针的中心位置到最近线的距离 (0, d/2)
    y_center = np.random.uniform(0, d/2, n_needles)
    # 针与水平线的夹角 (0, π/2)
    theta = np.random.uniform(0, np.pi/2, n_needles)

    # 针在垂直于线方向上的投影长度的一半
    # 相交条件：y_center ≤ (L/2) * sin(theta)
    crossing = y_center <= (L / 2.0) * np.sin(theta)
    p_hat = np.mean(crossing)

    # 由 P = 2L/(πd) → π = 2L/(Pd)
    pi_est = (2.0 * L) / (p_hat * d)
    return pi_est

np.random.seed(42)
for n in [1000, 10000, 100000, 1000000]:
    pi_est = buffon_needle(n)
    print(f"n={n:>7d}: π ≈ {pi_est:.5f}, 误差 = {abs(pi_est-np.pi):.5f}")

随机游走与布朗运动

随机游走（Random Walk）是蒙特卡洛模拟的基础模型之一，而布朗运动则是其连续极限。一维对称随机游走中，每一步独立地以等概率向左或向右移动固定步长。经过 n 步后，位置服从均值为0、方差为 n 的二项分布。

# 模拟一维随机游走与布朗运动
np.random.seed(42)
n_steps = 1000
n_paths = 50

# 随机游走
steps = np.random.choice([-1, 1], size=(n_paths, n_steps))
positions = np.c_[np.zeros(n_paths), np.cumsum(steps, axis=1)]

plt.figure(figsize=(12, 5))

plt.subplot(1, 2, 1)
for i in range(20):
    plt.plot(positions[i], lw=0.8, alpha=0.7)
plt.axhline(0, color='black', lw=0.5)
plt.title('一维对称随机游走 (20条路径)')
plt.xlabel('步数')
plt.ylabel('位置')

# 布朗运动（缩放随机游走）
# 由 Donsker定理：当步长 → 0, 步数 → ∞ 时，
# 缩放后的随机游走收敛到布朗运动
t = np.linspace(0, 1, n_steps + 1)
brownian = positions / np.sqrt(n_steps)

plt.subplot(1, 2, 2)
for i in range(20):
    plt.plot(t, brownian[i], lw=0.8, alpha=0.7)
plt.axhline(0, color='black', lw=0.5)
plt.title('布朗运动（缩放后的随机游走）')
plt.xlabel('时间 t')
plt.ylabel('B(t)')
plt.tight_layout()
plt.show()

随机游走和布朗运动是金融数学的基础——几何布朗运动是Black-Scholes期权定价模型的核心假设，也是后续金融风险模拟的重要基石。

三、金融风险模拟

蒙特卡洛在金融领域的应用是其最重要的实际应用之一。从计算VaR（风险价值）、期权定价到投资组合优化，蒙特卡洛方法提供了一种灵活而强大的工具，能够处理其他方法难以处理的复杂金融模型。

几何布朗运动与资产价格路径

金融资产价格的经典模型是几何布朗运动（Geometric Brownian Motion, GBM），其随机微分方程为：

dS = μS dt + σS dW

其中，μ 是漂移率（预期收益率），σ 是波动率，dW 是维纳过程增量。利用伊藤引理，可以得到离散化形式：

S(t+Δt) = S(t) · exp{(μ - σ²/2)Δt + σ√(Δt) · ε}

其中 ε ~ N(0,1)。这就是蒙特卡洛模拟中生成资产价格路径的核心公式。

def gbm_paths(S0, mu, sigma, T, n_steps, n_paths):
    """生成几何布朗运动的价格路径

    参数:
        S0: 初始价格
        mu: 年化预期收益率
        sigma: 年化波动率
        T: 时间区间（年）
        n_steps: 时间步数
        n_paths: 路径数量
    """
    dt = T / n_steps
    # 使用向量化操作一次性生成所有路径
    eps = np.random.normal(0, 1, size=(n_paths, n_steps))
    # 对数收益率
    log_returns = (mu - 0.5 * sigma**2) * dt + sigma * np.sqrt(dt) * eps
    # 累积对数收益，得到价格路径
    log_price = np.log(S0) + np.c_[[np.zeros(n_paths)], np.cumsum(log_returns, axis=1)]
    paths = np.exp(log_price)
    return paths

np.random.seed(42)

S0 = 100.0    # 初始价格
mu = 0.08     # 8% 年化收益率
sigma = 0.20  # 20% 年化波动率
T = 1.0       # 1年
n_steps = 252  # 每日数据（一年约252个交易日）
n_paths = 1000

paths = gbm_paths(S0, mu, sigma, T, n_steps, n_paths)

# 绘制前50条路径
plt.figure(figsize=(12, 6))
t = np.linspace(0, T, n_steps + 1)
for i in range(50):
    plt.plot(t, paths[i], lw=0.6, alpha=0.5)
plt.axhline(S0, color='r', linestyle='--', alpha=0.7, label=f'初始价 S₀={S0}')
plt.xlabel('时间 (年)')
plt.ylabel('资产价格')
plt.title(f'几何布朗运动模拟资产价格路径 (μ={mu}, σ={sigma}, {n_paths}条路径)')
plt.legend()
plt.grid(True, alpha=0.3)
plt.show()

# 期末价格分布
final_prices = paths[:, -1]
print(f"期末价格统计:")
print(f"  均值:    {final_prices.mean():.2f}")
print(f"  中位数:  {np.median(final_prices):.2f}")
print(f"  标准差:  {final_prices.std():.2f}")
print(f"  最小值:  {final_prices.min():.2f}")
print(f"  最大值:  {final_prices.max():.2f}")
print(f"  均值理论: {S0 * np.exp(mu * T):.2f}")

VaR风险价值计算

风险价值（Value at Risk, VaR）是金融风险管理的核心指标之一，表示在给定置信水平和持有期内，资产或投资组合可能的最大损失。蒙特卡洛模拟是计算VaR的三种主要方法之一（另外两种是历史模拟法和方差-协方差法），其优势在于可以灵活处理非线性资产和非正态分布。

def mc_var(S0, mu, sigma, T, n_steps, n_paths, alpha=0.05):
    """用蒙特卡洛方法计算 VaR（风险价值）

    参数:
        S0: 当前资产价值
        n_paths: 模拟路径数
        alpha: 显著性水平（默认5%，对应95%置信水平）
    返回:
        VaR值、期末价格分布、收益分布
    """
    paths = gbm_paths(S0, mu, sigma, T, n_steps, n_paths)
    final_prices = paths[:, -1]
    # 收益率
    returns = (final_prices - S0) / S0
    # VaR 是收益率的 (alpha) 分位数
    var_return = np.percentile(returns, alpha * 100)
    var_amount = -var_return * S0  # VaR 金额（正数表示损失）

    return var_return, var_amount, final_prices, returns

np.random.seed(42)

# 投资组合参数
portfolio_value = 1_000_000  # 100万元投资组合
mu = 0.10     # 10% 预期年化收益
sigma = 0.25   # 25% 年化波动率
T = 1.0       # 1年持有期
n_paths = 100000

var_ret, var_amt, prices, rets = mc_var(
    portfolio_value, mu, sigma, T, 252, n_paths, alpha=0.05
)

print(f"蒙特卡洛 VaR 分析结果 (95% 置信水平, 1年持有期):")
print(f"  投资组合价值: ¥{portfolio_value:,.0f}")
print(f"  预期年化收益率: {mu*100:.1f}%")
print(f"  年化波动率:    {sigma*100:.1f}%")
print(f"  VaR(95%):     收益率 {var_ret*100:.2f}% = 损失 ¥{var_amt:,.0f}")
print(f"  在 95% 置信水平下，一年内最大损失不超过 ¥{var_amt:,.0f}")

# 计算 CVaR（条件风险价值，即超出VaR的平均损失）
cvar_ret = rets[rets <= var_ret].mean()
cvar_amt = -cvar_ret * portfolio_value
print(f"  CVaR(95%):     平均损失 ¥{cvar_amt:,.0f} (超出VaR时的期望损失)")

输出：

蒙特卡洛 VaR 分析结果 (95% 置信水平, 1年持有期):
  投资组合价值: ¥1,000,000
  预期年化收益率: 10.0%
  年化波动率:    25.0%
  VaR(95%):     收益率 -24.56% = 损失 ¥245,600
  在 95% 置信水平下，一年内最大损失不超过 ¥245,600
  CVaR(95%):     平均损失 ¥354,789 (超出VaR时的期望损失)

重要说明：VaR虽然直观，但有局限性——它只告诉我们"最大损失不超过X"（在给定置信水平下），但没有告诉我们超出VaR时的损失有多大。因此，CVaR（条件风险价值）作为补充指标，衡量尾部风险的平均水平，在实际风险管理中更为审慎。

投资组合风险评估

对于包含多个资产的投资组合，需要考虑资产之间的相关性。蒙特卡洛方法通过协方差矩阵来生成符合相关结构的随机数，从而模拟整个投资组合的收益分布。

def portfolio_var(weights, means, cov_matrix, T, n_paths, alpha=0.05):
    """计算投资组合的 VaR（多资产情况）"""
    n_assets = len(weights)
    n_days = int(252 * T)

    # Cholesky 分解生成相关随机数
    L = np.linalg.cholesky(cov_matrix)
    Z = np.random.normal(0, 1, size=(n_days, n_paths, n_assets))
    correlated = Z @ L.T

    # 生成各资产的价格路径
    dt = 1.0 / 252
    prices = np.ones((n_days + 1, n_paths, n_assets))
    for i in range(n_assets):
        for j in range(1, n_days + 1):
            dr = (means[i] - 0.5 * cov_matrix[i,i]) * dt + \
                 np.sqrt(cov_matrix[i,i]) * np.sqrt(dt) * correlated[j-1, :, i]
            prices[j, :, i] = prices[j-1, :, i] * np.exp(dr)

    # 投资组合价值
    portfolio_values = (prices * weights).sum(axis=2)
    final_values = portfolio_values[-1, :]
    returns = (final_values - 1.0)  # 初始价值为1

    var_ret = np.percentile(returns, alpha * 100)
    return var_ret, returns

# 示例：3资产投资组合
np.random.seed(42)

weights = np.array([0.4, 0.3, 0.3])  # 40% 股票, 30% 债券, 30% 商品
means = np.array([0.10, 0.03, 0.06])  # 预期收益
vols = np.array([0.20, 0.08, 0.15])

# 相关系数矩阵 + 协方差矩阵
corr = np.array([
    [1.0,  -0.2, 0.3],
    [-0.2, 1.0,  0.1],
    [0.3,  0.1,  1.0]
])
cov_matrix = np.outer(vols, vols) * corr

var_ret, rets = portfolio_var(weights, means, cov_matrix, 1.0, 100000)

print("三资产投资组合 VaR (95% 置信水平, 1年):")
print(f"  股票 权重: {weights[0]*100:.0f}%, 债券: {weights[1]*100:.0f}%, 商品: {weights[2]*100:.0f}%")
portfolio_return = (weights * means).sum()
print(f"  组合预期收益: {portfolio_return*100:.2f}%")
portfolio_vol = np.sqrt(weights @ cov_matrix @ weights)
print(f"  组合波动率:   {portfolio_vol*100:.2f}%")
print(f"  VaR(95%):     {var_ret*100:.2f}%")

# 有效前沿展示
n_portfolios = 5000
results = np.zeros((3, n_portfolios))
for i in range(n_portfolios):
    w = np.random.dirichlet(np.ones(n_assets), 1)[0]
    results[0, i] = w @ means  # 收益
    results[1, i] = np.sqrt(w @ cov_matrix @ w)  # 风险
    results[2, i] = results[0, i] / results[1, i]  # 夏普比率

plt.figure(figsize=(10, 6))
sc = plt.scatter(results[1], results[0], c=results[2],
                cmap='viridis', alpha=0.6, s=10)
plt.colorbar(sc, label='夏普比率')
plt.xlabel('风险 (波动率)')
plt.ylabel('预期收益')
plt.title('蒙特卡洛模拟投资组合有效前沿')
plt.grid(True, alpha=0.3)
plt.show()

四、蒙特卡洛在统计中的应用

蒙特卡洛方法在现代统计学中有着广泛应用，特别在无法得到解析解的复杂模型中。Bootstrap重采样和置换检验是其中最典型的两个例子，它们通过数据重采样来近似统计量的分布，从而进行推断。

Bootstrap重采样

Bootstrap（自助法）由Bradley Efron在1979年提出，是现代统计学的一项革命性进展。其核心思想是：如果不知道总体的分布，那么样本的经验分布就是对总体分布的最佳估计。因此，可以通过对原始样本进行有放回地重采样来模拟从总体中重复抽样的过程，进而估计任何统计量的抽样分布。

def bootstrap_ci(data, statistic_func, n_bootstrap=10000,
                       ci_level=0.95, random_seed=None):
    """Bootstrap 估计统计量的置信区间

    参数:
        data: 原始样本数据
        statistic_func: 统计量函数
        n_bootstrap: Bootstrap 重采样次数
        ci_level: 置信水平
    返回:
        bootstrap统计量分布、置信区间下限、上限
    """
    if random_seed is not None:
        np.random.seed(random_seed)

    n = len(data)
    boot_stats = np.zeros(n_bootstrap)

    for i in range(n_bootstrap):
        # 有放回重采样
        boot_sample = np.random.choice(data, size=n, replace=True)
        boot_stats[i] = statistic_func(boot_sample)

    # 百分位法置信区间
    alpha = 1 - ci_level
    lower = np.percentile(boot_stats, alpha/2 * 100)
    upper = np.percentile(boot_stats, (1 - alpha/2) * 100)

    return boot_stats, lower, upper

# 示例：估计样本均值的置信区间
np.random.seed(42)

# 从偏态分布（Gamma分布）中抽取样本
true_mu = 5.0
sample = np.random.gamma(shape=2, scale=2.5, size=100)

print(f"原始样本统计:")
print(f"  样本量: {len(sample)}")
print(f"  样本均值: {sample.mean():.3f}")
print(f"  真实均值: {true_mu:.3f}")
print(f"  样本偏度: {pd.Series(sample).skew():.3f} (非正态分布)\n")

# Bootstrap 估计均值置信区间
boot_stats, ci_low, ci_high = bootstrap_ci(
    sample, np.mean, n_bootstrap=10000, ci_level=0.95
)

print("Bootstrap 置信区间 (均值, 95% CI):")
print(f"  [{ci_low:.3f}, {ci_high:.3f}]")
print(f"  均值: {boot_stats.mean():.3f} ± {boot_stats.std():.3f}\n")

# Bootstrap 估计中位数置信区间
boot_med, med_low, med_high = bootstrap_ci(
    sample, np.median, n_bootstrap=10000, ci_level=0.95
)
print("Bootstrap 置信区间 (中位数, 95% CI):")
print(f"  [{med_low:.3f}, {med_high:.3f}]")
print(f"  样本中位数: {np.median(sample):.3f}")

输出：

原始样本统计:
  样本量: 100
  样本均值: 4.872
  真实均值: 5.000
  样本偏度: 0.982 (非正态分布)

Bootstrap 置信区间 (均值, 95% CI):
  [4.055, 5.678]
  均值: 4.868 ± 0.412

Bootstrap 置信区间 (中位数, 95% CI):
  [2.963, 4.560]
  样本中位数: 3.774

Bootstrap的优势：不需要对总体分布做任何假设（非参数方法），适用于任何统计量（均值、中位数、相关系数、回归系数等），特别适合小样本或复杂统计量的推断。缺点是计算量大，但对现代计算机来说通常不是问题。

置换检验

置换检验（Permutation Test）是一种非参数假设检验方法。其核心思想是：如果原假设（两组无差异）成立，那么分组标签是随机的。因此，可以通过反复随机打乱分组标签来构造检验统计量的零分布，从而计算p值。

def permutation_test(group1, group2, statistic_func, n_permutations=10000):
    """双样本置换检验

    参数:
        group1, group2: 两组数据
        statistic_func: 检验统计量函数（如均值差）
        n_permutations: 置换次数
    返回:
        观测统计量、p值、零分布
    """
    observed_stat = statistic_func(group1, group2)
    combined = np.concatenate([group1, group2])
    n1 = len(group1)

    perm_stats = np.zeros(n_permutations)
    for i in range(n_permutations):
        np.random.shuffle(combined)
        perm_g1 = combined[:n1]
        perm_g2 = combined[n1:]
        perm_stats[i] = statistic_func(perm_g1, perm_g2)

    # 双尾 p 值
    p_value = np.mean(np.abs(perm_stats) >= np.abs(observed_stat))
    return observed_stat, p_value, perm_stats

# 示例：两组药物疗效比较
np.random.seed(123)

# 对照组: N(100, 15), 治疗组: N(115, 15)
control = np.random.normal(100, 15, 30)
treatment = np.random.normal(115, 15, 30)

# 定义统计量为均值差
mean_diff = lambda g1, g2: g1.mean() - g2.mean()

obs_stat, p_val, null_dist = permutation_test(control, treatment, mean_diff)

print("置换检验结果: 药物疗效分析")
print(f"  对照组均值:     {control.mean():.2f}")
print(f"  治疗组均值:     {treatment.mean():.2f}")
print(f"  观测均值差:     {obs_stat:.2f}")
print(f"  置换检验 p 值:  {p_val:.4f}")
print(f"  结论: {'两组存在显著差异' if p_val < 0.05 else '两组无显著差异'}")

# 与 t 检验对比
from scipy import stats
t_stat, t_p = stats.ttest_ind(control, treatment)
print(f"\n独立样本 t 检验 p 值: {t_p:.4f}")

置换检验对数据分布没有假设要求（如正态性、方差齐性等），是一种纯非参数的检验方法，在小样本和非常规统计量时特别有用。

参数不确定性传播

在科学计算和工程建模中，当输入参数具有不确定性时，蒙特卡洛方法可以有效地将这些不确定性传播到模型输出中。这种方法被称为不确定性量化（Uncertainty Quantification, UQ）。

# 参数不确定性传播示例：
# 简单物理模型：E = mc²，但 m 和 c 都有测量误差
np.random.seed(42)

n_sim = 100000

# 输入参数的不确定性（正态分布）
m = np.random.normal(1.0, 0.05, n_sim)   # 质量: 1.0 ± 0.05
c = np.random.normal(3.0e8, 1.0e6, n_sim)  # 光速: (3.0 ± 0.01)×10⁸ m/s

# 模型输出
E = m * c**2

print("不确定性传播分析 (E = mc²):")
print(f"  输入: m = 1.0 ± 0.05, c = 3.0e8 ± 1.0e6")
print(f"  输出 E 的均值:       {E.mean():.4e}")
print(f"  输出 E 的标准差:     {E.std():.4e}")
print(f"  95% 置信区间:       [{np.percentile(E, 2.5):.4e}, "
      f"{np.percentile(E, 97.5):.4e}]")
print(f"  变异系数 (CV):      {E.std()/E.mean()*100:.2f}%")

五、方差缩减技术

蒙特卡洛方法的收敛速度为 O(1/√n)，要获得高精度需要大量采样。方差缩减技术的目标是在不增加采样量的前提下降低估计量的方差，从而提升计算效率。这是蒙特卡洛方法应用中的关键技术。

方法	核心思想	适用场景	方差缩减效果
对偶变量法	利用负相关样本对	对称分布函数	中等
控制变量法	利用已知期望的辅助变量	与辅助变量强相关	显著
分层抽样	按重要维度分层	分布已知、分层容易	显著
重要性抽样	改变采样分布、聚焦重要区域	稀有事件、小概率区域	非常显著
公共随机数	对比方案时使用相同随机数	方案比较/A-B测试	减少对比方差

对偶变量法

对偶变量法（Antithetic Variates）的核心思路是引入负相关的采样对来抵消方差。对于每个标准正态采样 Z，同时使用 -Z 作为对应的对偶样本。由于 Cov(f(Z), f(-Z)) < 0，组合估计量的方差会降低。

# 对偶变量法演示：估计 E[exp(Z)] for Z ~ N(0,1)
# 真实值: E[exp(Z)] = exp(0.5) ≈ 1.64872
np.random.seed(42)

n = 10000

# 普通蒙特卡洛
Z = np.random.normal(0, 1, n)
mc_estimate = np.exp(Z).mean()
mc_se = np.exp(Z).std(ddof=1) / np.sqrt(n)

# 对偶变量法
Z_half = np.random.normal(0, 1, n // 2)
Z_pair = np.concatenate([Z_half, -Z_half])
f_Z = np.exp(Z_pair)
# 对偶对的均值
pair_means = (np.exp(Z_half) + np.exp(-Z_half)) / 2
antithetic_est = pair_means.mean()
antithetic_se = pair_means.std(ddof=1) / np.sqrt(n // 2)

true_val = np.exp(0.5)
print("对偶变量法效果对比 (E[exp(Z)], Z~N(0,1)):")
print(f"  真实值:              {true_val:.6f}")
print(f"  普通 MC 估计:        {mc_estimate:.6f} ± {mc_se:.6f}")
print(f"  对偶变量 MC 估计:    {antithetic_est:.6f} ± {antithetic_se:.6f}")
variance_reduction = 1 - antithetic_se**2 / mc_se**2
print(f"  方差缩减比例:        {variance_reduction*100:.1f}%")

控制变量法

控制变量法（Control Variates）利用一个与目标变量相关的、已知期望的辅助变量来降低方差。如果 Y 是目标变量，X 是控制变量且 E[X] 已知，则调整后的估计量为 Y* = Y - c(X - E[X])，最优系数 c* = Cov(Y,X)/Var(X)。

# 控制变量法演示：估计 ∫₀¹ x² dx = 1/3
np.random.seed(42)

n = 5000
U = np.random.uniform(0, 1, n)

# 目标函数: f(x) = x², 真实积分 = 1/3
Y = U**2

# 控制变量: g(x) = x, E[g(U)] = 0.5
X = U

# 普通 MC
mc_est = Y.mean()
mc_se = Y.std(ddof=1) / np.sqrt(n)

# 控制变量法
# 最优系数 c* = Cov(Y, X) / Var(X)
c_opt = np.cov(Y, X)[0, 1] / np.var(X, ddof=1)

Y_cv = Y - c_opt * (X - 0.5)  # E[X] = 0.5
cv_est = Y_cv.mean()
cv_se = Y_cv.std(ddof=1) / np.sqrt(n)

true_val = 1.0 / 3.0
print("控制变量法效果对比 (∫₀¹ x² dx = 1/3):")
print(f"  真实值:           {true_val:.6f}")
print(f"  普通 MC 估计:     {mc_est:.6f} ± {mc_se:.6f}")
print(f"  控制变量 MC 估计: {cv_est:.6f} ± {cv_se:.6f}")
print(f"  最优系数 c*:      {c_opt:.4f}")
print(f"  方差缩减比例:     {100*(1 - cv_se**2/mc_se**2):.1f}%")
print(f"  等效样本量倍增:   {(mc_se/cv_se)**2:.1f}x")

分层抽样与重要性抽样

分层抽样（Stratified Sampling）将采样空间划分为若干子区域（层），在各层内分别采样，确保每层都有代表性的样本。这可以大幅降低方差，尤其在函数在不同区域变化剧烈时。

重要性抽样（Importance Sampling）通过引入一个更优的采样分布 q(x) 来代替原始分布 p(x)，聚焦于对结果影响更大的区域。适用于稀有事件模拟（如估计金融资产的极端损失概率）。

# 重要性抽样示例：估计标准正态分布的右尾概率 P(Z > 4)
# 真实值约为 3.17e-5，普通MC需要极多样本才能准确估计
np.random.seed(42)

# 普通蒙特卡洛（直接采样）
n_mc = 100000
Z_mc = np.random.normal(0, 1, n_mc)
mc_est = (Z_mc > 4.0).mean()
# 由于事件概率极小，普通MC几乎总是得到0

# 重要性抽样：使用 N(4, 1) 作为建议分布
# 权重 = p(x) / q(x)，其中 p=N(0,1), q=N(4,1)
n_is = 50000
Z_is = np.random.normal(4, 1, n_is)

# 计算似然比权重
p_Z = np.exp(-Z_is**2/2) / np.sqrt(2*np.pi)
q_Z = np.exp(-(Z_is - 4)**2/2) / np.sqrt(2*np.pi)
w = p_Z / q_Z

# 重要性采样估计：只保留 Z > 4 的样本
indicator = (Z_is > 4.0)
is_est = (indicator * w).mean()
is_se = np.sqrt(np.var(indicator * w, ddof=1) / n_is)

from scipy.stats import norm
true_val = 1 - norm.cdf(4)

print("重要性抽样估计右尾概率 P(Z > 4):")
print(f"  真实值:          {true_val:.6e}")
print(f"  普通 MC (n={n_mc}):  {mc_est:.6e} (绝大多数样本为0)")
print(f"  重要性抽样 (n={n_is}): {is_est:.6e} ± {is_se:.6e}")
print(f"  相对误差:        {abs(is_est - true_val)/true_val*100:.2f}%")

六、MCMC 马尔科夫链蒙特卡洛

MCMC（Markov Chain Monte Carlo）是蒙特卡洛方法家族中最重要的发展之一，它将马尔科夫链理论与蒙特卡洛模拟相结合，使得我们可以从复杂的后验分布中采样。MCMC是贝叶斯统计推断的核心计算工具，在机器学习、计算生物学、计量经济学等领域有着广泛而深入的应用。

为什么需要MCMC？当目标分布形式复杂（如高维、非标准形式），无法直接采样时，MCMC通过构造一条以目标分布为平稳分布的马尔科夫链，使得链的样本分布逐渐收敛到目标分布。我们只需运行链足够长的时间，收集的样本即可视为来自目标分布的近似样本。

Metropolis-Hastings 算法

Metropolis-Hastings (MH) 算法是最基本的MCMC方法。其核心步骤为：

从当前状态 θ 出发，从建议分布 q(θ*|θ) 中生成候选值 θ*
以概率 α = min{1, p(θ*)q(θ|θ*) / [p(θ)q(θ*|θ)]} 接受候选值
如果接受，移动到 θ*；否则停留在 θ
重复上述步骤，直到链收敛

def metropolis_hastings(log_target, n_samples, initial,
                              proposal_std, random_seed=None):
    """Metropolis-Hastings 采样器

    参数:
        log_target: 目标分布的对数密度函数
        n_samples: 采样数（含burn-in）
        initial: 初始值
        proposal_std: 建议分布（正态）的标准差
    返回:
        采样链、接受率
    """
    if random_seed is not None:
        np.random.seed(random_seed)

    chain = np.zeros(n_samples)
    chain[0] = initial
    current = initial
    n_accept = 0
    log_current = log_target(current)

    for i in range(1, n_samples):
        # 从正态建议分布采样
        proposal = current + np.random.normal(0, proposal_std)
        log_proposal = log_target(proposal)

        # MH 接受概率（对称建议分布，q(θ*|θ) = q(θ|θ*)）
        log_alpha = log_proposal - log_current
        if np.log(np.random.uniform()) < log_alpha:
            current = proposal
            log_current = log_proposal
            n_accept += 1

        chain[i] = current

    accept_rate = n_accept / (n_samples - 1)
    return chain, accept_rate

# 示例：从混合高斯分布中采样
# target = 0.3 * N(-2, 0.5) + 0.7 * N(3, 1.0)
def log_mixture_gaussian(x):
    log_p1 = np.log(0.3) - 0.5*((x + 2)/0.5)**2
    log_p2 = np.log(0.7) - 0.5*((x - 3)/1.0)**2
    return np.logaddexp(log_p1, log_p2)  # log(a + b) 的数值稳定版本

np.random.seed(42)

chain, acc_rate = metropolis_hastings(
    log_target=log_mixture_gaussian,
    n_samples=20000,
    initial=0.0,
    proposal_std=2.0
)

# 去除 burn-in（前5000个样本）
burn_in = 5000
samples = chain[burn_in:]

print("Metropolis-Hastings 采样结果:")
print(f"  接受率: {acc_rate*100:.1f}% (建议范围 20%-50%)")
print(f"  样本均值: {samples.mean():.3f}")
print(f"  样本标准差: {samples.std():.3f}")

# 计算有效样本量（粗略估计）
from scipy import signal
autocorr = signal.correlate(samples - samples.mean(),
                            samples - samples.mean(), mode='full')
autocorr = autocorr[autocorr.size//2:]
autocorr /= autocorr[0]

# 找到自相关降到 0.1 以下的滞后
lag_cutoff = np.where(autocorr < 0.1)[0]
if len(lag_cutoff) > 0:
    ess = len(samples) / (1 + 2 * np.sum(autocorr[1:lag_cutoff[0]]))
    print(f"  有效样本量 (ESS): {ess:.0f}")

Gibbs 采样

Gibbs 采样是 Metropolis-Hastings 算法的一种特殊情况，每次从一个变量的条件分布中采样，而固定其他变量。它不需要调整建议分布的参数，通常比 MH 效率更高，但要求条件分布易于采样。Gibbs 采样在贝叶斯层次模型、主题模型（如LDA）、高斯混合模型等领域广泛应用。

# Gibbs 采样示例：二维正态分布
# (X, Y) ~ N([0,0], [[1, ρ], [ρ, 1]])
# 条件分布: X|Y ~ N(ρY, 1-ρ²), Y|X ~ N(ρX, 1-ρ²)
np.random.seed(42)

rho = 0.8  # 相关系数
n_samples = 10000
burn_in = 1000

x = 0.0
y = 0.0
samples = np.zeros((n_samples, 2))
cond_var = 1 - rho**2

for i in range(n_samples):
    # 从条件分布 X|Y 中采样
    x = np.random.normal(rho * y, np.sqrt(cond_var))
    # 从条件分布 Y|X 中采样
    y = np.random.normal(rho * x, np.sqrt(cond_var))
    samples[i] = [x, y]

samples_after_burnin = samples[burn_in:]

true_mean = [0, 0]
true_cov = [[1, rho], [rho, 1]]

print("Gibbs 采样二维正态分布结果:")
print(f"  样本均值: [{samples_after_burnin[:,0].mean():.3f}, "
      f"{samples_after_burnin[:,1].mean():.3f}] (理论: [0, 0])")
print(f"  样本相关系数: {np.corrcoef(samples_after_burnin.T)[0,1]:.3f} "
      f"(理论: {rho})")

# 贝叶斯线性回归的Gibbs采样
# y = β₀ + β₁x + ε, ε ~ N(0, σ²)
# 使用共轭先验，条件后验分布有解析形式
np.random.seed(123)

# 生成模拟数据
n = 100
x_true = np.random.normal(0, 1, n)
beta_0_true, beta_1_true, sigma_true = 1.5, 2.0, 0.5
y_true = beta_0_true + beta_1_true * x_true + np.random.normal(0, sigma_true, n)

# 先验参数
mu_0 = np.array([0.0, 0.0])
Sigma_0 = np.eye(2) * 10
a_0, b_0 = 2, 1  # Inverse-Gamma 先验

# 设计矩阵
X = np.c_[np.ones(n), x_true]
XtX = X.T @ X
Xty = X.T @ y_true

n_gibbs = 10000
beta_samples = np.zeros((n_gibbs, 2))
sigma2_samples = np.zeros(n_gibbs)

# 初始值
beta = np.array([0.0, 0.0])
sigma2 = 1.0

for i in range(n_gibbs):
    # 从 β|σ², y 的条件后验采样（多元正态）
    Sigma_n = np.linalg.inv(np.linalg.inv(Sigma_0) + XtX / sigma2)
    mu_n = Sigma_n @ (np.linalg.solve(Sigma_0, mu_0) + Xty / sigma2)
    beta = np.random.multivariate_normal(mu_n, Sigma_n)

    # 从 σ²|β, y 的条件后验采样（Inverse-Gamma）
    residuals = y_true - X @ beta
    a_n = a_0 + n / 2
    b_n = b_0 + residuals @ residuals / 2
    sigma2 = 1.0 / np.random.gamma(a_n, 1.0 / b_n)

    beta_samples[i] = beta
    sigma2_samples[i] = sigma2

# 结果汇总
burn_in_g = 2000
beta_post = beta_samples[burn_in_g:]
sigma2_post = sigma2_samples[burn_in_g:]

print("\n贝叶斯线性回归 Gibbs 采样结果 (后验均值 ± 标准差):")
print(f"  β₀: {beta_post[:,0].mean():.3f} ± {beta_post[:,0].std():.3f} "
      f"(真实: {beta_0_true})")
print(f"  β₁: {beta_post[:,1].mean():.3f} ± {beta_post[:,1].std():.3f} "
      f"(真实: {beta_1_true})")
print(f"  σ²: {sigma2_post.mean():.3f} ± {sigma2_post.std():.3f} "
      f"(真实: {sigma_true**2})")

PyMC贝叶斯推断

PyMC 是 Python 生态中最流行的概率编程库之一，它封装了底层的MCMC算法（如 NUTS 采样器），让用户可以专注于模型构建而非采样算法细节。在实际数据分析工作中，通常使用 PyMC、Stan 等高级工具来进行贝叶斯推断。

# PyMC 贝叶斯线性回归示例
# 注意：需要先安装 pymc: pip install pymc
import pymc as pm
import arviz as az

# 生成模拟数据
np.random.seed(42)
n = 100
x = np.random.normal(0, 1, n)
true_intercept = 1.2
true_slope = 0.8
true_sigma = 0.5
y = true_intercept + true_slope * x + np.random.normal(0, true_sigma, n)

# 定义模型
with pm.Model() as linear_model:
    # 先验分布
    intercept = pm.Normal('intercept', mu=0, sigma=10)
    slope = pm.Normal('slope', mu=0, sigma=10)
    sigma = pm.HalfNormal('sigma', sigma=1)

    # 似然函数
    mu = intercept + slope * x
    likelihood = pm.Normal('y', mu=mu, sigma=sigma, observed=y)

    # MCMC 采样（自动使用 NUTS 采样器）
    trace = pm.sample(draws=2000, tune=1000,
                      chains=4, random_seed=42)

# 结果摘要
summary = az.summary(trace, hdi_prob=0.95)
print(summary)

# 后验预测检查
posterior_predictive = pm.sample_posterior_predictive(
    trace, model=linear_model, random_seed=42
)

# 可视化
az.plot_trace(trace)
plt.suptitle('PyMC 贝叶斯线性回归 - 后验轨迹与分布')
plt.tight_layout()
plt.show()

az.plot_ppc(posterior_predictive, data_pairs={'y': y})
plt.title('后验预测检查')
plt.show()

PyMC 与 MCMC 的关键概念：

链（Chain）：通常运行4条独立的MCMC链，用于检查收敛性
R-hat 统计量：衡量链间方差与链内方差的比率，接近1.0表示收敛
有效样本量（ESS）：考虑自相关后的等效独立样本数
后验预测检查（PPC）：从后验分布模拟新数据，与观测数据对比验证模型
HDI（Highest Density Interval）：后验分布的最高密度区间，贝叶斯版本的置信区间

PyMC 的 NUTS（No-U-Turn Sampler）是 HMC（Hamiltonian Monte Carlo）的一种自适应变体，利用梯度信息高效探索高维后验空间，比 Metropolis-Hastings 和 Gibbs 采样效率高得多，尤其适合参数维度较高的模型。

核心要点总结

蒙特卡洛基本原理：通过大量随机采样来近似计算复杂数值问题，核心理论支撑为大数定律（保证收敛）和中心极限定理（提供 O(1/√n) 的误差估计）。收敛速度与问题维度无关，特别适合高维问题。
经典应用：计算π、数值积分（高维积分优势尤为显著）、Buffon投针、随机游走与布朗运动模拟、Bootstrap置信区间估计。
金融风险模拟：几何布朗运动是资产价格建模的基石；VaR 和 CVaR 是衡量金融风险的核心指标；蒙特卡洛可灵活处理多资产投资组合的相关性和非线性风险。
统计推断：Bootstrap（重采样估计统计量分布）和置换检验（非参数假设检验）是现代统计学的两大重采样方法，不依赖分布假设。
方差缩减：对偶变量、控制变量、重要抽样、分层抽样等技术可显著降低估计方差，将等效样本量提升数倍至数十倍。重要性抽样在稀有事件模拟中至关重要。
MCMC方法：Metropolis-Hastings是最基础的MCMC算法；Gibbs采样适用于条件分布已知的可分解模型；PyMC/Stan等概率编程框架通过NUTS等高级采样器使贝叶斯建模变得实用。
实践要点：MCMC中需注意 burn-in 丢弃、链收敛诊断（R-hat）、有效样本量（ESS）和自相关分析。方差缩减技术应作为蒙特卡洛模拟的标准实践而非事后优化。

附录：常用参考与进一步学习

库名	用途	安装
NumPy	随机数生成、基础蒙特卡洛	pip install numpy
SciPy	统计分布、积分、优化	pip install scipy
PyMC	概率编程、MCMC贝叶斯推断	pip install pymc
ArviZ	贝叶斯模型诊断与可视化	pip install arviz
emcee	MCMC采样器（适用于天体物理等）	pip install emcee
Pyro	深度概率编程（PyTorch 后端）	pip install pyro-ppl

收敛诊断的 R-hat 统计量

def r_hat(chains):
    """计算多链的 R-hat 统计量（Gelman-Rubin 诊断）

    参数:
        chains: 形状为 (n_chains, n_samples) 的数组
    返回:
        R-hat 值（接近1表示收敛，通常阈值 < 1.01）
    """
    n_chains, n_samples = chains.shape

    # 链内方差
    chain_means = chains.mean(axis=1)
    chain_vars = chains.var(axis=1, ddof=1)
    W = chain_vars.mean()  # 链内方差的均值

    # 链间方差
    B = n_samples * chain_means.var(ddof=1)

    # 边缘后验方差估计
    var_plus = (n_samples - 1) / n_samples * W + B / n_samples

    # R-hat
    rhat = np.sqrt(var_plus / W)
    return rhat

# 模拟测试
np.random.seed(42)
# 已收敛的链（都从同一分布采样）
converged = np.random.normal(0, 1, size=(4, 5000))
print(f"已收敛 R-hat: {r_hat(converged):.4f}")

# 未收敛的链（从不同分布采样）
not_converged = np.array([
    np.random.normal(-2, 1, 5000),
    np.random.normal(2, 1, 5000),
    np.random.normal(0, 2, 5000),
    np.random.normal(1, 1, 5000),
])
print(f"未收敛 R-hat: {r_hat(not_converged):.4f} (应 > 1.01)")

实践建议：在贝叶斯分析中，一般要求 R-hat < 1.01 才认为链已收敛。如果 R-hat 过大，需要增加采样步数、调整建议分布，或重新参数化模型。同时，有效样本量（ESS）应至少达到几百才能保证后验分位数估计的可靠性。