蒙特卡洛方法初步

1. 蒙特卡洛方法概述

📖 核心思想

通过统计采样（随机抽样）解决传统方法难以处理的数学/物理问题——积分、优化、多自由度系统模拟等。

🕰 历史起源

时间	人物	贡献
1946	Ulam, von Neumann, Metropolis	曼哈顿计划中发明，分析中子输运
后	各领域	统计物理、分子模拟、量子多体、粒子物理…

💡 物理直觉

中子与原子核作用受量子力学制约——只能知道发生概率，无法精确获得位置、速度、方向。用随机抽样模拟大量中子的统计行为，就能推知宏观传输性质。大量随机→确定结论。

🔄 方法流程

产生随机数 → 按概率分布采样 → 模拟/计算 → 统计平均 → 结果

2. 概率统计基础

2.1 概率密度函数（PDF）

概念	记号	说人话
概率密度函数	$p(x)$	随机变量取某值的"相对可能性"
累积分布函数	$P(x) = \int_{-\infty}^{x} p(t) dt$	从 $-\infty$ 积分到 $x$ 的累计概率
期望值	$\langle f\rangle = \int f(x)p(x)dx$	按概率加权平均
$n$ 阶中心矩	$\mu_n = \langle(x-\langle x\rangle)^n\rangle$	描述分布形状（2阶=方差）

💡 数值例子

某物理量 $x$ 服从 $p(x) = 2x, \; x\in[0,1]$ （归一化）：

$\langle x\rangle = \int_0^1 x \cdot 2x \, dx = \int_0^1 2x^2 dx = \frac{2}{3} \approx 0.667$
$\sigma^2 = \langle x^2\rangle - \langle x\rangle^2 = \int_0^1 2x^3 dx - (\frac{2}{3})^2 = \frac{1}{2} - \frac{4}{9} = \frac{1}{18} \approx 0.056$

2.2 中心极限定理 ⭐

不管原分布 $p(x)$ 长什么样，$m$ 个独立样本的均值 $\bar{x}$ 的分布 → 正态分布。

\bar{x} \sim \mathcal{N}\!\left(\mu,\; \frac{\sigma^2}{m}\right)

$\mu$ ：原分布的均值
$\sigma^2/m$ ：均值的方差 → 样本越多，精度越高

💡 说人话

你测量某个物理量 100 次求平均，这个"平均值"的误差是单次测量误差的 $1/\sqrt{100}=1/10$ 。样本量翻 4 倍，精度翻 2 倍。

3. 随机数生成

3.1 伪随机数要求

要求	说明
均匀性	在 $[0,1]$ 上均匀分布
去关联	相邻随机数之间相关性可忽略
长周期	重复之前能产生的序列尽可能长
速度快	算法效率高

💡 实践建议：使用系统自带随机数生成器（如 numpy.random），不自己造轮子。

3.2 高斯分布采样：Box-Muller 法

把两个在（0，1）上独立均匀分布的随机数（u1， u2），通过极坐标变换，变成两个独立的标准正态分布随机数！

因为电脑自带的随机数生成器/编程中，只能生成均匀分布的随机数，每个数字出现的概率一样大。但我们这里Guass想要得到的是标准正态分布的随机数，假设为 X和Y，有已知的pdf公式。

假设我们要同时生成两个独立的正态分布随机数 $X$ 和 $Y$ 。把它们的概率公式乘起来，看看它们作为一个整体（在二维平面上的一个点 $(X,Y)$ ）的概率是多少：

f(x, y) = f(x) \times f(y) = \frac{1}{2\pi} e^{-\frac{x^2 + y^2}{2}}

有平方和，自然想到极坐标来表示可能会简化这个事情

\frac{1}{2\pi} e^{-\frac{x^2 + y^2}{2}} dx dy = \frac{1}{2\pi} e^{-\frac{r^2}{2}} r dr d \theta

f(\theta) = \frac{1}{2\pi} ,\space \theta \in (0, 2\pi) \\ \space \\ f(r) = r e^{-\frac{r^2}{2}} , \space r \in (0, +\infty)

二维分布被拆解成了两个完全独立的部分可以尝试用两个独立的均匀分布随机数 u1，u2 来生成它们

选择用 “逆变换采样法” 生成 $\theta$ 和 r

逆变换采样法的原理：如果知道一个随机变量的累积分布函数 CDF为 F（x），那么令 $F(x) = u$ （其中 $u$ 是 $(0,1)$ 的均匀分布随机数），反解出 $x = F^{-1}(u)$ ，就能得到符合该分布的随机数。

就是说：假设你要生成的随机变量为 X，其CDF为 F（x）（严格单调递增）

逆变换法说：令 U = F（x），那么 U 一定服从（0， 1）上的均匀分布！

proof：

对于任意 $u∈(0,1)$ ，看 U 的累积分布函数：

P(U≤u)=P(F(X)≤u)

因为 F 是单调递增函数，不等式 $F(X)≤u$ 两边同时作用反函数 $F^{-1}$ ，等价于：

P(X≤F−1(u))

根据分布函数 F的定义，

P(X≤F^{−1}(u))=F(F^{−1}(u))=u

结论出来了：

P(U≤u)=u

这正是 (0,1)均匀分布的累积分布函数！

既然 U=F(X) 服从均匀分布，那么反过来： 如果我从均匀分布中取一个随机数 u，令 $x=F^{-1}(u)$ ，那么 x的分布函数一定就是 F(x)。

因为：

P(X≤x)=P(F^{−1}(U)≤x)=P(U≤F(x))=F(x)

第一部分：角度 $\theta$ 在 $0$ 到 $2\pi$ （即 360 度）之间是完全均匀分布的。那怎么用均匀随机数 $u_2$ 来生成角度呢？很简单，直接按比例放大：

它的概率密度是 $f(\theta) = \frac{1}{2\pi}$ 。求累积分布函数（CDF）：

F(\theta) = \int_{0}^{\theta} \frac{1}{2\pi} dt = \frac{\theta}{2\pi}

令 $F(\theta) = u_2$ ：

\frac{\theta}{2\pi} = u_2 \implies \theta = 2\pi u_2

