蒙特卡洛与马尔可夫链：从随机采样到复杂问题求解

前言

蒙特卡洛方法是一类基于随机采样进行数值计算的方法，常用于解决复杂的数学、物理、金融等问题。

蒙特卡洛是摩纳哥公国的一个行政区，也是其著名的赌博和旅游中心。蒙特卡洛以其豪华的赌场、美丽的海岸线和高档的生活方式而闻名于世。

蒙特卡洛方法的基本思想

蒙特卡洛方法依赖于随机采样 和统计估计 来解决计算问题。其核心思想是通过大量的随机样本来估算问题的解。一般步骤如下：

1. 定义问题 ：明确你要解决的问题，以及哪些变量是未知的或随机的。
2. 随机采样 ：根据问题中的随机变量生成大量的随机样本。
3. 统计估计 ：通过对这些随机样本的统计分析来估算问题的解。

示例说明

计算圆周率

问题描述：

假设我们想要估算圆周率 π 的值。我们可以利用一个几何方法，在一个单位正方形内随机撒点，通过计算落在内切圆中的点的比例来估算 π 的值。

具体步骤：

1. 定义几何模型 ：

   • 我们考虑一个边长为 1 的正方形，其内切圆的半径也为 1。
   • 圆的面积是 π × r²，由于 $r=1$，所以圆的面积是 $\pi$。
   • 正方形的面积是 1。

2. 随机采样 ：

   • 我们在正方形内随机生成大量的点，假设生成点的总数为 $N$。
   • 计算每个点到圆心的距离，判断该点是否落在圆内。如果点 $(x,y)$ 满足 $x^2+y^2\le 1$，则该点在圆内。

3. 统计比例 ：

   • 计算落在圆内的点的数量，记为 $N_{circle}$。
   • 圆的面积和正方形面积的比例应该接近于 $\frac{\pi }{4}$，因此有：$$\frac{N_{\text{circle}}}{N}\approx \frac{\pi }{4}$$
   • 由此可以估算 $\pi$ 的值：$$\pi \approx 4\times \frac{N_{\mathrm{circle}}}{N}$$

Python 代码

import random

def estimate_pi(num_samples):
    inside_circle = 0

    for _ in range(num_samples):
        # 生成随机点 (x, y)，范围在 [0, 1] 之间
        x = random.uniform(0, 1)
        y = random.uniform(0, 1)
        
        # 判断该点是否在单位圆内
        if x**2 + y**2 <= 1:
            inside_circle += 1

    # 根据公式估算π
    pi_estimate = 4 * (inside_circle / num_samples)
    return pi_estimate

# 使用1000000个随机样本来估算π
pi_approx = estimate_pi(1000000)
print(f"估算的π值: {pi_approx}")

计算积分

假设我们要计算函数 $f(x)=sin(x)$ 在区间 $[0,\pi ]$ 上的定积分。
问题描述：
定积分的数学表达式为：$$I={\int }_0^{\pi }\sin (x)dx$$

具体步骤

蒙特卡洛方法计算定积分的步骤：

