第三讲：蒙特卡洛方法——大数定律、重要性采样、蒙特卡洛积分

好，我们进入正题。

蒙特卡洛方法，说白了就是用随机数去解决确定性问题。你可能会问：为什么非要用随机？

我刚开始做定位算法时，也这么想。后来遇到一个非高斯噪声的场景——激光雷达在走廊里被玻璃反射干扰，噪声分布根本不是高斯。卡尔曼滤波直接崩了。那时候我才意识到，蒙特卡洛方法不是「花架子」，它是处理复杂概率问题的利器。

3.1 大数定律：蒙特卡洛的数学根基

大数定律告诉我们：当样本数量足够多时，样本均值会收敛到期望值。

用公式表达就是：

设 X₁, X₂, ..., Xₙ 是独立同分布的随机变量
E[X] = μ
那么 (X₁ + X₂ + ... + Xₙ) / n → μ  (当 n → ∞)

嗯，这里要注意：大数定律不保证「每次都能收敛」，它保证的是「概率意义上的收敛」。我在项目中遇到过有人拿10个样本就敢说「均值就是期望」，结果定位偏差大到离谱。

核心要点：

  • 样本量越大，估计越准
  • 但收敛速度是 O(1/√n)——想提高一位精度，样本量要翻100倍
  • 大数定律不要求噪声是高斯分布——这正是它在非高斯场景下的优势

我的经验：

在实际SLAM系统中，我一般取500-2000个粒子。太少，方差太大；太多，计算扛不住。这个范围是我在多个项目里试出来的平衡点。

3.2 蒙特卡洛积分：用随机点算面积

蒙特卡洛积分的思路很直观：

你想算一个复杂函数的积分，但解析解太难求。那就随机撒点，统计落在函数曲线下方的比例，再乘以总面积。

数学形式：

∫ f(x) dx ≈ (1/N) * Σ f(xᵢ)  其中 xᵢ ~ p(x)

举个例子：算圆周率π

import random
import math

def estimate_pi(num_samples):
    inside = 0
    for _ in range(num_samples):
        x = random.uniform(-1, 1)
        y = random.uniform(-1, 1)
        if x*x + y*y <= 1:
            inside += 1
    return 4 * inside / num_samples

# 试试不同样本量
for n in [100, 1000, 10000, 100000]:
    pi_est = estimate_pi(n)
    print(f"n={n:6d}, π≈{pi_est:.4f}, 误差={abs(pi_est-math.pi):.4f}")

输出结果：

n=   100, π≈3.1600, 误差=0.0184
n=  1000, π≈3.1480, 误差=0.0064
n= 10000, π≈3.1416, 误差=0.0000
n=100000, π≈3.1417, 误差=0.0001

看到了吗？样本量从100到10000，精度提升了一个数量级。但再往上翻10倍，精度提升就不明显了。这就是蒙特卡洛方法的「瓶颈」——收敛慢。

避坑指南：

我曾经在项目中直接用均匀采样去算一个尖峰分布函数的积分，结果样本全落在平坦区域，峰值区域一个点都没采到。算出来的积分值比真实值小了一个数量级。这就是「采样分布」选错的下场。

3.3 重要性采样：让采样更聪明

为什么要重要性采样？

你想想看，如果目标分布 p(x) 在大部分区域值都很小，只有一小块区域值很大。你用均匀采样，大部分点都浪费在「不重要」的区域。

重要性采样的思路：找一个容易采样的分布 q(x)，让它在「重要区域」多采点，然后通过权重修正偏差。

公式：

E[f(x)] = ∫ f(x) p(x) dx
        = ∫ f(x) [p(x)/q(x)] q(x) dx
        ≈ (1/N) * Σ f(xᵢ) * w(xᵢ)

其中 w(xᵢ) = p(xᵢ) / q(xᵢ) 是重要性权重

代码示例：

import numpy as np

# 目标分布：标准正态分布（但假设我们不知道如何直接采样）
def target_pdf(x):
    return np.exp(-x**2/2) / np.sqrt(2*np.pi)

# 提议分布：均匀分布（容易采样）
def proposal_pdf(x):
    return 1/10 if -5 <= x <= 5 else 0

# 重要性采样
np.random.seed(42)
n_samples = 1000
samples = np.random.uniform(-5, 5, n_samples)
weights = target_pdf(samples) / proposal_pdf(samples)

# 估计期望值
f = lambda x: x**2  # 想算 E[X²]
estimate = np.mean(f(samples) * weights) / np.mean(weights)
print(f"重要性采样估计 E[X²] = {estimate:.3f}")
print(f"真实值 E[X²] = 1.000")

关键点：

  • 提议分布 q(x) 的尾部要比目标分布 p(x) 更「胖」——否则权重会爆炸
  • 权重归一化很重要：实际中常用 wᵢ = wᵢ / Σ wⱼ
  • 有效样本量 ESS = 1 / Σ wᵢ² —— 用来衡量采样质量

我的习惯：

每次做重要性采样，我都会先算一下ESS。如果ESS小于样本量的10%，说明提议分布选得太差了，得重新设计。这个习惯帮我避免了好几次「看似收敛实则发散」的坑。

3.4 蒙特卡洛方法在粒子滤波中的角色

现在我们把三块内容串起来：

概念 在粒子滤波中的作用
大数定律 保证粒子数足够多时，粒子集能近似真实后验分布
蒙特卡洛积分 用粒子加权和来近似状态估计（如位置均值）
重要性采样 从提议分布采样粒子，用权重修正偏差

说白了，粒子滤波就是：

  1. 用重要性采样生成粒子（预测步骤）
  2. 用观测数据更新权重（更新步骤）
  3. 用蒙特卡洛积分计算状态估计（输出步骤）
  4. 大数定律保证：粒子越多，估计越准

我记得第一次在非高斯噪声下跑通粒子滤波时，看到定位轨迹从「飘忽不定」变成「稳定跟踪」，那种感觉真的很爽。蒙特卡洛方法不是万能的，但在非高斯、非线性场景下，它确实是最实用的工具之一。

最后提醒：

蒙特卡洛方法不是银弹。如果提议分布选得不好，或者粒子数太少，结果可能比卡尔曼滤波还差。我见过有人拿50个粒子去跑高维SLAM，结果定位完全发散。嗯，这种坑踩过一次就记住了。

下一讲，我们会把这三块知识整合起来，正式进入粒子滤波的核心算法。到时候你会看到，这些看似抽象的数学概念，是如何在代码中落地生根的。

← 上一章 📖 返回目录 下一章 →