第二部分：距离中心的半径 $R$ 点落在距离圆心为 $R$ 的圆环上的概率，经过微积分推导，满足一个叫做“指数分布”的规律。数学上有一个很成熟的方法（叫做逆变换采样），把均匀随机数 $u_1$ 变成指数分布。

它的概率密度是 $f(r) = r e^{-\frac{r^2}{2}}$ 。求累积分布函数（CDF）：

F(r) = \int_{0}^{r} s e^{-\frac{s^2}{2}} ds

为了积这个分，我们令 $v = \frac{s^2}{2}$ ，那么 $dv = s ds$ ：

F(r) = \int_{0}^{\frac{r^2}{2}} e^{-v} dv = \left[ -e^{-v} \right]_{0}^{\frac{r^2}{2}} = 1 - e^{-\frac{r^2}{2}}

令 $F(r)$ 等于一个均匀随机数。这里有个小技巧：如果 $u_1$ 在 $(0,1)$ 上均匀分布，那么 $1 - u_1$ 同样在 $(0,1)$ 上均匀分布。为了方便，我们直接令 $1 - e^{-\frac{r^2}{2}} = 1 - u_1$ （相当于令 $e^{-\frac{r^2}{2}} = u_1$ ）：

e^{-\frac{r^2}{2}} = u_1

两边取自然对数：

-\frac{r^2}{2} = \ln(u_1) \\ \space \\ r = \sqrt{-2 \ln(u_{1})}

现在我们已经用 $u_1$ 和 $u_2$ 求出了 $r$ 和 $\theta$ ，最后一步就是用极坐标定义把它塞回最初的 $X$ 和 $Y$ 里：

X = r \cos\theta = \sqrt{-2 \ln(u_1)} \cos(2\pi u_2) \\ \space \\ Y = r \sin\theta = \sqrt{-2 \ln(u_1)} \sin(2\pi u_2)

推导完成。通过这套严密的微积分坐标变换和积分反解，我们证明了只要输入两个均匀分布的 $u_1$ 和 $u_2$ ，输出的 $X$ 和 $Y$ 就必定是严格独立的标准正态分布。

import numpy as np

def box_muller(n):
    """产生 n 个标准正态分布随机数"""
    u1 = np.random.uniform(0, 1, n)
    u2 = np.random.uniform(0, 1, n)
    r = np.sqrt(-2 * np.log(u1))
    theta = 2 * np.pi * u2
    return r * np.cos(theta)  # 另一组为 r * np.sin(theta)

💡 直观验证

对生成的 $10^4$ 个随机数：均值 ≈ 0，标准差 ≈ 1。画直方图应呈钟形曲线。

4. 蒙特卡洛积分

核心：把积分写成期望值，再用随机样本平均来估计期望值！！！

4.1 空间随机投点法（打靶法）

计算 $\displaystyle I = \int_a^b f(x)dx$ ，用矩形区域包住函数曲线，然后"撒点"。
假设有 f(x) 是正定的，且大于0的同时有最大值 fmax
其积分就是曲线下方的面积！！！

用矩阵包住曲线，于是有 $\text{矩形面积}=(b-a)f_{\max}$

在矩形中随机撒点：

x_i\sim U(a,b) \\ \space \\ y_i\sim U(0,f_{\max})

如果：

y_i\le f(x_i)

说明点落在曲线下方。

设总点数 $N$ ，落在曲线下方点数 $n$ ，则：

\frac{n}{N} \approx \frac{\text{曲线下面积}}{\text{矩形面积}}

所以：

\frac{n}{N} \approx \frac{I}{(b-a)f_{\max}}

因此：

\boxed{ I\approx \frac{n}{N}(b-a)f_{\max} }

此估计是无偏估计！！！

┌─────────────────────┐
│    ○  ×  ○           │  f_max ↑
│  ○  ×  ×  ○  ×      │
│○  ×  ×  ○  ×  ×  ○ │          f(x) 曲线
│  ×  ×  ×  ○  ×  ×    │
│○  ×  ×  ○  ×  ×  ○ │
└─────────────────────┘
   a                  b

缺点在于如果f（x）大部分区域很小，而 fmax很大，那么会造成随机撒的点中，很多点都浪费了，虽然这个方法很直观，但是效率通常是不如期望值法！

4.2 期望值法 ⭐

I = \int_a^b f(x)dx \approx \frac{b-a}{N} \sum_{i=1}^{N} f(x_i), \quad x_i \in [a,b] \text{ 均匀采样}

方差估计

定义：

\bar f = \frac1N\sum_{i=1}^N f(x_i)

估计量：

I_N=(b-a)\bar f

假设样本独立同分布，则：

\operatorname{Var}(\bar f) = \frac{\operatorname{Var}(f)}{N}

其中：

\operatorname{Var}(f) = \langle f^2\rangle-\langle f\rangle^2

所以：

\operatorname{Var}(I_N) = (b-a)^2\operatorname{Var}(\bar f)

其中：

\sigma_f^2 = \frac1N\sum_{i=1}^N f(x_i)^2 - \left( \frac1N\sum_{i=1}^N f(x_i) \right)^2

误差：

\boxed{ \sigma_I = (b-a)\frac{\sigma_f}{\sqrt N} }

蒙特卡洛误差的核心结论：

误差与 N的平方根成反比！！！

🔢 MC 算 π 例子

import numpy as np

N = 100000
x = np.random.uniform(-1, 1, N)
y = np.random.uniform(-1, 1, N)
inside = (x**2 + y**2 <= 1).sum()
pi_est = 4 * inside / N
sigma = np.sqrt(pi_est * (4 - pi_est) / N)  # 二项式方差

print(f"π ≈ {pi_est:.5f} ± {sigma:.5f}")
# 典型输出: π ≈ 3.14159 ± 0.00518

4.3 多维积分

MC 积分的精度与维度无关——这是最大的优势！

考虑 $d$ 维积分：

I= \int_\Omega f(\mathbf x)\,d^d x

若区域体积是：

V_\Omega

均匀采样：

\mathbf x_i\sim U(\Omega)

则：

I = V_\Omega \langle f\rangle

估计：

\boxed{ I_N= \frac{V_\Omega}{N} \sum_{i=1}^N f(\mathbf x_i) }