1. 定义积分区间 ：我们知道积分的上下限是 $[0,\pi ]$。
2. 随机采样 ：在积分区间内生成大量的随机数。
3. 计算函数值 ：对于每个随机数 $x_i$​，计算 $f(x_i)$，即 $sin(x_i)$。
4. 估算积分值 ：[[#平均值和大数定律|通过这些随机样本的平均值乘以积分区间的长度来估算积分值]]。

根据蒙特卡洛积分的原理，积分 $I$ 可以近似为：$$I\approx \frac{\pi -0}{N}\sum _{i=1}^N\sin (x_i)$$
其中，$N$ 是随机样本的数量，$x_i$ 是在区间 $[0,\pi ]$ 内均匀分布的随机数。

Python 代码

import random
import math

def monte_carlo_integration(num_samples):
    sum_fx = 0.0
    
    for _ in range(num_samples):
        x = random.uniform(0, math.pi)  # 在 [0, π] 区间内生成随机数
        sum_fx += math.sin(x)            # 计算 sin(x) 并累加

    # 估算积分值
    integral_estimate = (math.pi - 0) / num_samples * sum_fx
    return integral_estimate

# 使用1000000个随机样本来估算积分
integral_approx = monte_carlo_integration(1000000)
print(f"蒙特卡洛方法估算的积分值: {integral_approx}")

总结

蒙特卡洛方法是一种强大的数值计算工具，通过随机采样和统计分析来解决复杂的数学和物理问题。其核心思想是通过大量随机样本来估算问题的解，尤其适用于高维问题和具有不确定性的系统。
优缺点：

• 优点 ：简单易实现，适用于高维问题，通用性强。
• 缺点 ：收敛速度慢，存在随机误差，不适合所有问题。

你可以把蒙特卡洛看成 “随机化的穷举法”

补充说明

平均值和大数定律

假设我们在区间 $[a,b]$ 内生成了 $N$ 个随机点 $x_1,x_2,\dots ,x_N$，这些点均匀分布在 $[a,b]$ 上。对于每个点 $x_i$，我们计算其函数值 $f(x_i)$。

根据概率论的大数定律，当 $N$ 足够大时，这些随机点上的函数值的平均值 $\overline{f}$ 会接近于函数 $f(x)$ 在区间 $[a,b]$ 上的平均值：
$$\overline{f}=\frac{1}{N}\sum _{i=1}^{N}f(x_i)\approx \frac{1}{b-a}{\int }_a^{b}f(x)dx$$
因此，函数 $f(x)$ 在区间 $[a,b]$ 上的平均值可以近似表示为：
$$\text{平均值}(f)\approx \frac{1}{N}\sum _{i=1}^{N}f(x_i)$$

马尔可夫蒙特卡洛

马尔可夫蒙特卡洛（Markov Chain Monte Carlo，MCMC）是一类结合了马尔可夫链和蒙特卡洛方法的算法，主要用于从复杂概率分布中进行采样，或估计积分和优化问题。

核心思想

MCMC方法的核心思想是通过构建一个马尔可夫链，使其平稳分布（stationary distribution）与目标概率分布一致。然后通过蒙特卡洛模拟，即从该马尔可夫链中进行采样，从而近似得到目标分布的样本。

基本步骤

1. 定义目标分布 ：首先，定义你想要从中采样的目标概率分布，通常是高维的且难以直接采样。

2. 构建马尔可夫链 ：设计一个马尔可夫链，其状态空间对应于目标分布的定义域，并使其转移概率满足细致平衡条件（detailed balance），从而保证马尔可夫链的平稳分布是目标分布。

3. 采样过程 ：从某个初始状态开始，按照马尔可夫链的转移概率进行状态转移，生成一系列样本。这些样本在经过一定的“[[#烧入期]]”（burn-in period）后，可以近似认为是来自目标分布的独立样本。

4. 估计与推断 ：利用生成的样本来估计目标分布的特性，比如期望、方差，或者进行贝叶斯推断。

总结

总结来说，马尔科夫蒙特卡洛是对于目标状态我们没有办法直接采样，所以我们才用马尔可夫链来一直运行使用蒙特卡洛来采用实现完成随机采样，在这个过程中，对每次采样进行判断是否采用采样。

1. 无法直接采样 ：对于复杂的目标分布 p(x)，我们通常无法直接从中进行采样。这些分布可能非常复杂，比如是高维的、多峰的、或者没有已知的解析形式，这使得传统的采样方法无法直接应用。

2. 构建马尔可夫链 ：为了处理这个问题，MCMC方法通过构建一个马尔可夫链，其平稳分布（stationary distribution）是我们的目标分布 p(x)。马尔可夫链的每一个状态代表一个可能的样本点。

3. 蒙特卡洛采样 ：通过使用蒙特卡洛方法，我们从马尔可夫链中进行采样。具体来说，我们从一个初始状态开始，按照马尔可夫链的转移概率进行状态转移，生成一系列样本。

4. 判断是否接受采样 ：在每一步生成候选样本后，MCMC算法（例如Metropolis-Hastings）会通过计算接受概率来决定是否接受该候选样本。这个过程确保了生成的样本逐渐符合目标分布的特性。

可以这样总结整个过程：

1. 初始化：选择一个初始状态 $x_0$。

2. 提议新样本：从提议分布 $q(x^{\prime }\mid x_t)$ 中生成一个候选样本 $x^{\prime }$。

3. 计算接受概率：根据目标分布 $p(x)$ 和提议分布 $q$，计算接受概率 $\alpha$。

4. 决定是否接受：

   • 如果 $\alpha \ge 1$，则接受候选样本，设置 $x_{t+1}=x^{\prime }$。
   • 如果 $\alpha <1$，则根据接受概率 $\alpha$ 决定是否接受候选样本。具体来说，从均匀分布 $U(0,1)$ 中生成一个随机数 $u$，如果 $u\le \alpha$，则接受候选样本，否则拒绝候选样本，保持 $x_{t+1}=x_t$。

5. 重复：重复上述步骤，生成一系列样本 { x t } \{x_t\} {xt​}。

烧入期

烧入期（Burn-in Period）指的是在从马尔可夫链中收集样本的初期阶段，这段时间内的样本会被丢弃，不纳入后续的统计分析。

当我们开始运行MCMC算法时，初始状态可能离目标分布较远，因此在最开始的若干步中，生成的样本可能并不能很好地代表目标分布的特性。这是因为在MCMC的早期阶段，马尔可夫链还处于“探索”阶段，尚未达到其平稳分布（stationary distribution）。

通过指定一个烧入期，我们可以在开始采样之前，让马尔可夫链有时间“稳定”到目标分布。在此之后，我们才开始收集样本用于后续的分析和推断。

End

关注唯一微信公众号：码上地球！🌹

参考文献

• 👉 蒙特卡洛