蒙特卡洛方法

2026-07-22

目录

1. 蒙特卡洛方法概述
2. 概率统计基础
3. 随机数生成
4. 蒙特卡洛积分
5. 重要性采样
6. 随机行走与马尔可夫链
7. Metropolis 算法与 MCMC
8. 二维 Ising 模型
9. 量子蒙特卡洛
10. 总结对比表

1. 蒙特卡洛方法概述

📖 核心思想

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

🕰 历史起源

💡 物理直觉

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

🔄 方法流程

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

2. 概率统计基础

2.1 概率密度函数（PDF）

💡 数值例子

某物理量 $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 伪随机数要求

💡 实践建议：使用系统自带随机数生成器（如 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$。

💡 说人话

算 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$

$$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]的 均匀随机分布！

📊 对比

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$ 步后位置：$X_n = \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} w_1(t)\\ w_2(t)\\ w_3(t)\\ w_4(t) \end{pmatrix} $$ 其中： $$ w_i(t)=P(\text{系统在状态 }i) $$ 转移矩阵 $W$ 定义为： $$ W_{ij}=P(j\to i) $$ 注意你的讲义采用的是：

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

因此每列之和为 1： $$ \sum_i 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}\, p_i = W_{j\to i}\, p_j $$

意思是：$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)$，那么样本： $$ x_1,x_2,\dots,x_N $$ 虽然不是独立的，但长时间平均仍然可以估计： $$ \langle A\rangle_p=\int A(x)p(x)dx $$ 即： $$ \boxed{ \langle A\rangle \approx \frac1N\sum_{i=1}^N A(x_i) } $$

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)

📊 三类采样方法对比

💡 物理直觉

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采样！！！

观测量：

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$：涨落最大，比热和磁化率出现临界行为。

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 $$

框架

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

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

1. 求能量期望值：

$$ \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) $$

1. 优化 $\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$

局域能量： $$ E_L(x) = \frac{\hat{H}\Psi_T}{\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 $$ 这有两部分：

1. 扩散项：

$$ \frac12\nabla^2\Psi $$

1. 生灭项：

$$ -[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)

如果：

$$ V(R)>E_T $$ 该区域权重小，walker 死亡。

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

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

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

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

$$ \Psi(...,r_i,...,r_j,...) = -\Psi(...,r_j,...,r_i,...) $$ 概率解释困难，正负 walker 会相互抵消，噪声指数增长。

这就是 sign problem。

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

💡 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. 总结对比表

📊 蒙特卡洛方法族谱

🧠 核心要点速记

Akuiro 整理补充 2026-07-05