方差：

\boxed{ \sigma_I^2 = \frac{V_\Omega^2}{N} \left( \langle f^2\rangle-\langle f\rangle^2 \right) }

关键是：

\sigma_I\propto \frac1{\sqrt N}

这个 $1/\sqrt N$ 不显式依赖维度 $d$ 。

定积分方法	误差 vs 维度
梯形法/辛普森法	误差 ∝ $N^{-k/d}$ （$d$=维度，指数衰减！）
蒙特卡洛法	误差 ∝ $1/\sqrt{N}$ （与 $d$ 无关）

💡 说人话

算 100 维积分时，网格法需要天文数字的网格点，因为维度d越大，收敛越慢。
MC 只需要足够多样本——维度灾难的克星。

5. 重要性采样（Importance Sampling）

5.1 动机

均匀采样低效——被积函数 $f(x)$ 在大部分区域接近零时，大量采样点"浪费"了。
在大的地方多采样，在小的地方少采样！！！

5.2 原理

引入与 $f(x)$ 形状相似的采样函数 $F(x)$ ，要求$F(x)$ 在 $f(x)\ne0$ 的地方不能为0：

I = \int \frac{f(x)}{F(x)} \cdot F(x) dx = \left\langle \frac{f}{F} \right\rangle_F

离散形式 ：

I \approx \frac{1}{N}\sum_{i=1}^{N} \frac{f(x_i)}{F(x_i)}, \quad x_i \sim F(x)

估计量 p就是上文中的F：

Y(x)=\frac{f(x)}{p(x)}

那么：

I_N=\frac1N\sum_i Y(x_i)

方差：

\operatorname{Var}(I_N) = \frac1N \left[ \langle Y^2\rangle_p-\langle Y\rangle_p^2 \right]

其中：

\langle Y\rangle_p=I

而：

\langle Y^2\rangle_p = \int \left( \frac{f(x)}{p(x)} \right)^2p(x)\,dx

因此：

\boxed{ \operatorname{Var}(I_N) = \frac1N \left[ \int\frac{f(x)^2}{p(x)}dx - I^2 \right] }

这个公式告诉我们：

如果 $p(x)$ 在 $f(x)$ 大的地方也大，那么 $\frac{f^2}{p}$ 不会太大，方差就小。

最优重要性分布

如果 $f(x)\ge 0$ ，理论最优选择是：

p_{\text{opt}}(x)=\frac{f(x)}{\int f(x)dx}

因为此时：

\frac{f(x)}{p_{\text{opt}}(x)} = \int f(x)dx = I

每个样本都给出同样结果，方差为 0。

但问题是：

p_{\text{opt}}

需要知道 $I$ ，而 $I$ 正是我们要求的。

所以实际中选择一个和 $f(x)$ 形状相似、但容易采样的 $p(x)$ 。

💡 说人话

把"平均撒点"改成"在 $f(x)$ 大的地方多撒、小的地方少撒"——花同样的计算量，精度更高。

🔢 数值示例

求 $\displaystyle I = \int_0^3 x^2 e^{-x} dx$ ，已知 5 个 $[0,1]$ 随机数：0.033, 0.45, 0.51, 0.83, 0.95。

方法一：均匀采样

$x_i = a + (b-a)t_i = 3t_i$

$t_i$	$x_i$	$f(x_i) = x_i^2 e^{-x_i}$
0.033	0.099	0.009
0.450	1.350	0.472
0.510	1.530	0.508
0.830	2.490	0.514
0.950	2.850	0.470

I \approx \frac{3}{5}(0.009+0.472+0.508+0.514+0.470) = \frac{3}{5} \times 1.973 = 1.184

\sigma_f^2 = \langle f^2\rangle - \langle f\rangle^2 \approx 0.298 - 0.156 = 0.142

\sigma_I = (b-a)\sqrt{\frac{\sigma_f^2}{N}} = 3\sqrt{\frac{0.142}{5}} \approx 0.506

可以看出这个值本身接近精确值，但是误差估计很大因为样本量太少了！

方法二：重要性采样选取 $p(x) \propto e^{-x}$

因为被积分函数有类似 e的负x的形式分布，所以自然的为重要性采样选取

p(x) = C e^{-x} , x \in [0,3] \\ \space \\ 归一化后，可以求的 \space C = \frac{1}{1-e^{-3}}

现在：

\frac{f(x)}{p(x)} = \frac{x^2e^{-x}}{e^{-x}/(1-e^{-3})}

所以估计量变成：

\boxed{ I_N= \frac1N \sum_{i=1}^N (1-e^{-3})x_i^2 }

其中：

x_i\sim p(x)=\frac{e^{-x}}{1-e^{-3}}

这比直接算 $x^2e^{-x}$ 波动小很多。

如何从 $p(x)\propto e^{-x}$ 抽样？

CDF：

F(x)=\int_0^x\frac{e^{-t}}{1-e^{-3}}dt

令：

u=F(x)

即：

u= \frac{1-e^{-x}}{1-e^{-3}}

于是：

u(1-e^{-3})=1-e^{-x}

取负对数：

\boxed{ x= -\ln\left[ 1-u(1-e^{-3}) \right] }

这才是标准的反函数采样公式, u是[0,1]的均匀随机分布！

📊 对比

方法	$N$	估值	方差 $\sigma_I$	精度
均匀采样	5	1.184	0.506	低
重要性采样	5	≈ 1.155	远小于均匀	高

6. 随机行走与马尔可夫链

6.1 一维随机行走

每一步以概率 $R$ 向右、$L$ 向左，步长 $\Delta x = l$ , R + L = 1

第 $k$ 步的随机变量：$\xi_k = \begin{cases} +l & \text{概率 } R \ -l & \text{概率 } L \end{cases}$
$n$ 步后位置：$Xn = \sum{k=1}^{n} \xi_k$

\langle X_n\rangle = n \cdot [R \cdot l + L \cdot (-l)] = nl(R-L)

\sigma^2(X_n) = 4nRL \cdot l^2

特殊的，如果 R = L = 0.5

标准差：

\sigma(X_n)=l\sqrt n

这就是随机行走的典型结论：

