蒙特卡洛方法¶

前言

起源：曼哈顿计划
- 计算中子在非均匀介质中的散射和运输问题
Metropolis, Ulam

Note

随机试验
- 可以在相同条件下重复进行
- 每次试验的结果不止一个，但事先知道所有可能的结果
- 每次试验的结果不确定
样本点 \(e\)：随机试验的每个可能结果
样本空间 \(S\)：随机试验的所有可能结果的集合
随机事件 \(A\)：样本空间 \(S\) 的子集

对于样本空间的子集（随机事件） \(A\)，赋予一个实数 \(P(A)\)，称为事件 \(A\) 的概率，满足

非负性：\(P(A) \ge 0\)
规范性：\(P(S) = 1\)
可列可加性：对于互补相容的事件列 \(A_1, A_2, \ldots\)，\(A_i \cap A_j = \varnothing\)，有

\[ P\left(\bigcup_{i=1}^{n} A_i\right) = \sum_{i=1}^{n} P(A_i) \]

全概率公式
贝叶斯公式
随机变量：将样本空间 \(S\) 的每个样本点 \(e\) 单值映射到实数域的函数 \(X = X(e)\)，称 \(X\) 为随机变量

例

3MeV 的光子入射到铅板中，可能会发生：

光电效应
康普顿散射
正负电子对产生

分别赋予离散值 \(X = 1,2,3\)，相应的定义不同随机变量值下的概率 \(P(X = 1), P(X = 2), P(X = 3)\).

随机数产生¶

真随机与准随机¶

真随机过程：微观物理过程
- 量子尺度下，不可预测每次事件的结果
- 不稳定同位素的衰变
- 粒子间碰撞产生的末态
准随机过程：可以用经典的物理定律描述，但是对于初始值的微小差异及外界干扰非常敏感，以至于无法精确预测
- 混沌系统
- 大气环流
- 布朗运动

伪随机数¶

真随机过程可以产生随机数，但是无法快速产生大量随机数，与程序接口难以匹配。

伪随机数：通过确定性算法产生的数列，具有类似随机数的统计特性。

要求：

长的循环周期
良好的随机性
- 随机数序列 \(\{x_n\}\) 中的所有随机数之间没有大的关联，或自相关函数 \(\rho(k)\) 快速衰减
速度快

线性同余法¶

LCG（Linear Congruential Generator）：用于生成 \([0,1]\) 之间的均匀分布的伪随机数。

\[ \begin{aligned} X_{n+1} &= (aX_n + c) \mod m \\ \xi_n &= \frac{X_n}{m} \end{aligned} \]

\(a\)：乘子
\(c\)：增量
\(m\)：模数
初始值 \(X_0\)：种子

\(X_n = 0,1,2,\ldots,m-1\)，故 \(\xi_n < 1\)。

Remark

LCG 的周期最大为 \(m\)，但大部分情况会小于 \(m\)。为了使 LCG 达到最大周期和更好的随机性：

\(m\) 足够大
\((c,m) = 1\)
合理选择乘子 \(a\)，例如 \(m\) 的所有质因数 \(p\,|\,a-1.\)

Python 实现

import numpy as np
import matplotlib.pyplot as plt

rnd1 = []
rnd2 = []

# LCG parameters
a = 69069
c = 1
m = 419430409200233

# seed
rand = 100
for i in range(10000):
    rand = (a * rand + c) % m
    if i % 2 == 0:
        rnd1.append(rand / m)
    else:
        rnd2.append(rand / m)

plt.scatter(rnd1, rnd2, s=0.05)
plt.show()

任意分布的抽样¶

离散随机变量¶

离散随机变量 \(X\) 的概率分布：\(P(X = x_i) = p_i, \, i = 1,2,\ldots,n.\)
定义累积概率 \(\xi\)

应用¶

模拟退火算法¶

如何快速找到函数的全局最小值？

模拟退火算法（Simulated Annealing） 借鉴了 Metropolis-Hastings 算法和物理中固体物质的退火过程。

如果有一个物理系统处在热平衡态，各状态的分布符合玻尔兹曼分布：

\[ P(E_1) = \frac{e^{-\beta E_1}}{Z} \]

系统有一个最低能级，不妨设 \(E_0 = 0\)，其他能级 \(E_i > 0\)。假设系统冷却到绝对零度 \(T \to 0\)，那么除了最低能级外，其他能级的概率 \(e^{-\beta E_i} \to 0\)，所以