位移平均可能为 0，但扩散宽度随 $\sqrt n$ 增长。

💡 物理直觉

设每步时间为：

\tau

则：

t=n\tau

无偏随机行走方差：

\langle X^2\rangle=nl^2

写成时间：

\langle X^2\rangle=\frac{t}{\tau}l^2

扩散方程中一维均方位移：

\langle X^2\rangle=2Dt

比较：

2Dt=\frac{l^2}{\tau}t

所以：

\boxed{ D=\frac{l^2}{2\tau} }

这就是扩散系数。

如果有偏随机行走，还会出现漂移速度：

v=\frac{l(R-L)}{\tau}

大量粒子独立随机行走 → 粒子分布演化为扩散方程： $$\frac{\partial P(x,t)}{\partial t} = D\frac{\partial^2 P(x,t)}{\partial x^2}$$

其中扩散系数 $D = \dfrac{l^2}{2\tau}$ ，$\tau$ 为每步时间。

6.2 马尔可夫过程

定义：下一步只取决于当前状态，与历史无关。行走路径 = 马尔可夫链。

马尔可夫过程的定义：

P(X_{n+1}=j|X_n=i,X_{n-1},\dots,X_0) = P(X_{n+1}=j|X_n=i)

意思是：

下一步只依赖当前状态，不依赖过去历史。

随机行走就是典型马尔可夫链。

转移概率矩阵

设状态概率向量： $$ \mathbf w(t)= \begin{pmatrix} w1(t)\ w2(t)\ w3(t)\ w4(t) \end{pmatrix} $$ 其中： $$ wi(t)=P(\text{系统在状态 }i) $$ 转移矩阵 $W$ 定义为： $$ W{ij}=P(j\to i) $$ 注意你的讲义采用的是：

第 $j$ 列表示从状态 $j$ 出发，转移到各状态的概率。

因此每列之和为 1： $$ \sumi W{ij}=1 $$

W = \begin{pmatrix} 1-R & L & 0 & 0 \\ R & 0 & L & 0 \\ 0 & R & 0 & L \\ 0 & 0 & R & 1-L \end{pmatrix}

$W_{ij}$ ：从状态 $j$ 转移到状态 $i$ 的概率
每列之和 = 1（归一化）

时间演化（矩阵形式）

\mathbf{w}(t+\Delta t) = W \cdot \mathbf{w}(t)

\mathbf{w}(n\Delta t) = W^n \cdot \mathbf{w}(0)

🔢 数值例子

4 状态系统，$L=R=0.5$，初始 $\mathbf{w}(0) = (1,0,0,0)^T$ ：

t=0: [1.000, 0.000, 0.000, 0.000]
t=1: [0.500, 0.500, 0.000, 0.000]
t=2: [0.500, 0.250, 0.250, 0.000]
t=3: [0.500, 0.125, 0.250, 0.125]
...
平衡: (4/15, 8/15, 1/5) ≈ (0.267, 0.533, 0.200)  (3态例)

💡 自然结论

如果经过很久后分布不再变化：

\mathbf w_* = W\mathbf w_*

这说明 $\mathbf w_*$ 是 $W$ 的本征矢量，对应本征值：

\lambda=1

所以：

\boxed{ W\mathbf w_*=\mathbf w_* }

这就是平衡分布。

Perron-Frobenius 定理告诉我们，在合适条件下：

转移矩阵非负；
马尔可夫链不可约；
非周期；

则长期分布会收敛到唯一平衡分布。

马尔可夫链演化到最后,分布 $\mathbf{w}(t)$ 收敛到转移矩阵 $W$ 的最大本征值对应的本征矢量（Perron-Frobenius 定理）。

7. Metropolis 算法与 MCMC

7.1 平衡态与细致平衡

目标：从 Boltzmann 分布 $p(E) \propto e^{-\beta E}$ 采样。

我们想从某个复杂分布采样：

p_i

在统计物理中最常见的是 Boltzmann 分布：

p_i=\frac{e^{-\beta E_i}}{Z}

其中：

Z=\sum_i e^{-\beta E_i}

是配分函数。

问题：

$Z$ 很难算，但 Metropolis 只需要概率比值。

因为：

\frac{p_j}{p_i} = \frac{e^{-\beta E_j}/Z}{e^{-\beta E_i}/Z} \\ \space \\ = e^{-\beta(E_j-E_i)} \\ \space \\ = e^{-\beta\Delta E}

配分函数 $Z$ 抵消了。

系统趋于平衡态的条件——细致平衡（Detailed Balance）： $$ W{i\to j}, pi = W{j\to i}, pj $$

意思是：$i \to j$ 的流量 = $j \to i$ 的流量 → 分布不再变化。

总流入状态 $j$ 的概率：

p_j'=\sum_i p_i W_{i\to j}

如果细致平衡成立：

p_iW_{i\to j}=p_jW_{j\to i}

所以：

p_j' = \sum_i p_j W_{j\to i}

从状态 $j$ 出发到所有状态的概率和为 1：

\sum_i W_{j\to i}=1

所以：

p_j'=p_j

因此分布稳定。

💡 说人话

好比教室里的座位：有人从 A 换到 B，有人从 B 换到 C ……当每个方向的人流相等时，座位分布稳定了。这就是"热平衡"的统计图像。

7.2 Metropolis 算法 ⭐（20世纪十大算法之一）

转移概率可以写成：

W_{i\to j}=T_{i\to j}A_{i\to j}

其中：

$T_{i\to j}$ ：提议概率；
$A_{i\to j}$ ：接受概率。

标准 Metropolis 假设提议对称：

T_{i\to j}=T_{j\to i}

细致平衡要求：

p_iT_{i\to j}A_{i\to j} = p_jT_{j\to i}A_{j\to i}

提议概率抵消：

p_iA_{i\to j}=p_jA_{j\to i}

所以：

\frac{A_{i\to j}}{A_{j\to i}} = \frac{p_j}{p_i}

Metropolis 选择，在上述中选择接受概率最大的那个，因为接受概率越大，链的“惰性”越小，随机游走探索的越快，收敛效率最高！概率的上线又是1，所以就把一个方向的概率顶到天花板1！

\boxed{ A_{i\to j} = \min\left(1,\frac{p_j}{p_i}\right) }

这就是在满足物理约束下，最高效的走法！并不是“随便取的”，而是最优解

对于 Boltzmann 分布：

\frac{p_j}{p_i}=e^{-\beta(E_j-E_i)} =e^{-\beta\Delta E}

所以：

\boxed{ A_{i\to j} = \min(1,e^{-\beta\Delta E}) }

如果：

\Delta E=E_j-E_i\le 0

说明新状态能量更低。

那么：

e^{-\beta\Delta E}\ge 1

所以：

A=1

总是接受。

如果：

\Delta E>0

说明能量升高。

那么：

A=e^{-\beta\Delta E}<1

以一定概率接受。

温度越高，$\beta=1/k_BT$ 越小，接受升能态概率越大。

算法：Metropolis 采样

1. 从当前构型 i，试探一个新构型 j（对称提议 T(i→j) = T(j→i)）
2. 计算能量差 ΔE = Ej - Ei
3. 接受概率：
   A(i→j) = min(1, e^{-β·ΔE})

   即：
   - 若 ΔE ≤ 0（能量降低）→ 总是接受
   - 若 ΔE > 0（能量升高）→ 以概率 e^{-β·ΔE} 接受

4. 生成 [0,1] 均匀随机数 r：
   - 若 r ≤ e^{-β·ΔE} → 接受
   - 否则 → 拒绝，保留旧构型

5. 回到步骤 1

取 min：是因为我们要让算法跑得最快（接受率最大化），同时必须遵守详细平衡。

接受升能态：是为了让链能翻越势垒，不困在局部最优，从而获得正确的物理分布。

消掉 Z：是让这个算法在超大规模组合空间中变得可行（这是它横扫各个学科的根本原因）。

🔢 数值例子

设 $\beta=1$ ，当前能量 $E_i = 2$ ，试探态 $E_j = 3$ ：

$\Delta E = 3 - 2 = +1 > 0$
接受概率 $A = e^{-1} \approx 0.368$ 这个是接受去到新的试探态的概率上界因为比这个小的代表能量差值会更大
生成随机数 $r = 0.5$ → 拒绝（$0.5 > 0.368$）
生成随机数 $r = 0.2$ → 接受（$0.2 \leq 0.368$）

7.3 Metropolis-Hastings 算法

推广到非对称提议分布：

A(i\to j) = \min\!\left(1,\; \frac{p_j \cdot T(i|j)}{p_i \cdot T(j|i)}\right)

$T(j|i)$ ：提议从 $i$ 产生 $j$ 的概率
当 $T$ 对称时，退化为标准 Metropolis

7.4 MCMC 积分

一次完整的 MCMC 期望值估计大概是这样展开的：先随便选一个初始构型 $x_0$ ，用 Metropolis（或其他满足细致平衡的规则）反复更新,产生 $x_1,x_2,\dots$ ；观察能量等量随步数的曲线 , 找到大致的热化拐点，将拐点之前的 $n_{\text{equil}}$ 步全部丢弃；对剩下的样本 $x_{n_{\text{equil}}+1},\dots,x_N$ ，计算感兴趣物理量的自相关函数 $C(t)$ ，从中拟合出 $\tau_{\text{int}}$ ；用剩余样本算出 $\langle A\rangle\approx\frac1N\sum A(x_i)$ 作为最佳估计,并用 $\sigma_{\bar A}\approx\sqrt{2\tau_{\text{int}}\operatorname{Var}(A)/N}$ 给出这个估计值的误差范围；如果发现 $\tau_{\text{int}}$ 太大（比如在相变点附近，系统会出现"临界慢化"，关联时间急剧变长），说明要么增加总步数 $N$ ，要么换用效率更高的算法（比如针对伊辛模型专门设计的 cluster algorithm，它能大幅缩短关联时间）。

MCMC 估计期望值

如果马尔可夫链收敛到目标分布 $p(x)$ ，那么样本： $$ x1,x2,\dots,xN $$ 虽然不是独立的，但长时间平均仍然可以估计： $$ \langle A\ranglep=\int A(x)p(x)dx $$ 即： $$ \boxed{ \langle A\rangle \approx \frac1N\sum{i=1}^N A(xi) } $$

MCMC 的两个重要问题

1. Burn-in 平衡段

一开始链还没到平衡分布，需要扔掉前面一段： $$ n_{\text{equil}} $$ 这叫 thermalization 或 burn-in。

2. 样本相关性

MCMC 样本不是独立的。

如果当前状态是 $x_i$ ，下一个状态 $x_{i+1}$ 通常和它很像。

所以有效样本数小于总样本数。

误差估计应考虑自相关时间： $$ N{\text{eff}}\approx \frac{N}{2\tau{\text{int}}} $$ 误差大约： $$ \sigma{\bar A} \approx \sqrt{ \frac{\operatorname{Var}(A)}{N{\text{eff}}} } $$

# 伪代码：MCMC 积分框架

def mcmc_integral(target_pdf, proposal, n_equil, n_samples):
    """
    target_pdf: 目标分布 p(x)
    proposal: 提议函数，给定 x 产生 x'
    n_equil: 平衡步数（扔掉）
    n_samples: 有效采样步数
    """
    x = initial_state()
    samples = []

    for t in range(n_equil + n_samples):
        x_new = proposal(x)
        # Metropolis 接受/拒绝
        A = min(1, target_pdf(x_new) / target_pdf(x) * proposal_ratio)
        if np.random.uniform() < A:
            x = x_new
        if t >= n_equil:
            samples.append(x)

    return np.mean(samples), np.std(samples) / np.sqrt(n_samples)

📊 三类采样方法对比

特性	直接采样	传统MC积分	MCMC
采样点是否独立	✅ 独立	✅ 独立	❌ 关联
能处理的分布	简单分布	均匀/简单加权	几乎任意分布
关键需求	已知采样算法	计算 $f(x)$	知道分布比例
收敛速度	快	中等	慢（需平衡+关联）

💡 物理直觉

MCMC 的核心矛盾：牺牲"样本独立"换取"能处理任意分布"。就像你没法直接生成复杂波函数，但可以让体系按 Boltzmann 概率自己"游走"到各个态。

8. 二维 Ising 模型

8.1 模型定义

    ↑       ↑
    ↓   ↑   ↓   ← 自旋 σ_i = ±1
    ↑   ↓   ↑
        ↑

L × L 格点，总自旋数 N = L²

能量（最近邻相互作用）：

E = -J \sum_{\langle i,j\rangle} \sigma_i \sigma_j

$J>0$ ：铁磁耦合（邻居同向 → 能量低）同向，贡献 -J ；反向，贡献 +J 。因此同向更低，因此低温下系统倾向于全部向上或全部向下
构型总数：$2^N$（指数级！必须 MC）总自旋数为 L的平方，其中每个自旋有两种取值，所以总构型数为 $2^{N}$ ，如果 L为20，那么最后构型数为 $2^{400}$ 天文数字，必须用 Monte Carlo采样！！！

观测量：

物理量	公式
磁化强度	$M = \frac{1}{N}\sum_i \sigma_i$
能量	$E = -J\sum_{\langle i,j\rangle} \sigma_i \sigma_j$
比热	$C_V = \frac{1}{k_B T^2}(\langle E^2\rangle - \langle E\rangle^2)$
磁化率	$\chi = \frac{1}{k_B T}(\langle M^2\rangle - \langle M\rangle^2)$

8.2 Metropolis Ising 算法

一次 sweep 通常指尝试更新 $N=L^2$ 次。

算法：

1、初始化自旋：

\sigma_{ij}=\pm 1

2、随机选择一个格点 $(i,j)$ 。

3、计算它四个邻居和：

S=\sigma_{i+1,j}+\sigma_{i-1,j}+\sigma_{i,j+1}+\sigma_{i,j-1}

4、计算：

\Delta E=2J\sigma_{ij}

5、接受规则：

若 $\Delta E\le 0$ ，接受；
若 $\Delta E>0$ ，以概率 $e^{-\beta\Delta E}$ 接受。

6、重复。

为了减少边界效应，常用周期边界： $$ \sigma{i+L,j}=\sigma{i,j} $$ 也就是格子左右相连，上下相连，拓扑上像一个环面。

代码中常用模运算：

Python

up = spins[(i+1) % L, j]

down = spins[(i-1) % L, j]

right = spins[i, (j+1) % L]

left = spins[i, (j-1) % L]

Metropolis MC for 2D Ising Model:

Parameters: L (格点大小), T (温度), β = 1/(k_B T), J (耦合常数)

1. 随机初始化 L×L 格点自旋 σ[i][j] = ±1

2. FOR mc_cycle = 1 to n_cycles:
     FOR 每个格点 (随机顺序遍历一次 = 1 sweep):
       a. 翻转该自旋，计算能量变化 ΔE
       b. 2D Ising 中 ΔE 只有 5 种可能：
          ΔE ∈ { -8J, -4J, 0, +4J, +8J }
       c. 若 ΔE ≤ 0 → 接受翻转
          若 ΔE > 0 → 以 p = e^{-β·ΔE} 接受
     ENDFOR
     记录 E, M 等观测量

3. 输出 ⟨E⟩, ⟨M⟩, C_V, χ 及其统计误差

💡 为什么 ΔE 只有 5 种值？

Metropolis 算法中，每次尝试翻转一个自旋：

\sigma_k\to -\sigma_k

只有与这个自旋相关的能量项会变化。

原来和邻居的相互作用能：

E_{\text{old}} = -J\sigma_k\sum_{\text{nn}}\sigma_j

翻转后：

\sigma_k'=-\sigma_k

所以：

E_{\text{new}} = -J(-\sigma_k)\sum_{\text{nn}}\sigma_j

能量变化：

\Delta E = E_{\text{new}}-E_{\text{old}}

所以：

\boxed{ \Delta E = 2J\sigma_k \sum_{\text{邻居}}\sigma_j }

二维正方格子每个格点有 4 个邻居。

邻居和可能是：

-4,-2,0,2,4

因此：

\Delta E = 2J\sigma_k (-4,-2,0,2,4)

所以只可能是：

\boxed{ \Delta E\in\{-8J,-4J,0,4J,8J\} }

2D 正方格子中每个自旋有 4 个邻居（上下左右）。翻转一个 $\sigma_k \to -\sigma_k$ ：
$\Delta E = 2J \sigma_k \sum_{\text{邻居}} \sigma_{\text{邻居}}$
4 个邻居的和只能是 $-4, -2, 0, 2, 4$ → $\Delta E = -8J, -4J, 0, 4J, 8J$ （5 种）。

🔢 典型结果

二维 Ising 模型有精确临界温度：

\boxed{ k_BT_c = \frac{2J}{\ln(1+\sqrt2)} }

数值：

\boxed{ T_c\approx 2.269\frac{J}{k_B} }

现象：

$T：铁磁相，自发磁化；
$T>T_c$ ：顺磁相，平均磁化为 0；
$T=T_c$ ：涨落最大，比热和磁化率出现临界行为。

物理现象	说明
一维 Ising	无相变（$T_c=0$）
二维 Ising	二级相变，$Tc \approx 2.269 J/kB$
低于 $T_c$	自发磁化 $M \neq 0$
高于 $T_c$	$M=0$ ，顺磁态

8.3 KT 相变（2016 诺贝尔物理奖）

二维 XY 模型自旋不是 $\pm1$ ，而是角度：

\mathbf S_i=(\cos\theta_i,\sin\theta_i)

低温下会出现涡旋-反涡旋束缚对。

高温时涡旋对解离。

这就是 Kosterlitz-Thouless 相变。

特点：

没有普通意义上的长程磁化；
是拓扑缺陷驱动的相变；
2016 年诺贝尔物理奖相关。

二维 XY 模型中的 Kosterlitz-Thouless 相变：低温下涡旋-反涡旋成对束缚 → 高温下解耦 → 拓扑相变（无自发磁化，但有超流/超导性质的突变）。

9. 量子蒙特卡洛（QMC）

💡 为什么需要 QMC？

主要用于量子多体问题！！！

量子多体波函数 $\Psi(R_1, \ldots, R_N)$ 不是简单概率分布，维度是3N。如果N很大，直接网格法解析求解几乎不可能，所以用Monte Carlo在高维构型空间采样。QMC 是唯一能精确计算大体系的方法。

9.1 变分量子蒙特卡洛（VMC）

能量期望：

E_T = \frac{ \int \Psi_T^*(R)\hat H\Psi_T(R)\,dR } { \int |\Psi_T(R)|^2dR }

假设波函数实数，写成：

E_T = \frac{ \int |\Psi_T(R)|^2 \frac{\hat H\Psi_T(R)}{\Psi_T(R)} dR } { \int |\Psi_T(R)|^2dR }

定义概率分布：

P(R)= \frac{|\Psi_T(R)|^2} { \int |\Psi_T(R)|^2dR }

定义局域能量：

\boxed{ E_L(R)= \frac{\hat H\Psi_T(R)}{\Psi_T(R)} }

于是：

\boxed{ E_T= \int P(R)E_L(R)dR }

用 MC 估计：

\boxed{ E_T\approx \frac1N\sum_{i=1}^N E_L(R_i) }

其中：

R_i\sim |\Psi_T(R)|^2

框架

构造试探波函数 $\Psi_T(R; \alpha)$ ，含变分参数 $\alpha$
用 Metropolis 从 $P(R)\propto |\Psi_T(R;\alpha)|^2$ 采样空间构型 $R$
对每个样本计算局域能量：

E_L(R) = \frac{\hat{H}\Psi_T(R)}{\Psi_T(R)}

求能量期望值：

\langle E_{L} \rangle = \frac{\int |\Psi_T|^2 E_L dR}{\int |\Psi_T|^2 dR} \approx \frac{1}{N}\sum_{i=1}^{N} E_L(R_i) = E(\alpha)

优化 $\alpha$ 使 $E(\alpha)$ 最小

💡 说人话

VMC = 聪明的猜谜：猜一个波函数形状，MC 采样→算能量→调参数让能量尽可能低。如果 $E_L$ 不随构型变化，说明波函数就是精确的！

为什么局域能量恒定说明是精确本征态？

如果：

\Psi_T

正好是 Hamiltonian 的本征态：

\hat H\Psi_T=E\Psi_T

那么：

E_L(R)= \frac{\hat H\Psi_T(R)}{\Psi_T(R)} = \frac{E\Psi_T(R)}{\Psi_T(R)} = E

它不依赖 $R$ 。

所以局域能量方差为零。

这就是 VMC 中常说的：

精确波函数的局域能量没有涨落。

🔢 谐振子 VMC 例子

试探波函数：$\Psi_T(x) = e^{-\alpha x^2/2}$

哈密顿量：$\hat{H} = -\frac{1}{2}\frac{d^2}{dx^2} + \frac{1}{2}x^2$

局域能量： $$ EL(x) = \frac{\hat{H}\PsiT}{\Psi_T} = \frac{\alpha}{2} + \frac{1-\alpha^{2}}{2} x^{2} $$

$\alpha=1$ 时：$E_L = 1/2$（与 $x$ 无关！）= 精确解标准一维谐振子基态能量！
对 $\alpha=1$ 做 VMC：$\langle E\rangle \to 1$，方差 $\to 0$

import numpy as np

def vmc_harmonic_oscillator(alpha, N_mc=100000):
    """VMC for 1D harmonic oscillator"""
    # Metropolis 从 |psi|^2 ∝ exp(-alpha*x^2) 采样
    x = 0.0
    delta = 1.0 / np.sqrt(alpha)
    samples = []

    for i in range(N_mc):
        x_new = x + delta * np.random.uniform(-1, 1)
        # 接受率 = |psi(x_new)/psi(x)|^2
        ratio = np.exp(-alpha * (x_new**2 - x**2))
        if np.random.uniform() < min(1, ratio):
            x = x_new
        samples.append(x)

    # 局域能量
    x_arr = np.array(samples[N_mc//2:])  # 扔掉平衡段
    e_local = alpha**2 + x_arr**2 * (1 - alpha**4)

    return np.mean(e_local), np.std(e_local)/np.sqrt(len(x_arr))

E, err = vmc_harmonic_oscillator(alpha=1.0)
print(f"<E> = {E:.4f} ± {err:.4f}")  # 应接近 1.0

9.2 扩散蒙特卡洛（DMC）

DMC 的思想来自虚时间 Schrödinger 方程。

真实时间 Schrödinger 方程：

i\frac{\partial \Psi}{\partial t} = \hat H\Psi

做虚时间变换：

t=-i\tau

则方程变成类似：

-\frac{\partial \Psi}{\partial \tau} = (\hat H-E_T)\Psi

或者：

\frac{\partial \Psi}{\partial \tau} = -(\hat H-E_T)\Psi

设：

\hat H= -\frac12\nabla^2+V(R)

则：

\frac{\partial \Psi}{\partial \tau} = \frac12\nabla^2\Psi - [V(R)-E_T]\Psi

这有两部分：

扩散项：
$\frac12\nabla^2\Psi$
生灭项：
$-[V(R)-E_T]\Psi$

为什么虚时间会投影到基态？

把初始波函数按本征态展开：

\Psi(R,0)=\sum_n c_n\phi_n(R)

其中：

\hat H\phi_n=E_n\phi_n

虚时间演化：

\Psi(R,\tau) = e^{-(\hat H-E_T)\tau}\Psi(R,0)

如果取 $E_T$ 接近基态能量，那么高能态：

E_n>E_0

会衰减得更快。

当：

\tau\to\infty

只剩基态：

\Psi(R,\tau)\propto \phi_0(R)

所以 DMC 可以投影出基态。

DMC 的 walker 图像

用一群 walkers 表示波函数分布。

每个 walker 做两件事：

1. 扩散

R\to R+\sqrt{\Delta\tau}\,\eta

其中 $\eta$ 是高斯随机数。

2. 分支 branching

根据势能决定 walker 生灭。

如果：

V(R)<E_T

该区域权重大，walker 增殖。

如果：

V(R)>E_T

该区域权重小，walker 死亡。

最后 walker 分布趋向基态波函数相关分布。

DMC的难点：费米子符号问题

玻色子基态可以取正，DMC 比较自然。

但费米子波函数反对称，会有正负号：

\Psi(...,r_i,...,r_j,...) = -\Psi(...,r_j,...,r_i,...)

概率解释困难，正负 walker 会相互抵消，噪声指数增长。

这就是 sign problem。

常用近似是 fixed-node approximation：用试探波函数的节点面固定符号结构。

对比	VMC	DMC
原理	变分原理 + Metropolis 采样	虚时间 Schrödinger 方程演化
波函数	需要好的试探函数	可以精确得到基态（给定节点）
实现	随机行走在坐标空间	walkers 的扩散 + 分支过程
难点	试探波函数选择	费米子符号问题（Sign Problem）

💡 DMC 物理直觉

把 Schrödinger 方程的 $t \to -i\tau$ （虚时间），它就变成了扩散方程 + 势能项。Walkers（采样粒子）像布朗粒子一样扩散，同时在高势能区域"死亡"、低势能区域"增殖"。最终 walker 分布 → 基态波函数。

9.3 路径积分蒙特卡洛（PIMC）

PIMC 用于有限温度量子系统。

目标是计算配分函数：

Z=\operatorname{Tr}(e^{-\beta \hat H})

把：

\beta

分成 $M$ 个小片：

\Delta\tau=\frac{\beta}{M}

则：

e^{-\beta H} = \left(e^{-\Delta\tau H}\right)^M

插入 $M$ 次完备坐标基：

I=\int dR |R\rangle\langle R|

得到路径积分形式。

物理图像：

一个量子粒子在有限温度下等价于一条闭合的虚时间路径。多粒子系统变成很多“聚合物环”的统计采样。

PIMC 适合：

有限温度；
玻色体系；
超流；
量子液体；
量子晶体。

但费米子仍有严重 sign problem。

对多个虚时间路径采样，可直接计算配分函数 $Z = \text{Tr}(e^{-\beta H})$ ，适用于有限温度量子体系。

10. 总结对比表

📊 蒙特卡洛方法族谱

方法	核心思想	适用场景	关键技巧
MC 积分	均匀/加权随机采样	高维积分、算 π	投点法、期望值法
重要性采样	在函数大值区多采样	$f(x)$ 有尖峰	选择好的 $F(x)\approx f(x)$
MCMC	马尔可夫链趋于平衡分布	任意复杂分布	Metropolis-Hastings
Metropolis	能量差决定接受概率	统计物理、Boltzmann 分布	$A=\min(1,e^{-\beta\Delta E})$
VMC	$	\Psi_T	^2$ 采样 + 变分优化	量子多体基态	试探波函数是关键
DMC	虚时间扩散 + 分支	精确基态能量	克服符号问题
PIMC	虚时间路径积分采样	有限温度量子体系	配分函数计算

🧠 核心要点速记

概念	一句话
中心极限定理	$m$ 个样本均值 ~ 正态分布，方差为 $\sigma^2/m$
MC 积分的优势	误差与维度无关（高维积分的王牌）
细致平衡	$W_{ij} p_i = W_{ji} p_j$ — 平衡态的充分条件
Metropolis 接受率	$A = \min(1, e^{-\beta\Delta E})$ — 无比简洁
MCMC 核心矛盾	牺牲独立性 ↔ 换取任意分布的采样能力
VMC 变分原理	$\langle E\rangle \geq E_0$ （永远不低于真实基态能）
DMC 符号问题	费米子交换 → 负概率 → 计算困难

Akuiro 整理补充 2026-07-05

蒙特卡洛方法初步

蒙特卡洛方法初步

目录

1. 蒙特卡洛方法概述

📖 核心思想

🕰 历史起源

💡 物理直觉

🔄 方法流程

2. 概率统计基础

2.1 概率密度函数（PDF）

💡 数值例子

2.2 中心极限定理 ⭐

💡 说人话

3. 随机数生成

3.1 伪随机数要求

3.2 高斯分布采样：Box-Muller 法

💡 直观验证

4. 蒙特卡洛积分

4.1 空间随机投点法（打靶法）

4.2 期望值法 ⭐

方差估计

🔢 MC 算 π 例子

4.3 多维积分

💡 说人话

5. 重要性采样（Importance Sampling）

5.1 动机

5.2 原理

最优重要性分布

💡 说人话

🔢 数值示例

方法一：均匀采样

方法二：重要性采样 选取 p(x)∝e−xp(x) \propto e^{-x}p(x)∝e−x

如何从 p(x)∝e−xp(x)\propto e^{-x}p(x)∝e−x 抽样？

📊 对比

6. 随机行走与马尔可夫链

6.1 一维随机行走

💡 物理直觉

6.2 马尔可夫过程

转移概率矩阵

时间演化（矩阵形式）

🔢 数值例子

💡 自然结论

7. Metropolis 算法与 MCMC

7.1 平衡态与细致平衡

💡 说人话

7.2 Metropolis 算法 ⭐（20世纪十大算法之一）

🔢 数值例子

7.3 Metropolis-Hastings 算法

7.4 MCMC 积分

MCMC 估计期望值

MCMC 的两个重要问题

1. Burn-in 平衡段

2. 样本相关性

📊 三类采样方法对比

💡 物理直觉

8. 二维 Ising 模型

8.1 模型定义

8.2 Metropolis Ising 算法

💡 为什么 ΔE 只有 5 种值？

🔢 典型结果

8.3 KT 相变（2016 诺贝尔物理奖）

9. 量子蒙特卡洛（QMC）

💡 为什么需要 QMC？

9.1 变分量子蒙特卡洛（VMC）

框架

💡 说人话

为什么局域能量恒定说明是精确本征态？

🔢 谐振子 VMC 例子

9.2 扩散蒙特卡洛（DMC）

为什么虚时间会投影到基态？

DMC 的 walker 图像

1. 扩散

2. 分支 branching

DMC的难点：费米子符号问题

💡 DMC 物理直觉

9.3 路径积分蒙特卡洛（PIMC）

10. 总结对比表

📊 蒙特卡洛方法族谱

🧠 核心要点速记

方法二：重要性采样选取 $p(x) \propto e^{-x}$

如何从 $p(x)\propto e^{-x}$ 抽样？