1. 两个路径积分的例子1.1. Ito 积分形式的路径积分例子1.2. 关于随机热力学的例子2. 一般形式的路径积分一些前置概念请见随机分析、Langevin 方程与 Fokker-Planck 方程。1. 两个路径积分的例子1.1. Ito 积分形式的路径积分例子考虑过阻尼 Langevin 方程$$ \dot{x}(t) = f(x(t)) + g(x(t))\xi(t), $$其中 $\xi(t)$ 为白噪声，满足 $\langle \xi(t) \rangle = 0$，$\langle \xi(t)\xi(t') \rangle = 2D\delta(t-t')$.考虑将时间做离散化处理，将时间区间 $[0, \tau]$ 等分为 $K$ 份，即 $\tau = K\Delta t$，简记 $t^k \equiv k\Delta t$，$x^k \equiv x(t^k)$. 将离散化后的 Langevin 方程记为$$ x^{k+1} - x^{k} \overset{(\alpha)}{=} f(\bar{x}^k)\Delta t + g(\bar{x}^k)\Delta \omega^k, $$其中 $\alpha$ 涉及到随机积分的表达形式（稍后会提到）， $\Delta \omega^k \equiv \omega^{k+1} - \omega^k = \omega(t^{k+1}) - \omega(t^k)$， $\omega(t)$ 为强度为 $2D$ 的维纳过程。维纳过程由高斯独立增量给出，因此 $\Delta\omega^k$ 满足如下性质$$ \langle \Delta\omega^k \rangle = 0, \quad \langle \Delta\omega^k \Delta\omega^{k'} \rangle = 2D \Delta t \delta_{kk'}. $$如何处理函数 $f$ 和 $g$ 的离散化形式涉及到选择随机积分表达形式的问题。定义 $\bar{x}^k = \alpha x^{k+1} + (1-\alpha) x^k$，则常用的 Ito 和 Stratonovich 两种取法分别对应 $\alpha = 0$ 和 $\alpha = 1/2$.简单起见，我们首先讨论 Ito 形式的随机积分，即 $\alpha = 0$ 的情况。此时离散化的 Langevin 方程为$$ x^{k+1} - x^{k} \overset{(\alpha = 0)}{=} f(x^k, t^k)\Delta t + g(x^k, t^k)\Delta \omega^k. $$对于给定的初始位置 $x^0$ ，当我们指定了维纳过程的具体实现时，即 $\omega(t)$ 在每个离散时刻都被确定下来 $W_\tau \equiv \{ \Delta\omega^0, \Delta\omega^1, \cdots, \Delta\omega^{K-1} \}$，随机轨迹也成为一个确定性的轨迹 $X_\tau \equiv \{x^1, x^2, \cdots, x^K \}$. 生成这样的维纳过程的概率为$$ \begin{aligned} \mathcal{P}[W_\tau] &= \prod_{k=0}^{K-1} P(\Delta \omega^k) \\ &= \prod_{k=0}^{K-1} \frac{1}{\sqrt{4\pi D \Delta t}} \exp\left[ -\frac{(\Delta\omega^k)^2}{4D\Delta t} \right]， \end{aligned} $$这是因为 $\{ \Delta\omega^0, \Delta\omega^1, \cdots, \Delta\omega^{K-1} \}$ 是独立同分布的高斯增量。由离散化的 Langevin 方程可得$$ \Delta\omega^k \overset{(\alpha = 0)}{=} \frac{1}{g(x^k, t^k)} \left( x^{k+1} - x^{k} - f(x^k, t^k)\Delta t \right), $$带入 $\mathcal{P}[W_\tau]$ 可得$$ \mathcal{P}[W_\tau] = \prod_{k=0}^{K-1} \frac{1}{\sqrt{4\pi D \Delta t}} \exp\left[ -\frac{\Delta t}{4D g^2(x^k, t^k)} \left( \frac{x^{k+1} - x^{k}}{\Delta t} - f(x^k, t^k) \right)^2 \right]. $$现在距离路径积分的表达式只剩最后一步，即 $\mathcal{P}[W_\tau]$ 到 $\mathcal{P}[X_\tau | x^0]$ 的变换，它们之间相差了一个 Jacobi 行列式。由离散化的 Langevin 方程可得 Jacobi 行列式为$$ |J| \equiv \left| \frac{\partial(x^1, x^2, \cdots, x^K)}{\partial(\Delta\omega^0, \Delta\omega^1, \cdots, \Delta\omega^{K-1})} \right| = \prod_{k=0}^{K-1} \left| g(x^k, t^k) \right| = \prod_{k=0}^{K-1} \sqrt{g^2(x^k, t^k)}. $$由此可得 Langevin 方程生成的轨迹概率为$$ \begin{aligned} \mathcal{P}[X_\tau] ={}& p(x^0) \cdot \mathcal{P}[W_\tau] \cdot \left| J^{-1} \right| \\ ={}& p(x^0) \cdot \prod_{k=0}^{K-1} \frac{1}{\sqrt{4\pi D g^2(x^k, t^k) \Delta t}} \\ &\times \exp\left[ -\frac{\Delta t}{4D g^2(x^k, t^k)} \left( \frac{x^{k+1} - x^{k}}{\Delta t} - f(x^k, t^k) \right)^2 \right] \\ ={}& p(x^0) \cdot \left( \prod_{k=0}^{K-1} \frac{1}{\sqrt{4\pi D g^2(x^k, t^k) \Delta t}} \right) \\ &\times \exp\left[ -\frac{1}{4D} \sum_{k=0}^{K-1} \Delta t \left\{ \frac{1}{g^2(x^k, t^k)} \left( \frac{x^{k+1} - x^{k}}{\Delta t} - f(x^k, t^k) \right)^2 \right\} \right] \\ ={}& \underbrace{\mathcal{N}[X_\tau]}_{\text{normalizing factor}} \cdot \exp\{ -\underbrace{\mathcal{S}[X_\tau]}_{\text{path action}} \}, \end{aligned} $$其中归一化因子为$$ \mathcal{N}[X_\tau] = p(x^0) \cdot \prod_{k=0}^{K-1} \frac{1}{\sqrt{4\pi D g^2(x^k, t^k) \Delta t}}, $$轨迹作用量为$$ \mathcal{S}[X_\tau] = \frac{1}{4D} \sum_{k=0}^{K-1} \Delta t \left\{ \frac{1}{g^2(x^k, t^k)} \left( \frac{x^{k+1} - x^{k}}{\Delta t} - f(x^k, t^k) \right)^2 \right\}, $$路径积分的轨迹测度为$$ \int \mathcal{D}[X_\tau] \equiv \lim_{\substack{\Delta t \to 0, ~ K \to \infty \\K\Delta t = \tau \text{ fixed}}} \prod_{k=0}^{K} \int_{-\infty}^{+\infty} \mathrm{d}x^k. $$归一化条件为$$ \int \mathcal{D}[X_\tau] \mathcal{P}[X_\tau] = \int \mathcal{D}[X_\tau] p(x^0) \mathcal{P}[X_\tau|x^0] = 1. $$在连续时间极限下，即 $\Delta t \to 0, K \to \infty, K\Delta t = \tau$ 时，作用量可以写成时间积分的形式：$$ \mathcal{S}[X_\tau] = \frac{1}{4D}\int_0^\tau \left[ \frac{\dot{x} - f(x, t)}{g(x, t)} \right]^2 \mathrm{d}t. $$这看起来像是 $\xi^2(t)$ 的时间积分，因为 Langevin 方程可以重排为 $\xi(t) = \frac{\dot{x} - f(x, t)}{g(x, t)}$，但下一个例子会说明这个想法并不总是正确。1.2. 关于随机热力学的例子在随机热力学中，通常考虑过阻尼 Langevin 方程$$ \dot{x}(t) = \mu F(x(t), \lambda(t)) + \xi(t), $$其中 $\lambda(t)$ 为外部控制参数， $\xi(t)$ 为白噪声，满足 $\langle \xi(t) \rangle = 0$，$\langle \xi(t)\xi(t') \rangle = 2D\delta(t-t')$.现在做与上一节相同的离散化处理。唯一的不同点在于，在随机热力学中，人们通常取 Stratonovich 形式的随机积分，因为这样的形式更符合实际物理情形。因此，离散化的 Langevin 方程为$$ \begin{aligned} x^{k+1} - x^{k} &= \mu F(\bar{x}^k)\Delta t + \Delta \omega^k \\ &\approx \frac{\mu}{2} \left[ F(x^{k+1}, \lambda^{k+1}) + F(x^{k}, \lambda^{k}) \right] \Delta t + \Delta\omega^k. \end{aligned} $$此时唯一发生变化的是 Jacobi 行列式。Jacobi 矩阵为$$ \begin{aligned} J^{-1} &= \frac{\partial(\Delta\omega^0, \Delta\omega^1, \cdots, \Delta\omega^{K-1})}{\partial(x^1, x^2, \cdots, x^K)} \\ &= \begin{pmatrix} 1 - \frac{\mu\Delta t}{2}\partial_{x^1} F(x^1, \lambda^1) & 0 & \cdots & 0 \\ -1 - \frac{\mu\Delta t}{2}\partial_{x^1} F(x^1, \lambda^1) & 1 - \frac{\mu\Delta t}{2}\partial_{x^2}F(x^2, \lambda^2) & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & 1 - \frac{\mu\Delta t}{2}\partial_{x^K} F(x^K, \lambda^K) \end{pmatrix}. \end{aligned} $$因此 Jacobi 行列式为$$ \begin{aligned} \left| J^{-1} \right| &= \prod_{k=0}^{K-1} \left[ 1 - \frac{\mu\Delta t}{2} \frac{\partial F(x^{k+1}, \lambda^{k+1})}{\partial x^{k+1}} \right] \\ &= \exp\left[ \sum_{k=0}^{K-1} \ln\left( 1 - \frac{\mu\Delta t}{2} \frac{\partial F(x^{k+1}, \lambda^{k+1})}{\partial x^{k+1}} \right) \right] \\ &\approx \exp\left[ -\sum_{k=0}^{K-1} \frac{\mu\Delta t}{2} \frac{\partial F(x^{k+1}, \lambda^{k+1})}{\partial x^{k+1}} \right]. \end{aligned} $$带入可得轨迹概率表达式$$ \begin{aligned} \mathcal{P}[X_\tau] ={}& p(x^0) \cdot \mathcal{P}[W_\tau] \cdot \left| J^{-1} \right| \\ ={}& \frac{p(x^0)}{(4\pi D \Delta t)^{K/2}} \cdot \exp\left[ -\frac{\Delta t}{4D} \sum_{k=0}^{K-1} \left( \frac{x^{k+1} - x^k}{\Delta t} - \mu F(\bar{x}^k) \right)^2 \right. \\ & \left. - \frac{\mu\Delta t}{2} \sum_{k=0}^{K-1} \frac{\partial F(x^{k+1}, \lambda^{k+1})}{\partial x^{k+1}} \right] \\ ={}& \underbrace{\mathcal{N}[X_\tau]}_{\text{normalizing factor}} \cdot \exp\{ -\underbrace{\mathcal{S}[X_\tau]}_{\text{path action}} \}, \end{aligned} $$其中归一化因子为$$ \mathcal{N}[X_\tau] = \frac{p(x^0)}{(4\pi D \Delta t)^{K/2}}, $$作用量为$$ \mathcal{S}[X_\tau] = \frac{\Delta t}{4D} \sum_{k=0}^{K-1} \left( \frac{x^{k+1} - x^k}{\Delta t} - \mu F(\bar{x}^k) \right)^2 + \frac{\mu\Delta t}{2} \sum_{k=0}^{K-1} \frac{\partial F(x^{k+1}, \lambda^{k+1})}{\partial x^{k+1}}. $$在连续时间极限下，即 $\Delta t \to 0, K \to \infty, K\Delta t = \tau$ 时，作用量 $\mathcal{S}$ 同样可以写成时间积分的形式：$$ \mathcal{S}[X_\tau] = \underbrace{ \frac{1}{4D}\int_0^\tau \left( \dot{x} - \mu F(x, t) \right)^2 \mathrm{d}t }_{\text{Gaussian noise weight}} + \underbrace{ \frac{\mu}{2}\int_0^\tau \frac{\partial F(x, t)}{\partial x} \mathrm{d}t }_{\text{Jacobian weight}}. $$可以发现轨迹作用量由两部分组成。一部分是高斯白噪声的贡献，另一部分是 Jacobi 行列式的贡献。实际上第一个例子中的 Jacobi 行列式同样对轨迹概率有贡献，只是被整合到了归一化因子中。白噪声的贡献在连续时间极限下总是相同的，与取什么样的随机积分形式无关，而 Jacobi 行列式的贡献与选取的形式有关，不同的积分形式会给出不同的项。2. 一般形式的路径积分考虑一般的过阻尼 Langevin 方程（也称为线性 Langevin 方程，因为最高只出现了 $x$ 的一阶时间导数）$$ \dot{x} = f(x(t)) + g(x(t)) \xi(t), $$其中 $\xi(t)$ 为高斯白噪声，满足 $\langle \xi(t) \rangle = 0$，$\langle \xi(t)\xi(t') \rangle = 2D\delta(t-t')$. 高斯白噪声的轨迹概率 $\mathcal{P}[W_\tau]$ 总是相同的，因为它是白噪声自身的性质，与 $f(x(t))$ 和 $g(x(t))$ 的形式无关：$$ \mathcal{P}[W_\tau] = \prod_{k=0}^{K-1} \frac{1}{\sqrt{4\pi D \Delta t}} \exp\left[ -\frac{(\Delta\omega^k)^2}{4D\Delta t} \right]. $$将 Langevin 方程做离散化处理得到$$ x^{k+1} - x^k \overset{(\alpha)}{=} f(\bar{x}^k) \Delta t + g(\bar{x}^k) \Delta\omega^k, $$其中 $\bar{x}^k = \alpha x^{k+1} + (1-\alpha) x^k$. 离散化的方式最终不影响 $\mathcal{P}[W_\tau]$ 的表达式，因为在连续极限下 $f(\bar{x}^k) \to f(x, t)$，$g(\bar{x}^k) \to g(x, t)$ 对任意 $\alpha \in [0, 1]$ 成立。变量代换可得$$ \mathcal{P}[W_\tau] = \prod_{k=0}^{K-1} \frac{1}{\sqrt{4\pi D \Delta t}} \exp\left\{ -\frac{1}{4D\Delta t} \left[ \frac{x^{k+1} - x^k - f(\bar{x}^k)\Delta t}{g(\bar{x}^k)} \right]^2 \right\}. $$轨迹不同离散化给出的轨迹概率的区别在于 Jacobi 行列式 $\left| J^{-1} \right|$：$$ \mathcal{P}[X_\tau] = p(x^0) \cdot \mathcal{P}[W_\tau] \cdot \left| J^{-1} \right|. $$在一般的 $\alpha$ 离散化下，Jacobi 行列式的计算结果为$$ \begin{aligned} \left| J^{-1} \right| ={}& \prod_{k=0}^{K-1} \frac{1}{|g(\bar{x}^k)|} \exp\left\{ -2\alpha\Delta\omega^k g'(\bar{x}^k) - \alpha\Delta t f'(\bar{x}^k) - \alpha\Delta t f(\bar{x}^k)\frac{g'(\bar{x}^k)}{g(\bar{x}^k)}\right. \\ &- \left. D\alpha^2\Delta t \left[ 2(g'(\bar{x}^k)^2 - g(\bar{x}^k) g''(\bar{x}^k) \right] \right\}. \end{aligned} $$当 $\alpha = 0$ 时，Jacobi 行列式的计算结果与第一个例子相吻合。对于可加的噪声，有 $g'(x) = 0$，此时 Jacobi 行列式简化为$$ \left| J^{-1} \right| = \prod_{k=0}^{K-1} \exp\left[ -\alpha f'(\bar{x}^k)\Delta t \right], $$这与第二个例子相吻合。Stratonovich 形式下的轨迹泛函有时也被称为 Onsager-Machlup 泛函，因为 Lars Onsager 和 Stefan Machlup 最早做出了相关的工作。取 $\alpha = 1/2$，此时 Onsager-Machlup 作用量为$$ \begin{aligned} \mathcal{S}[X_\tau] ={}& \underbrace{ \sum_{k=0}^{K-1} \frac{1}{4D\Delta t} \left[ \frac{x^{k+1} - x^k - f(\bar{x}^k)\Delta t}{g(\bar{x}^k)} \right]^2 }_{\text{Gaussian noise weight}} \\ &+ \underbrace{ \sum_{k=0}^{K-1} \frac{\Delta t}{2}\left[ f'(\bar{x}^k) - \frac{f(\bar{x}^k)g'(\bar{x}^k)}{g(\bar{x}^k)} \right] + \sum_{k=0}^{K-1} \frac{D\Delta t}{4}\left[ 2(g'(\bar{x}^k))^2 - g(\bar{x}^k) g''(\bar{x}^k) \right] }_{\text{Jacobian weight}}. \end{aligned} $$
考虑定义在同一个概率空间 $(\Omega, \mathscr{F}, P)$ 上的两个随机过程 $X$ 和 $Y$。若对于任意时间 $t\ge 0$ 和样本 $\omega \in \Omega$ ，都有 $X_t(\omega) = Y_t(\omega)$，那么我们称随机过程 $X$ 和 $Y$ 相等。在概率测度 $P$ 的意义下，我们还可以定义一些更弱的“相等”：若对任意 $t \ge 0$，有 $P[X_t = Y_t] = 1$，则称随机过程 $Y$ 是 $X$ 的修正；过程 $X$ 和 $Y$ 不可分辨，若它们几乎所有的轨迹相同：$$ P[X_t = Y_t; \forall 0\le t < \infty] = 1. $$“不可分辨”的条件比“修正”的条件更强，因为“不可分辨”要求几乎所有的轨迹相同，在时间尺度上是一种整体性质；而“修正”仅要求轨迹在某一时刻相同，在时间尺度上是点的性质。为了更进一步说明，我们可以考虑如下例子：考虑具有连续密度函数的正值随机变量 $T$，考虑如下随机过程：$$ X_t \equiv 0, \quad Y_t = \displaystyle{\begin{cases}0, t \neq T; \\ 1, t = T.\end{cases}} $$此时，$Y$ 是 $X$ 的修正，因为对于任意 $t \ge 0$，$P[Y_t = X_t] = P[t \neq T] = 1$。然而，$Y$ 与 $X$ 并非不可分辨，因为它们任意随机轨迹的实现都不一样： $P[Y_t = X_t; \forall t \ge 0] = 0$.如果我们进一步要求随机轨迹几乎处处右连续，那么“修正”可以推出“不可分辨”。
1. Gillespie 算法与连续时间 Markov 过程2. 算法原理3. 算法细节3.1. Python 生成器实现算法1. Gillespie 算法与连续时间 Markov 过程Gillespie 算法是一种可以用于采样连续时间 Markov 过程（CTMC）的随机轨迹的算法。这一算法被广泛用于模拟化学反应网络、生物物理过程等随机过程。在介观热力学中，Gillespie 算法被广泛应用于模拟主方程描述的系统，即连续时间 Markov 系统。对于连续时间、离散状态的马尔可夫过程，其概率演化可以由如下主方程表达$$ \frac{\mathrm d p(x, t)}{\mathrm dt} = R p(x, t). $$其中矩阵 $R$ 的矩阵元满足$$ \begin{cases} R_{ij} \ge 0, \quad i \neq j; \\ \\ R_{ii} = \displaystyle{-\sum_{j\neq i}}R_{ij}. \end{cases} $$列和为零的性质保证的概率的守恒（即归一性）。上式表明我们可以将连续时间 Markov 过程视为概率密度在各个状态上的流动过程。实际上，我们也可以考虑大量按初始概率分布的无相互作用粒子，粒子在不同状态之间的转移速率由转移矩阵 $R$ 给出。每个粒子在所有状态上的运动都构成一条随机轨迹，每一时刻所有粒子的分布情况即为该时刻系统概率密度的演化情况。Gillespie 算法即是基于上述的第二种视角，该算法可以模拟单个粒子在系统各个状态上的随机轨迹，由此得到的随机轨迹是上述主方程的解。随机轨迹本身给出系统的动力学性质，重复对轨迹进行采样，又可由轨迹平均得到系统的一系列统计性质。因此，Gillespie 算法之于介观热力学的重要性不言而喻。2. 算法原理对于主方程描述的系统，粒子在某一状态上向其相邻状态的转移概率于转移矩阵 $R$ 有关。一般地，定义矩阵元 $R_{ij}$ 表示状态 $j$ 向状态 $i$ 的转移速率。假设在极短的时间内，粒子仅可能发生一步跳跃，即发生多步跳跃的概率为时间的高阶无穷小量。这样的概率模型实际上是 Poisson 过程，即粒子每一步的跳跃过程都是一个 Poisson 过程。从状态 $j$ 到状态 $i$ 的 Poisson 过程的等待时间满足 Poisson 分布$$ P_{i\leftarrow j}(T \le t) = 1 - e^{-R_{ij}t}. $$对于一般的 Markov 系统，与某一状态 $i$ 相邻的状态可能有很多，粒子跳出状态 $i$ 的概率应当为跳到任意与 $i$ 相邻的状态的概率之和。由于 Poisson 分布之和依然是 Poisson 分布，因此总的跳出过程依然是 Poisson 过程，且其速率为所有可能的跳出速率之和。由矩阵 $R$ 的矩阵元定义可知，所有从状态 $i$ 跳出的速率之和即为 $-R_{ii}$。由上述分析可知，粒子在状态 $i$ 的等待时间满足分布$$ P(T \le t) = 1 - e^{R_{ii}t}. $$在等待时间之后，粒子跳到不同状态的概率与各个转移速率成正比。待粒子到了新的状态，再重复上述过程，即可得到一条随机轨迹。3. 算法细节每一步跳跃都涉及到两个随机的因素，即随机的等待时间和随机的跳跃方向。我们可以通过生成两个随机数来实现这两个过程。粒子在状态 $i$ 上的等待时间可以由一个在 $[0, 1)$ 上均匀分布的随机数 $r_1$ 生成。即$$ \text{waiting time} = \frac{1}{R_{ii}}\ln r_1. $$粒子下一步所处的位置可以由另一个在 $[0, 1)$ 上均匀分布的随机数 $r_2$ 生成。粒子下一步所处的位置为满足如下表达式的位置 $j'$$$ \frac{1}{-R_{ii}}\sum_j^{j'}R_{ji} < r_2 <\frac{1}{-R_{ii}}\sum_j^{j'+1}R_{ji}. $$这表明下一步处在某一状态的概率与其对应的转移速率成正比。算法有如下伪代码Step 1. Initialize t = t0, x = x0 Step 2. Calculate R_ji and R_ii Step 3. Generate waiting time with R_ii and r_1 Step 4. Generate next state with R_ii, R_ji and r_2 Step 5. Record results, return to step 1 or end the loop3.1. Python 生成器实现算法生成器是一种特殊的 Python 结构，它可以当作普通的可迭代对象使用，比如用在循环中。相比于一个普通的循环函数而言，生成器更加的节省内存空间，代码更加简洁优美，且具有更好的泛用性和拓展性。同时，Numba 等即时编译器也支持生成器的使用，为代码运行速度提供了保障。综上所述，我们选用生成器来实现 Gillespie 算法。对于任意矩阵 $R$ 和初始概率分布 $p_{\text{ini}}$，如下算法构建了一个随机轨迹生成器。import numpy as np import numba as nb @nb.jit(nopython=True) def traj(R, p_ini): # total number of states n = len(R[0, :]) # modify R with diagonal elements 0 for i in range(n): R[i, i] = 0 # choose a state based on the initial prob distribution r_1 = np.random.rand() state = np.searchsorted(np.cumsum(p_ini), r_1, side="right") # calculate waiting time with a random number a_i = R[:, state] sum_a_i = np.sum(a_i) r_2 = np.random.rand() waiting_time = -np.log(r_2) / sum_a_i # yield current state and waiting time yield state, waiting_time # start the loop while True: # generate the next state r_1 = np.random.rand() state = np.searchsorted(np.cumsum(a_i/sum_a_i), r_1, side="right") # calculate the waiting time on the next state a_i = R[:, state] sum_a_i = np.sum(a_i) r_2 = np.random.rand() waiting_time = -np.log(r_2) / sum_a_i # yield next state and waiting time yield state, waiting_time上述代码又两部分组成。第一部分在 while True 之前，这部分由初始概率分布生成处在某一位置的粒子，并计算粒子在该状态上的等待时间。第二部分在 while True 之后，由一个循环构成，每次使用 next() 函数调用该生成器都将进行一步迭代，生成一次跳跃过程及其对应的等待时间。函数开头的修饰器 @nb.jit(nopython=True) 调用了即时编译器 Numba，它通过即使编译该函数的方式大幅提高其运行速度，代价是函数内容需要仔细编写以满足 Numba 类型推断器的要求。可以使用 Python 内置的 next() 函数调用该生成器。例如可以用如下代码生成一段时间在 $[0, 100]$ 的随机轨迹T = 100 x = [] t = [] current_time = 0 trajectory = traj(R, p_ini) while current_time < T: state, waiting_time = next(trajectory) x.append(state) current_time += waiting_time t.append(current_time)
1. Wiener 过程2. 白噪声3. 随机积分4. Ito 微积分5. 随机微分方程6. Fokker-Planck 方程1. Wiener 过程首先考虑一个离散的一维随机漫步 (random walk)，记过程的起始点为 $x=0$，每一步行走距离为 $\Delta x$，向左向右走的概率均为 $1/2$，粒子每过 $\Delta t$ 的时间就移动一次。用 $X(t)$ 记录粒子在 $t$ 时刻的位置，则有$$ X(t) = \sum_{i=1}^{t/\Delta t} \eta_i, $$其中$$ \eta_i = \begin{cases} +\Delta x, \text{ with prob. 1/2 }; \\ -\Delta x, \text{ with prob. 1/2 }. \end{cases} $$由中心极限定理可知，当 $t$ 很大、或者说 $\Delta t$ 很小的时候，$X(t)$ 的分布趋近于高斯分布。此时，我们只需要知道 $X(t)$ 的均值和方差即可确定 $X(t)$ 的概率分布。均值很容易计算，即$$ \langle X(t) \rangle = \frac{t}{\Delta t}\langle\eta_i\rangle = 0; $$其方差为$$ \sigma^2_{X(t)} = \langle X(t)^2\rangle = \frac{t}{\Delta t}\sigma^2_{\eta_i} = \frac{(\Delta x)^2}{\Delta t}t. $$因此有 $X(t)\sim \mathcal N(0, \sqrt{(\Delta x)^2 t / \Delta t})$，由此我们可以计算很多与该过程有关的性质。实际上，不严谨地说，Wiener 过程可以看作上述离散随机游走的连续极限，即取 $\Delta x \to 0$，$\Delta t \to 0$ 。为了使 $(\Delta x)^2 / \Delta t$ 为一个有限值，$\Delta x$ 需要与 $\sqrt{\Delta t}$ 同阶。记 $\Delta x^2 / \Delta t = \sigma^2$，由此得到的连续过程即为 Wiener 过程 $W(t) \sim \mathcal N(0, \sigma\sqrt{t})$ 。从数学角度严格地说，Wiener 过程定义为满足下面三个性质的随机过程：$W(t) \sim \mathcal N(0, \sigma\sqrt{t})$；$W(t)$ 是马尔可夫过程；对于任意 $t' \ge 0$， $W(t+t') - W(t)$ 与 $W(t)$ 相互独立。条件 (3) 称为独立增量，它描述了过程的空间均一性，即过程与起点无关。由 Wiener 过程的定义即可计算 Wiener 过程的时间关联函数。由性质 (1)、(3) 可知 $\langle [W(t+t') - W(t)]W(t) \rangle = 0$，因此$$ \langle W(t+t')W(t)\rangle = \langle W(t)^2\rangle = \sigma^2t. $$故 Weiner 过程在时间 $t_1$ 、$t_2$ 的时间关联函数为$$ \langle W(t_1)W(t_2)\rangle = \sigma^2\min\{ t_1, t_2 \}. $$这一结论也可以从离散的随机游走看出来，即$$ \begin{aligned} \langle X(t_1)X(t_2) \rangle &= \sum_{i=1}^{t_1/\Delta t}\sum_{j=1}^{t_2/\Delta t}\langle \eta_i\eta_j\rangle \\ &= \Delta x^2\min\left\{ \frac{t_1}{\Delta t}, \frac{t_2}{\Delta t} \right\}. \end{aligned} $$2. 白噪声白噪声定义为 Weiner 过程的导数：$$ \xi(t) \equiv \frac{\mathrm d W(t)}{\mathrm d t}. $$这个导数在数学分析层面无法良好定义，在此仅将其形式上写作这样。计算它的均值可得$$ \langle \xi(t)\rangle = \frac{\mathrm d }{\mathrm d t}\langle W(t)\rangle = 0, $$时间关联函数为$$ \begin{aligned} \langle\xi(t_{1})\xi(t_{2})\rangle &= \left\langle\frac{\mathrm d}{\mathrm dt_{1}}W(t_{1})\frac{\mathrm d}{\mathrm dt_{2}}W(t_{2})\right\rangle \\ &= \frac{\mathrm d}{\mathrm dt_{1}}\frac{\mathrm d}{\mathrm dt_{2}}\sigma^{2} \min\{t_{1},t\} \\ &= \sigma^2\frac{\mathrm d }{\mathrm d t_1}\Theta(t_1 - t_2) \\ &= \sigma^2\delta(t_1 - t_2). \end{aligned} $$其中 $\Theta$ 是阶梯函数，它的导数是 $\delta$ 函数。由关联函数可知，$\xi(t_1)$ 与 $\xi(t_2)$ 相互独立 ($t_1 \neq t_2$)，但当 $t_1 = t_2$ 时，$\langle \xi (t)^2 \rangle$ 发散。这样的 $\xi(t)$ 我们称为广义随机过程。在离散情况下，我们可以算出$$ \frac{X(t') - X(t)}{t'-t} = \frac{1}{t'-t}\sum_{i=1}^{(t'-t)/\Delta t} \eta_i \sim \mathcal N\left(0, \frac{\sigma}{\sqrt{t'-t}}\right). $$这一形式可以用于模拟白噪声。我们将这一噪声称为白噪声，是因为其功率谱密度 $S_\xi(\omega)$ 为常数$$ S_\xi(\omega) = \int_{-\infty}^\infty\mathrm d t~e^{-\mathrm i\omega t}\langle\xi(t)\xi(0)\rangle = \sigma^2, $$这与白光的谱密度形式相同。我们可以将白噪声看作色噪声的特殊形式。若一个噪声 $Z(t)$ 的时间关联函数不是 $\delta$ 函数，那么该噪声就是色噪声。一般，我们将色噪声 $Z(t)$ 的时间关联函数写为$$ \langle Z(t)Z(t')\rangle = \frac{\sigma^2}{2\tau}e^{-\frac{|t-t'|}{\tau}}. $$白噪声可以视为色噪声在 $\tau \to 0$ 的极限情况。3. 随机积分考虑随机微分方程$$ \dot{x}(t) = f(x(t)) + g(x(t))\xi(t), $$其中 $\xi(t)$ 为高斯白噪声。由于噪声项的出现，这一方程为随机微分方程。这一方程有微分形式的表达式，即$$ \mathrm d x(t) = f(x(t))\mathrm d t + g(x(t))\mathrm d W(t), $$其中 $W(t)$ 是 Weiner 过程。对上述方程积分，形式上可以得到$$ x(t) = x(0) + \int_0^t\mathrm d t'~f(x(t')) + \int_0^t\mathrm d t'~g(x(t'))\mathrm d W(t'). $$这一形式告诉我们，要解随机微分方程，就要先定义随机积分。然而，随机积分并不容易定义。它的定义取决于不同表示方式的选取。考虑一个一般的随机积分$$ \int_0^\tau \psi(t)\mathrm d W(t). $$黎曼积分定义为如下求和的极限$$ \sum_{k=1}^{\tau/\Delta t}\psi(t_k^*)[W(t_k) - W(t_{k_1})], $$其中 $t_k = k\Delta t$，$t_k^* \in [t_{k-1}, t_k]$ 。黎曼积分的良好性质在于，积分的值与 $t_k^*$ 的选取无关，我们可以选区间内的任意一点。然而，对于随机积分，选取不同的 $t_k^*$ 会使得最终的积分结果不同，因此随机积分的结果取决于表示形式。有两种经常被选用的表示，一种是 Ito 表示，一种是 Stratonovich 表示：$$ \begin{cases} t_k^* = t_{k-1}: \text{ Ito}, \\ t_k^* = \frac{1}{2}(t_{k-1}+t_k) = \left(k-\displaystyle\frac{1}{2}\right)\Delta t: \text{Stratonovich}. \end{cases} $$我们用一个例子来展示两种表示的不同积分结果。考虑随机积分$$ \int_0^\tau \mathrm d t~W(t)\mathrm d W(t). $$对于 Ito 表示，我们可以计算积分结果的均值为$$ \begin{aligned} & \left\langle \sum_{k=1}^{\tau/\Delta t}W(t_{k-1})[W(t_{k})-W(t_{k-1})]\right\rangle \\ ={}& \sum_{k=1}^{\tau/\Delta t}\left\langle W(t_{k-1})\right\rangle\left\langle [W(t_{k})-W(t_{k-1})]\right\rangle \\ ={}& 0. \end{aligned} $$对于 Stratonovich 表示，可得$$ \begin{aligned} & \left\langle \sum_{k=1}^{\tau/\Delta t}W(t_k^*)[W(t_{k})-W(t_{k-1})]\right\rangle \\ ={}& \sum_{k=1}^{\tau/\Delta t}\left\langle W(t_k^*)W(t_{k})\right\rangle - \left\langle W(t_k^*)W(t_{k-1})\right\rangle \\ ={}& \sum_{k=1}^{\tau/\Delta t}(\sigma^2 t_k^* - \sigma^2 t_{k-1}) \\ ={}& \sigma^2\sum_{k=1}^{\tau/\Delta t}(t_k^* - t_{k-1}) \\ ={}& \frac{\sigma^2}{2}\tau. \end{aligned} $$其中第二个等号用了 $W(t)$ 的时间关联函数表达式。容易看出，这两种表示给出的计算结果是不同的。现在可以与微积分中的计算方式做一个对比。对于确定函数的微积分，积分结果应当为$$ \begin{gathered} \int_0^\tau W(t)\mathrm d W(t) = \left.\frac{W(t)^2}{2}\right|_{0}^\tau = \frac{W(\tau)^2}{2} \\ \implies \left\langle \frac{W(t)^2}{2} \right\rangle = \frac{\sigma^2}{2}\tau. \end{gathered} $$实际上，可以证明，Stratonovich 表示与确定函数的微积分给出的结果相同，而 Ito 表示则有所不同。计算 Stratonovich 表示的随机微积分时，我们可以直接使用链式法则。4. Ito 微积分为什么传统的微积分对伊藤积分无效？从前面对于离散情况的分析我们已经可以略知一二，即 $(\Delta x)^2$ 与 $\Delta t$ 同阶。我们可以从离散情况出发推导 Ito 表示下的 $W(t)\mathrm d W(t)$ 是什么。$$ \begin{aligned} {}&W(t) [W(t+\Delta t) - W(t)] \\ ={}&\frac{1}{2}[W(t+\Delta t)^2 - W(t)^2] - \frac{1}{2}[W(t+\Delta t) - W(t)]^2, \end{aligned} $$而 $\mathrm d W(t)^2 / 2$ 的离散表达式为$$ \frac{\mathrm d W(t)^2}{2} = \frac{1}{2}[W(t+\Delta t)^2 - W(t)^2], $$若二阶项不能忽略，则它们不一样。我们已经计算出了二阶项的均值为$$ \langle [W(t+\Delta t) - W(t)]^2 \rangle = \sigma^2\Delta t, $$可以发现，Weiner 过程的二阶项与时间的一阶项同阶，这就是随机分析与普通微积分的不同之处。伊藤微积分的思想便是，保持式中的二阶项 $(\mathrm d W)^2$，并将其记为一阶项 $\sigma^2\mathrm d t$ 。由此，我们可以得到随机微分方程的变量代换公式，该公式由伊藤引理给出。令 $x(t)$ 为随机微分方程$$ \mathrm d x = f(x)\mathrm d t + g(x)\mathrm d W $$在 Ito 表示下的解，那么 $y(t) = y(x(t))$ 是如下方程的解：$$ \begin{aligned} \mathrm d y &= \frac{\mathrm d y}{\mathrm d x}\mathrm d x + \frac{1}{2}\frac{\mathrm d^2y}{\mathrm d x^2}(\mathrm d x)^2 \\ &= \frac{\mathrm d y}{\mathrm d x}[f(x)\mathrm d t + g(x)\mathrm d W] + \frac{1}{2}\frac{\mathrm d^2 y}{\mathrm d x^2}[f(x)\mathrm d t + g(x)\mathrm d W]^2 \\ &= \left[\frac{\mathrm d y}{\mathrm d x}f(x)+\frac{\sigma^2}{2}\frac{\mathrm d^2y}{\mathrm d x^2}g(x)^2\right]\mathrm d t + \frac{\mathrm d y}{\mathrm d x}g(x)\mathrm d W. \end{aligned} $$上述公式可以简记为$$ \mathrm d y = \left(y'f+\frac{\sigma^2}{2}y''g^2\right)\mathrm d t + y'g\mathrm d W. $$我们称 $\frac{\sigma^2}{2}y''g^2\mathrm d t$ 为 Ito 项。5. 随机微分方程在解随机微分方程之前，我们需要指明选用的表示形式，因为不同的表示会给出不同的积分结果。对于一个一般的随机微分方程$$ \dot x = f(x) + g(x)\xi, $$为了区分不同的表示方法，通常我们记 Ito 表示为 $g(x)\xi$，记 Stratonovich 表示为 $g(x)\circ \xi$ 。对应的微分表达为$$ \begin{cases} \mathrm d x = f(x)\mathrm d t + g(x)\mathrm d W: \text{Ito interpertation}; \\ \mathrm d x = f(x)\mathrm d t + g(x)\circ \mathrm d W: \text{Stratonovich interpertation.} \end{cases} $$若我们考虑色噪声，再令 $\tau \to 0$，那么将会得到 Stratonovich 表示下的白噪声。若 $g(x)$ 为常数，则 Ito 表示和 Stratonovich 表示给出的解相同。Langevin 方程就是 $g(x)$ 为常数的情况，因此以往的处理并不需要区分这两种不同的表示形式。然而，若我们想要在 Langevin 系统上定义热力学量，那么我们需要选择 Stratonovich 表示。对于 Ito 积分，我们需要使用伊藤引理，它非常适合描述离散时间系统，例如股票等；对于 Stratonovich 积分，我们可以直接使用与普通微积分相同的形式，它适用于描述连续时间系统，例如物理和化学系统。对于速度为 $v(t)$ 的布朗粒子，其运动方程为（欠阻尼） Langevin 方程$$ m\dot v(t) = -\gamma v(t) + \xi(t). $$其中 $\xi(t)$ 为高斯白噪声。此时两种表述给出的解相同。由于没有外力，系统的能量即为粒子的动能$$ E(t) = \frac{1}{2}mv(t)^2. $$在 Stratonovich 表示下能量满足的运动方程为$$ \dot E(t) = mv(t)\dot v(t) = -\gamma v(t)^2 + v(t)\circ\xi(t). $$在 Ito 表示下，使用 Ito 引理，运动方程可以写为$$ \dot E(t) = mv\dot v + \frac{\sigma^2}{2}m\left( \frac{1}{m} \right)^2 = -\gamma v^2+\frac{\sigma^2}{2m} + v\xi. $$可以看出，Stratonovich 表示给出的能量变化率的均值为零，而原 Langevin 方程对于这两种表示会给出相同的结果，因此 Ito 表示下的能量运动方程的均值也应该为零。故有$$ \langle \dot E(t)\rangle = -\frac{2\gamma}{m}\langle E(t)\rangle + \frac{\sigma^2}{2m} = 0. $$由此即可得到$$ \frac{\sigma^2}{2m} = \frac{2\gamma}{m}\frac{kT}{2} \implies \sigma^2 = 2\gamma kT. $$这就是涨落耗散定理的最简单形式。下面我们给出两种表示下随机微分方程之间的关系。若 $x(t)$ 是方程$$ \dot x(t) = f(x)+g(x)\xi $$的解，那么 $x(t)$ 也是方程$$ \dot{x}(t) = f(x) - \frac{\sigma^2}{2}g(x)g'(x) + g(x)\circ \xi $$的解。证明：令$$ y(x) = \int\frac{\mathrm d x'}{g(x')}, \quad \frac{\mathrm d y}{\mathrm dx} = \frac{1}{g(x)}; \quad \frac{\mathrm d^2 y}{\mathrm d x^2} = -\frac{g'(x)}{g(x)^2}. $$由伊藤引理可得$$ \begin{aligned} \dot y &= \frac{1}{g}(f+g\xi) + \frac{\sigma^2}{2}\left( -\frac{g'}{g^2} \right)g^2 \\ &= \frac{f}{g} - \frac{\sigma^2}{2}g' + \xi, \end{aligned} $$再由普通的微积分可得$$ \begin{aligned} \dot x(t) &= g(x)\dot y \\ &= g\left( \frac{f}{g} - \frac{\sigma^2}{2}g' + \xi \right) \\ &= f - \frac{\sigma^2}{2}gg' + g\circ \xi. \end{aligned} $$反过来也有类似的结论，若 $x(t)$ 是方程$$ \dot x(t) = f(x)+g(x)\circ\xi $$的解，那么 $x(t)$ 也是方程$$ \dot{x}(t) = f(x) + \frac{\sigma^2}{2}g(x)g'(x) + g(x)\xi $$的解。注意此定理与上述定理的正负号不同。6. Fokker-Planck 方程Fokker-Planck 方程是概率密度 $p(x, t)$ 的演化方程。首先处理 Ito 表示。对于一般的随机微分方程$$ \dot x(t) = f(x(t)) + g(x(t))\xi(t), $$考虑一个任意函数 $A(t) = A(x(t))$，使用伊藤引理，$A(t)$ 的运动方程为$$ \dot A(t) = A'(x(t))[f(x(t)) + g(x(t))\xi(t)] + \frac{\sigma^2}{2}A''(x(t))g(x(t))^2. $$对上述方程求期望可得$$ \langle \dot A(t)\rangle = \langle A' f\rangle + \frac{\sigma^2}{2}\langle A'' g^2\rangle. $$而 $A(t)$ 的期望为$$ \langle A(x(t)) \rangle = \int\mathrm d x~\rho(x, t)A(x). $$对上式求时间导数，再对照先前的式子可得$$ \int\frac{\partial \rho}{\partial t}A = \int\mathrm d x~\rho\left( fA' + \frac{\sigma^2}{2}g^2 A'' \right). $$对上式右边的两项做分部积分。第一项分部积分得到$$ \int_a^b \mathrm d x~\rho(x, t)f(x)A'(x) = -\int\mathrm d x~[\rho(x, t)f(x)]'A(x), $$其中分部积分得到的第一项由于边界条件而为零。对第二项做两次分部积分，同理可得$$ \int\mathrm d x~\rho(x, t)g^2(x) A''(x) = \int\mathrm d x~[\rho(x, t)g^2(x)]''A(x). $$将这两项放在一起可得$$ \int \mathrm d x~A(x)\left[ \frac{\partial \rho}{\partial t} + \frac{\partial f\rho}{\partial x} - \frac{\sigma^2}{2}\frac{\partial^2g^2\rho}{\partial x^2} \right] = 0. $$此式对任意 $A(x)$ 都成立，因此有关于 $\rho(x, t)$ 的方程$$ \frac{\partial \rho(x, t)}{\partial t} = -\frac{\partial}{\partial x}[f(x)\rho(x, t)] + \frac{\sigma^2}{2}\frac{\partial^2}{\partial x^2}[g(x)^2 \rho(x, t)]. $$此式即为伊藤表示下的一维 Fokker-Planck 方程。Stratonovich 表示下的 Fokker-Planck 方程可以由上面的方程变换得到，其结果为$$ \begin{aligned} \frac{\partial \rho(x, t)}{\partial t} =&{} -\frac{\partial}{\partial x}[f(x)\rho(x, t)] + \frac{\sigma^2}{2}\frac{\partial}{\partial x}[g(x)g'(x)\rho(x, t)] \\ &{}+ \frac{\sigma^2}{2}\frac{\partial^2}{\partial x^2}[g(x)^2\rho(x, t)] \\ =&{} -\frac{\partial}{\partial x}[f(x)\rho(x, t)] + \frac{\sigma^2}{2}\frac{\partial}{\partial x}\left[ g(x)\frac{\partial}{\partial x}g(x)\rho(x, t) \right]. \end{aligned} $$若将 $\partial_x g$ 记为算符，那么上式可以写为$$ \partial_t \rho = -\partial_x(f\rho) + \frac{\sigma^2}{2}(\partial_x g)^2\rho. $$
1. 由 Langevin 动力学给出的随机热力学1.1 随机 Langevin 动力学1.2 熵1.3 涨落定理参考文献1. 由 Langevin 动力学给出的随机热力学1.1 随机 Langevin 动力学考虑过阻尼 Langevin 动力学，运动方程为$$ \dot{x} = \mu F(x, \lambda) + \zeta $$噪声$~\zeta~$满足$$ \langle \zeta(\tau) \rangle = 0, \quad \langle \zeta(\tau)\zeta(\tau') \rangle = 2D\delta(\tau - \tau') $$其中$~D~$与迁移率$~\mu~$满足爱因斯坦关系$~D = T\mu~$，$~T~$为介质的温度（方便起见这里将$~k_B~$记为$~1~$）。力$~F(x, \lambda)~$可分解为保守部分和非保守部分$$ F(d, \lambda) = -\partial_x V(x, \lambda) + f(x, \lambda). $$其中$~V(x, \lambda)~$是保守势场，$~f(x, \lambda)~$是直接作用于粒子的外力。$~\lambda~$是含时的外部参数，可以简记为$~\lambda(0) = \lambda_0, \lambda(t) = \lambda_t~$.Langevin 方程有等价的 Fokker-Planck 方程$$ \begin{aligned} \partial_\tau p(x, \tau) &= -\partial_x j(x, \tau) \\ &= \partial_x (\mu F(x, \lambda) - D\partial_x)p(x, \tau) \end{aligned} $$其中概率分布$~p(x, \tau)~$描述$~\tau~$时刻在$~x~$处找到粒子的概率，$~j(x, \tau)~$为概率流。在$~\tau = 0~$时刻，系统满足归一化的概率分布$~p(x, 0) \equiv p_0(x)~$.1.2 熵在随机轨迹的层面上定义系统的熵$$ s(\tau) = - \ln p(x(\tau), \tau). $$其中概率分布$~p(x, \tau)~$是 Fokker-Planck 方程的解。对所有路径求均值，得到宏观的系统熵为$$ S(\tau) \equiv -\int \mathrm d x p(x, \tau)\ln p(x, \tau). $$对于定义在轨迹层面的系统熵，其变化率为$$ \begin{aligned} \dot{s}(\tau) &= -\left.\frac{\partial_\tau p(x, \tau)}{p(x, \tau)}\right|_{x(\tau)} - \left.\frac{\partial_x p(x, \tau)}{p(x, \tau)}\right|_{x(\tau)}\dot{x} \\ &= -\left.\frac{\partial_\tau p(x, \tau)}{p(x, \tau)}\right|_{x(\tau)} + \left.\frac{j(x, \tau)}{D p(x, \tau)}\right|_{x(\tau)}\dot{x} - \left.\frac{\mu F(x, \lambda)}{D}\right|_{x(\tau)}\dot{x} \end{aligned} $$系统到介质的热耗散则与环境的熵（或称介质熵）有关$$ \dot{q}(\tau) = F(x, \lambda)\dot{x} \equiv T\dot{s}_\text{m}(\tau) $$根据$~D = \mu T~$，不难发现介质熵的变化率与系统熵的最后一项形式相同，符号相反。因此总熵的变化率为$$ \begin{aligned} \dot{s}_\text{tot}(\tau) &= \dot{s}_\text{m}(\tau) + \dot{s}(\tau) \\ &= -\left.\frac{\partial_\tau p(x, \tau)}{p(x, \tau)}\right|_{x(\tau)} + \left.\frac{j(x, \tau)}{D p(x, \tau)}\right|_{x(\tau)}\dot{x}. \end{aligned} $$第一项中$~p(x, \tau)~$随时间的变化可以由参数$~\lambda(\tau)~$引起，也可以是在$~\lambda~$不变的情况下，系统由非定态的初态向定态演化引起的。$~\dot{x}~$的平均为$$ \langle \dot{x} \mid x, \tau \rangle = \frac{j(x, \tau)}{p(x, \tau)}, $$由概率守恒有$$ \int\mathrm d x \partial_\tau p(x, \tau) = 0, $$因此总熵变化率的轨迹平均为$$ \dot{S}_\text{tot}(\tau) \equiv \langle\dot{s}_{\text{tot}}(\tau)\rangle = \int\mathrm dx\frac{j(x, \tau)^2}{D p(x, \tau)} \ge 0 $$这给出经典的热力学第二定律。上式当且仅当系统处于平衡态时取等。介质熵的变化率为$$ \dot{S}_\text{m}(\tau) \equiv \langle\dot{s}_\text{m}(\tau) \rangle = \frac{\langle F(x, \tau)\dot{x}\rangle}{T} = \int\mathrm d x\frac{1}{T}F(x, \tau)j(x, \tau). $$系统熵的变化率为$$ \dot{S}(\tau) \equiv \langle\dot{s}(\tau)\rangle = \dot{S}_\text{tot}(\tau) - \dot{S}_\text{m}(\tau). $$这些熵的变化率我们可以称为熵产生。1.3 涨落定理引入时间反演的参数$$ \lambda^\dagger(\tau) \equiv \lambda(t - \tau), $$逆向路径为$$ x^\dagger(\tau) \equiv x(t - \tau). $$起始位置为$$ x_0 \equiv x(0) = x^\dagger(t) \equiv x^\dagger_t, $$终止位置为$$ x_t \equiv x(t) = x^\dagger(0) \equiv x^\dagger_0. $$记正向路径的概率为$~p[x(\tau)\mid x_0]~$，逆向路径的概率为$~p^\dagger[x^\dagger(\tau)\mid x^\dagger_0]~$。对 Langevin 方程做路径积分得到$$ \ln\frac{p[x(\tau)\mid x_0]}{p^\dagger[x^\dagger(\tau)\mid x^\dagger_0]} = \int_0^t \frac{F(x, \tau)\dot{x}}{T}\mathrm d \tau = \Delta s_\text{m}. $$路径积分：考虑过阻尼 Langevin 方程$$ \dot{x}(\tau) = \mu F(x(\tau), \lambda(\tau)) + \zeta(\tau). $$将时间离散化，$~t \equiv i\epsilon,~ (i = 0, \cdots, N)~$，记时间格点上的位置和参数为$$ x_i = x(i\epsilon), \quad \lambda_i = \lambda(i\epsilon). $$因此离散化的 Langevin 方程为$$ \frac{x_i - x_{i-1}}{\epsilon} = \frac{\mu}{2}[F_i(x_i) + F_{i-1}(x_{i-1})] + \zeta_i. $$其中$~F_i(x_i) \equiv F(x_i, \lambda_i)~$，这是 Stratonovich 形式的随机微积分（即取区间中点的值）。噪声满足$$ \langle \zeta_i\rangle = 0, \quad \langle \zeta_i\zeta_j\rangle = 2\frac{D}{\epsilon}\delta_{ij}. $$路径权重为$$ p(\zeta_1, \cdots, \zeta_N) = \left( \frac{\epsilon}{4\pi D} \right)^{N/2} \exp\left[ -\frac{\epsilon}{4D}\sum_i \zeta_i^2 \right]. $$因此从$~x_0~$出发沿某条路径转移的概率为$$ p(x_1, \cdots, x_N\mid x_0) = \det\left( \frac{\partial \zeta_i}{\partial x_j} \right) p(\zeta_1, \cdots, \zeta_N). $$其中$~\partial\zeta_i / \partial x_j~$为 Jacobi 矩阵$$ \frac{\partial \zeta_i}{\partial x_j} = \begin{pmatrix} \frac{1}{\epsilon} - \frac{\mu}{2}F'_1(x_1) & 0 & \cdots & \cdots \\ -\frac{1}{\epsilon} - \frac{\mu}{2}F'_1(x_2) & \frac{1}{\epsilon} - \frac{\mu}{2}F'_2(x_2) & 0 & \cdots \\ \cdots & \cdots & \cdots & \cdots \end{pmatrix} $$相应的 Jacobi 行列式为$$ \begin{aligned} \frac{\partial \zeta_i}{\partial x_j} &= \left(\frac{1}{\epsilon}\right)^N \prod_{i=1}^N \left[ 1 - \frac{\epsilon\mu}{2}F'_i(x_i) \right] \\ &= \left( \frac{1}{\epsilon} \right)^N \exp\left[ \sum_{i=1}^N \ln\left( 1 - \frac{\epsilon\mu}{2}F'_i(x_i) \right) \right] \\ &\approx \left( \frac{1}{\epsilon} \right)^N \exp\left[ -\sum_{i=1}^N \frac{\epsilon\mu}{2}F'_i(x_i) \right]. \end{aligned} $$带入可得随机轨迹的概率为$$ \begin{aligned} p(x_1, \cdots, &{} x_N \mid x_0) = \frac{1}{(4\pi D\epsilon)^{N/2}} \\ &\times \exp\left\{ -\frac{1}{4D\epsilon}\left[ \sum_{i=1}^{N} (x_i - x_{i-1} - \epsilon\mu F_i(x_i))^2 \right] - \frac{\epsilon\mu}{2}\sum_{i=1}^N F'_i(x_i) \right\} \end{aligned} $$再取极限$~\epsilon \to 0, N \to \infty, N\epsilon = t~$，此时回到连续情形，得到路径概率为$$ p[x(\tau)\mid x_0] \equiv \exp\left[ -\frac{1}{4D}\int_0^t (\dot{x} - \mu F)^2\mathrm d \tau - \frac{\mu}{2}\int_0^t F'\mathrm d \tau \right] \equiv \exp\left\{ -\mathcal A[x(\tau)] \right\}. $$其中$~F = F(x(\tau), \lambda(\tau))~$，$~F' = F'(x(\tau), \lambda(\tau))~$。作用量$~\mathcal A~$为$$ \mathcal A[x(\tau)] \equiv \frac{1}{D}\int_0^t \mathrm d \tau L(x(\tau), \dot{x}(\tau); \lambda(\tau)). $$$~L~$为拉格朗日函数$$ L(x, \dot{x}, \lambda(\tau)) \equiv \frac{1}{4}(\dot{x} - \mu F)^2 + \frac{\mu D}{2}F'. $$对于所有从$~x_0~$开始的路径积分测度为$$ \int_{x_0} \mathrm d [x(\tau)] \equiv \lim_{\epsilon \to 0, N \to \infty, \epsilon N = t}\left( \frac{1}{4\pi D\epsilon} \right)^{N/2} \prod_{i=1}^N \int_{-\infty}^{+\infty}\mathrm d x_i. $$路径概率自然应该是归一化的，计算可得$$ \int_{x_0} \mathrm d [x(\tau)] p[x(\tau) \mid x_0] = \int \mathrm d [x(\tau)] p[x(\tau)\mid x_0] p_0(x_0) = 1. $$考虑逆向路径，有$$ p[x^\dagger(\tau)\mid x^\dagger_0] = \exp\left\{ -\mathcal A^\dagger [x^\dagger(\tau)] \right\}. $$其中$$ \mathcal A^\dagger [x^\dagger(\tau)] = \mathcal A[x(\tau)] + \frac{1}{T}\int_0^t \mathrm d \tau \dot{x}(\tau)F(x(\tau), \lambda(\tau)). $$因此可得$$ \frac{p[x(\tau)\mid x_0]}{p[x^\dagger(\tau)\mid x^\dagger_0]} = \exp\left\{\frac{q[x(\tau)]}{T}\right\} = e^{\Delta s_\text{m}}. $$正、逆向路径概率之比的对数正好与介质熵的变化量有关。现在对任意的归一化分布$~p_0(x_0), p_1(x_0^\dagger) = p_1(x_t)~$定义一个量$~R~$$$ \begin{aligned} R[x(\tau), \lambda(\tau); p_0, p_1] &\equiv \ln\frac{p[x(\tau)\mid x_0]p_0(x_0)}{p^\dagger[x^\dagger(\tau)\mid x^\dagger_0]p_1(x^\dagger_0)} \\ &= \Delta s_\text{m} + \ln\frac{p_0(x_0)}{p_1(x_t)}. \end{aligned} $$可以发现，$~e^{-R}~$对初始位置和所有轨迹求均值可得$$ \begin{aligned} \left\langle e^{-R} \right\rangle &\equiv \sum_{x(\tau), x_0} p[x(\tau)\mid x_0] p_0(x_0)e^{-R} \\ &= \sum_{x^\dagger(\tau), x^\dagger_0} p^\dagger[x^\dagger(\tau)\mid x^\dagger_0] p_1(x^\dagger_0) \\ &= 1 \end{aligned} $$若末态$~p_1(x_t) = p(x, t)~$为 Fokker-Planck 方程在初始条件为$~p_0(x_0)~$下的解，则有积分涨落定理$$ \left\langle e^{-\Delta s_{\text{tot}}} \right\rangle = 1. $$若非保守外力$~f = 0~$，且系统的初、末态都是平衡态$$ p_{0, 1}(x) = p^\text{s}(x, \lambda_{0, t}) = \exp[-V(x, \lambda_{0, t}) - \mathcal F(\lambda_0, t)] $$则$$ R = \Delta s_{\text m} + \frac{1}{T}[V(x_t, \lambda_t) - V(x_0, \lambda_0) - \mathcal F(\lambda_t) + \mathcal F(\lambda_0)] = \frac{w_\text d}{T}. $$带入得到 Jarzynski 恒等式$$ \left\langle e^{-w_\text d / T} \right\rangle = 1. $$对于有为常量的$~\lambda~$和不为零的非保守外力$~f~$的定态，令$~p_0(x) = p_1(x) = p^\text s(x)~$，则有$$ \frac{p(-R)}{p(R)} = e^{-R}, \quad \frac{p(-\Delta s_\text{tot})}{p(\Delta(s_{\text{tot}}))} = e^{-\Delta s_{\text{tot}}}. $$此关系式称为细致涨落定理。参考文献Seifert, Udo. "Stochastic Thermodynamics: Principles and Perspectives." European Physical Journal B, vol. 64, no. 3, 2008, pp. 423–431.Seifert, Udo. "Stochastic Thermodynamics, Fluctuation Theorems and Molecular Machines." Reports on Progress in Physics, vol. 75, no. 12, 2012, pp. 126001–126001.
1. 马氏过程1.1. 马氏链1.2. 多步转移概率1.3. 状态的分类1.4. 状态的判定方法1.5. 平稳分布1.6. 极限行为1.7. 一些例子1.8. 离出分布2. 连续时间马氏链2.1. 极限行为1. 马氏过程1.1. 马氏链考虑离散时间的随机过程$~X_n(n = 1, 2, \cdots)~$，$~X_n~$在有限集合内取值。我们称$~X_n~$的所有可能取值为系统的状态。定义转移概率为$$ p_{n+1}(i, j) = \mathbb P(X_{n+1} = j\mid X_n = i, X_{n-1} = i_{n-1}, \cdots, X_0 = i_0) $$转移概率是一个条件概率。这是一个离散时间，离散状态的随机过程。若$~p_{n+1}(i, j) = \mathbb P(X_{n+1} = j\mid X_n = i)~$，则称​这个过程具有马氏性（马尔可夫性）。马尔可夫性即是说，转移概率只与当前的状态有关，与先前任何时间的状态都无关。若时刻$~n~$​的取值与转移概率无关，即对任意$~n~$​都有$$ \mathbb P(X_{n+1} = j\mid X_n = i) = p(i, j) $$则称转移概率具有时齐性，相应的马氏链称作时齐马氏链。转移概率$~p(i, j)~$的相应数值，可以标注在矩阵$~\mathbf P~$的第$~i~$行，第$~j~$列。我们称这个矩阵为转移矩阵$$ \mathbf P = \begin{pmatrix} p(0, 0) & p(0, 1) & \cdots & p(0, N) \\ p(1, 0) & p(1, 1) & \cdots & p(1, N) \\ \vdots & \vdots & \ddots & \vdots \\ p(N, 0) & p(N, 1) & \cdots & p(N, N) \end{pmatrix} $$转移矩阵是一个$~(N+1)\times (N+1)~$矩阵。马氏链有如下性质：转移概率非负，即$~p(i, j) \ge 0~$.转移矩阵每行元素之和均为$~1~$，即$~\displaystyle\sum_{j}p(i, j) = 1~$.1.2. 多步转移概率转移概率$~p(i, j) = \mathbb P(X_{n+1} = j\mid X_n = i)~$给出了从状态$~i~$到状态$~j~$的一步转移概率，类似地，从$~i~$到$~j~$的$~m~$步$~(m > 0)~$转移概率定义为$$ p^m(i, j) = \mathbb P(X_{n+m} = j\mid X_n = i) $$可以证明，$~m~$步转移矩阵$~\mathbf P^m~$是$~1~$步转移矩阵$~\mathbf P~$的$~m~$次幂，$~p^m(i, j)~$是$~m~$步转移矩阵$~\mathbf P^m~$第$~i~$行，第$~j~$列的元素。我们有$~\text{C-K}~$​方程，即$~\text{Chapman-Kolmogorov}~$​​方程：$$ p^{m+n}(i, j) = \mathbb P(X_{m+n} = j\mid X_0 = i) = \sum_k p^m(i, k)\cdot p^n(k, j) $$由$~\text{C-K}~$方程，可以得到如下不等式$$ p^{m+n}(i, j) \ge p^m(i, k)\cdot p^n(k, j),\quad \forall k $$1.3. 状态的分类研究马氏链主要进行两类分析：瞬态分析和稳态分析。瞬态分析是研究在某一固定时刻$~n~$，马氏链对应系统的概率特征，即求$~n~$步转移概率。稳态分析则是研究当$~n\to \infty~$​时，马氏链对应系统的概率特征。记几个符号初始时刻状态为$~x~$的条件下，事件$~A~$发生的概率$~\mathbb P_x(A) = \mathbb P(A\mid X_0 = x)~$.首次返回状态$~y~$的最短时间$~T_y = \min\{n: n\ge 1, X_n = y\}~$​. $~T_y~$可看作一个停时。从初始时刻状态$~y~$经过有限步首次返回到状态$~y~$的概率$~\rho_{yy} = \mathbb P_y(T_y < \infty)~$.由于当$~n = 0~$时$~\rho_{y}\equiv 1~$，故为避免此情况，我们规定$~n \ge 1~$. 不难看出，从初始时刻状态$~y~$返回状态$~y~$的次数为$~k~$的概率为$~\rho_{yy}^k~$.当某事件在时刻$~n~$停止，记$~\{T = n\}~$为停时。若$~T_y~$是一个停时（首次返回状态$~y~$的最短时间），则可以表示为$$ \{T_y = n\} = \{X_n = y, X_{n-1}\neq y,\cdots, X_1 \neq y\} $$停时也可以满足马氏性。我们称这样的马氏性为强马氏性。设$~T~$是停时，假定$~T = n, X_T = y~$，则关于$~X_0, \cdots, X_T~$​的其他信息与未来的预测无关，并且$~X_{T + k}, k \ge 0~$的行为与初始状态为$~y~$的马氏链相同。即：$$ \begin{aligned} \mathbb P(X_{T+1} &= z\mid X_T = y, T = n) = \mathbb P(X_1 = z\mid X_0 = y) = p(y, z) \\ \mathbb P(X_{T+k} &= z\mid X_T = y, T = n) = \mathbb P(X_k = z\mid X_0 = y) = p^k(y, z), \quad k\ge 1 \end{aligned} $$马氏性中，时间$~n~$的取值是确定的；在强马氏性中，停时$~T~$​​​是随机的。对于随机的$~T~$，马氏性，仍成立，条件更强，故称强马氏性。定义符号$~T_y^k = \min\{n : n > T_y^{k-1}, X_n = y\}~$为第$~k~$次返回状态$~y~$的最短时间；$~\mathbb P_y(T_y^k < \infty) = \rho_{yy}^k~$为经过有限时间返回状态$~y~$次数为$~k~$的概率。若$~\rho_{yy} < 1~$，此时随着$~k\to \infty~$，有$~\rho_{yy}^k \to 0~$. 这意味着最终马氏链不再回到状态$~y~$，该状态$~y~$称为非常返态，也称暂态。若$~\rho_{yy} = 1~$，此时$~\rho_{yy}^k = 1~$. 这意味着马氏链回到状态$~y~$的次数有无穷多次，该状态$~y~$称为常返态。若$~\rho_{yy} < 1~$​，则系统可能无法从状态$~y~$​再次回到状态$~y~$​；若$~\rho_{yy} = 1~$​​，则系统一定能在有限时间从状态$~y~$​返回状态$~y~$​.常返态有一个特例，即吸收态：如果$~p(y, y) = 1~$，即$~\mathbb P_y(T_y = 1) = 1~$，则$~y~$是吸收态，此时$~y~$是非常强的常返态，因为链将会一直留在那里。这两种性质有等价的表达式：非常返：状态$~y~$有正的概率不返回$~y~$$$ \rho_{yy} < 1 \iff \mathbb P_y(T_y < \infty) < 1 \iff \mathbb P_y(T_y = \infty) > 0 $$常返：状态$~y~$在有限时间内返回$~y~$的概率为$~1~$$$ \rho_{yy} = 1 \iff \mathbb P_y(T_y < \infty) = 1 \iff \mathbb P_y(T_y = \infty) = 0 $$1.4. 状态的判定方法若从状态$~x~$，有正的概率到达状态$~y~$，则称$~x~$可达$~y~$，记作$~x \to y~$，即$$ \rho_{xy} = \mathbb P_x(T_y < \infty) > 0 $$也可以表示为$$ \exist ~n,\text{ s.t. } p^n(x, y) > 0 $$即从状态$~x~$经过$~n~$步到达状态$~y~$的概率为正。若$~x~$可达$~y~$，并且$~y~$可达$~x~$，则称$~x~$与$~y~$互通，记作$~x \leftrightarrow y~$. 也可以表示为$$ \exist~m, n, \text{ s.t. }p^m(x, y)>0, p^n(y, x) > 0 \implies x \leftrightarrow y $$即从状态$~x~$经过$~m~$步到达状态$~y~$的概率为正，且从状态$~y~$经过$~n~$步到达状态$~x~$的概率也为正。若$~\rho_{xy} > 0~$，且$~\rho_{yx} < 0~$，则$~x~$​​是非常返态。互通具有传递性，自反性，对称性，可作为一种分类。对于任何互通类$~C~$而言，其中的所有状态均为常返态或非常返态。对于状态集合$~A~$，若$~i\in A, j\notin A \implies p(i, j) = 0~$，则称$~A~$为闭集。状态空间的最小闭集称为不可约闭集。若马氏链只有一个互通类，则称这个链不可约。若马氏链中的互通类超过一个，则称这个链可约。若对有限状态空间中的任意两个状态$~x~$和$~y~$，存在$~n\in\mathbb Z^+~$，使得$~p^n(x, y) > 0~$，则由此构造的转移矩阵$~\mathbf P~$​不可约。记$~N(y)~$为访问状态$~y~$的次数；$~\mathbb P_x(X_n = y)~$表示初始状态为$~x~$、第$~n~$时刻状态为$~y~$的概率，则从状态$~x~$到$~y~$的访问次数期望值$~\mathbb E_x[N(y)]~$为：$$ \mathbb E_x[N(y)] = \frac{\rho_{xy}}{1 - \rho_{yy}} = \sum_{n = 1}^{\infty}p^n(x, y) $$1.5. 平稳分布如果增加非周期性的假定条件，则一个有限状态的不可约马氏链会收敛于一个平稳分布，即$$ \lim_{n\to\infty} p^n(x, y) = \pi(y) $$这个分布$~\pi (y)~$与期初状态是无关的。用条件概率的定义，有$$ \begin{aligned} \mathbb P(X_n = j) &= \sum_i \mathbb P(X_0 = i, X_n = j) \\ &= \sum_i\mathbb P(X_0 = i)\mathbb P(X_n = j\mid X_0 = i) \\ &= \sum_i q(i)p^n(i, j) \end{aligned} $$其中$~q(i) = \mathbb P(X_0 = i)~$是初始概率，由$~q(i), i\in S~$组成的向量$~\mathbf q~$构成初始概率分布。写成矩阵形式为$$ \begin{aligned} \mathbf q\mathbf P^n &= \begin{pmatrix} q(1) & q(2) & \cdots & q(k) \end{pmatrix}\begin{pmatrix} p^n(1, 1) & p^n(1, 2) & \cdots & p^n(1, k) \\ p^n(2, 1) & p^n(2, 2) & \cdots & p^n(2, k) \\ \vdots & \vdots & \ddots & \vdots \\ p^n(k, 1) & p^n(k, 2) & \cdots & p^n(k, k) \end{pmatrix} \\ &= \begin{pmatrix} \displaystyle\sum_{i = 1}^k q(i)p^n(i, 1) & \displaystyle\sum_{i = 1}^k q(i)p^n(i, 2) & \cdots & \displaystyle\sum_{i = 1}^k q(i)p^n(i, k) \end{pmatrix} \\ &= \mathbf q^n \end{aligned} $$记$~\mathbf q\mathbf P^n = \pi~$，若$~\pi \mathbf P = \pi~$，则称$~\pi~$为平稳概率向量，其中的各元素组成平稳概率分布。若$~\displaystyle\lim_{n\to\infty}\mathbb P(X_n = i) = \bar\pi(i)~$，则称此$~\bar\pi(i)~$组成的是极限分布。对于不可约、非周期的马氏链，其平稳分布和极限分布是相等的。设$~(k\times k)~$转移矩阵$~\mathbf P~$不可约，则$~\pi\mathbf P = \pi~$存在唯一解，其中$~\displaystyle\sum_k \pi_k = 1,~\forall x, \pi_x > 0~$.1.6. 极限行为如果$~y~$是一个非常返态，则对$~\forall x~$，均有$$ \sum_{n = 1}^\infty p^n(x, y) < \infty $$从而$$ p^n(x, y)\to 0 $$这意味着我们只需要将注意力集中在常返态上，集中在马氏链只包含一个不可约常返集的情形。一个状态$~x~$的周期是$~I_x = \{ n\ge 1: p^n(x, x) > 0 \}~$的最大公约数，即$~\mathrm{ged}(I_x)~$.集合$~I_x~$是加法封闭集。若$~x~$的周期为$~1~$，则$~\exist n_0~$，使得$~n \ge n_0~$时，$~n\in I_x~$，即：$~x~$的周期为$~1~$时，$~I_x~$包含某个数值$~n_0~$之后的所有整数。若$~p(x, x) > 0~$，则状态$~x~$的周期为$~1~$.若$~\rho_{xy} > 0~$，且$~\rho_{yx} > 0~$，则状态$~x~$和$~y~$具有相同的周期。收敛定理：假设$~\mathbf P~$不可约、非周期且具有平稳分布$~\pi~$，则有$$ \lim_{n\to\infty}p^n(x, y) = \pi(y) $$假设$~\mathbf P~$不可约，且所有状态均常返，则存在平稳测度$~\mu(x)~$，使得$~\forall x,~\mu(x) > 0~$.渐进频率：假设$~\mathbf P~$不可约，且所有状态均常返，记$~N_n(y)~$为在时刻$~n~$之前访问$~y~$的总次数，则$$ \frac{N_n(y)}{n} \to \frac{1}{\mathbb E_y[T_y]},\qquad \text{a.s.} $$假设$~\mathbf P~$不可约，且具有平稳分布$~\pi~$，则$$ \pi(y) = \frac{1}{\mathbb E_y[T_y]} $$这说明状态$~y~$下的平稳分布对应的概率$~\pi(y)~$等于首次返回状态$~y~$期望步数的倒数；并且在状态期间$~n\to\infty~$时等于返回状态$~y~$的次数占整个状态步数$~n~$的比例。假设$~\mathbf P~$不可约，且具有平稳分布$~\pi~$，并且对$~\forall f: S\to \mathbb R~$，有$~\displaystyle\sum_x|f(x)|\pi(x) < \infty~$（绝对收敛），则$$ \frac{1}{n}\sum_{m = 1}^n f(X_m) \to \sum_{x\in S}\pi(x)f(x),\qquad \text{a.s.} $$其中$~f(X_m)~$是状态空间$~S~$下的一个函数，对应的是状态$~m~$下的取值。上式可以理解为：$~n~$各状态下随机变量的平均值，等于平稳状态下随机变量的期望值。大数定律的前提是随机变量独立同分布，并且随机变量的期望值存在。包含两类：弱大数定律：采样次数$~n~$越大，样本平均值接近随机变量真实期望值的可能性越来越大（即依概率收敛）。强大数定律：采样次数$~n~$越大，样本平均值几乎一定越来越接近随机变量真实期望值（即几乎确定收敛）。1.7. 一些例子双随机链：弱转移矩阵格列元素之和均为$~1~$，则称该矩阵$~\mathbf P~$是双随机的。即$$ \sum_{y}p(x, y) = 1 $$若$~\mathbf P~$是$~N~$状态马氏链的双随机转移概率，则均匀分布$~\pi(x) = 1 / N,~\forall x~$是其平稳分布。若$~\mathbf P~$是$~N~$状态马氏链的转移概率，且均匀分布$~\pi(x) = 1/ N,~\forall x~$是其平稳分布，则$~\mathbf P~$是双随机的。细致平衡条件：如果$~\pi(x)p(x, y) = \pi(y)p(y, x)~$，则称$~\pi~$满足细致平衡条件。该条件比之前的$~\pi\mathbf P = \pi~$更严格。$~\pi\mathbf P = \pi~$说明在所有的转移结束后，每个状态的概率与初始时的概率相等。而在细致平衡条件下，从$~x~$状态一步转移到$~y~$状态的概率，刚好等于$~y~$状态一步转移至$~x~$状态的概率。满足细致平衡条件的马氏链一定存在平稳分布，但是具有平稳分布的马氏链不一定满足细致平衡条件。1.8. 离出分布离出分布考察的是：对于至少存在两个不同常返类的马氏链，从给定的非常返态$~j~$开始，该马氏链最终进入某一特定常返态$~i~$的概率是多少。换句话说，离出分布考察的是从一个非常返态$~j~$最终被某一状态$~i~$吸收的概率。非常返态$~x~$最终被状态$~z~$吸收的概率为$~h(x) = \mathbb P(X_T = z\mid X_0 = x)~$，根据$~\text{C-K}~$方程，有$$ \begin{aligned} \mathbb P(X_T = z \mid X_0 = x) &= \sum_y\mathbb P(X_T = z \mid X_1 = y)\mathbb P(X_1 = y\mid X_0 = x) \\ &= \sum_y p(x, y)\mathbb P(X_T = z\mid X_1 = y) \\ h(x) &= \sum_y p(x, y)h(y) \end{aligned} $$离出时刻考察的是：对于至少存在两个不同常返类的马氏链，从给定的非常返态$~j~$开始，该马氏链最终进入常返态$~i~$的期望步数是多少。换句话说，离出时刻考察的是从一个非常返态$~j~$最终被吸收态吸收的期望时间。非常返态$~x~$最终被吸收的期望时间为$~g(x) = \mathbb E_x(T)~$，根据马氏性，对任意非常返态$~y~$，有$$ \begin{aligned} \mathbb E_x(T) &= 1 + \sum_y \mathbb P(X_1 = y \mid X_0 = x)\mathbb E_y(T) \\ &= 1 + \sum_y p(x, y)\mathbb E_y(T) \\ g(x) &= 1 + \sum_y p(x, y)g(y) \end{aligned} $$无限状态空间马氏链遇到的新问题是常返并不能保证平稳分布的存在。无限状态空间马氏链也称可数状态马氏链。在无限状态空间中，常返分为两类：正常返：$~\mathbb P_x(T_x < \infty) = 1,~ \mathbb E_x(T_x) < \infty~$.零常返：$~\mathbb P_x(T_x < \infty) = 1,~\mathbb E_x(T_x) = \infty~$.正常返的马氏链，一定可以找到其对应的平稳分布；若马氏链不存在平稳分布，则其可能为零常返或非常返。对于一个不可约链，下面的结论等价：某状态是常返的；存在一个平稳分布$~\pi~$；所有的状态都是正常返的。2. 连续时间马氏链离散时间马氏链：对任意状态$~i, j, i_0, \cdots, i_{n-1}~$，有$$ \mathbb P(X_{n+1} = j\mid X_n = i, X_{n-1} = i_{n-1}, \cdots, X_0 = i_0) = \mathbb P(X_{n+1} = j\mid X_n = i) $$连续时间马氏链：对于任意状态$~i, j, i_0, \cdots, i_{n}~$以及任意时间$~0 \le s_0 < s_1 <\cdots < s_n < s~$，有$$ \begin{aligned} \mathbb P(X_{t+s} = j\mid X_s = i, X_{s_n} = i_n, \cdots, X_{s_0} = i_0) &= \mathbb P(X_{t+s} = j\mid X_s = i) \\ &= \mathbb P(X_t = j\mid X_0 = i) \end{aligned} $$这里的$~s_i~$可以取任意的正实数。从时刻$~s~$的状态$~i~$到时刻$~(t + s)~$的状态$~j~$的概率，只依赖于时间间隔$~t~$.在离散时间马氏链中，时间和状态均是离散的。从状态$~i~$一步转移到状态$~j~$的概率记为$~p(i, j)~$，即$$ p(i, j) = \mathbb P(X_{t + 1} = j \mid X_t = i) $$在连续时间马氏链中，时间是连续的，状态是离散的。在$~h > 0~$时间段，从状态$~i~$转移到状态$~j~$的概率记为$~p_h(i, j)~$，即$$ p_h(i, j) = \mathbb P(X_{t + h} = j\mid X_t = i) = \mathbb P(X_h = j \mid X_0 = i) $$连续时间马氏链的$~\text{Chapman-Kolmogorov}~$方程为$$ \sum_k p_s(i, k)p_t(k, j) = p_{s + t}(i, j) $$当$~h \to 0~$时，我们引入转移速率$~q(i, j)~$，即$$ q(i, j) = \lim_{h\to 0}\frac{p_h(i, j)}{h} = \lim_{h\to 0}\frac{\mathbb P(X_{t + h} = j\mid X_t = i)}{h},\qquad i\neq j $$在连续时间马氏链中，我们除了要考虑某一时刻马氏链处于什么状态，还关心它在离开这个状态之前会停留多久。这个停留时间具有无记忆性，服从指数分布。转移速率$~q(i, j)~$有如下性质$~q(i, i) \le 0,~i = 1, 2, \cdots, n~$$~q(i, j) \ge 0,~i \neq j,~~i,j = 1, 2, \cdots, n~$$~\displaystyle\sum_j q(i, j) = 0,~i = 1, 2, \cdots, n~$对于$~N~$各状态的连续时间马氏链，其转移速率往往以矩阵形式表示$$ \mathbf Q = \begin{pmatrix} q(1, 1) & q(1, 2) & \cdots & q(1, N) \\ q(2, 1) & q(2, 2) & \cdots & q(2, N) \\ \vdots & \vdots & \ddots & \vdots \\ q(N, 1) & q(N, 2) & \cdots & q(N, N) \end{pmatrix} $$降$~[0, t + h]~$拆分成$~[0, h]~$和$~[h, t + h]~$，则$$ \begin{aligned} p_{t + h}(i, j) - p_t(i, j) &= \left[ \sum_k p_h(i, k)p_t(k, j) \right] - p_t(i, j) \\ &= \left[ \sum_{k \neq i}p_h(i, k)p_t(k, j) \right] + [p_h(i, i) - 1]p_t(i, j) \end{aligned} $$两端同时除以$~h~$，并令$~h\to 0~$，得到$$ \lim_{h\to 0}\frac{p_{t + h}(i, j) - p_t(i, j)}{h} = \lim_{h\to 0}\left[ \sum_{k \neq i}\frac{p_h(i, k)}{h}p_t(k, j) \right] + \left[ \lim_{h\to 0}\frac{p_h(i, i) - 1}{h} \right]p_t(i, j) $$即$$ \begin{aligned} p'_t(i, j) &= \sum_{k\neq i}q(i, k)p_t(k, j) + \left[ \lim_{h\to 0}\frac{p_h(i, i) - 1}{h} \right]p_t(i, j) \\ &= \sum_{k\neq i} q(i, k)p_t(k, j) - \left[ \lim_{h\to 0}\frac{1 - p_h(i, i)}{h} \right]p_t(i, j) \\ &= \sum_{k\neq i} q(i, k)p_t(k, j) - \lambda_i p_t(i, j) \\ &= \sum_k Q(i, k)p_t(k, j) \end{aligned} $$其中$$ Q(i, k) = \begin{cases} q(i, k), & i\neq k \\ -\lambda_i, & i = k \end{cases} $$这个式子可以用矩阵表示，即$$ \frac{\mathrm d \mathbf P_t}{\mathrm d t} = \mathbf Q\mathbf P_t $$该方程称为$~\text{Kolmogorov}~$向后方程。其中$~\mathbf Q~$称为转移速率矩阵。在这里我们将$~[0, t + h]~$拆分成了$~[0, h]~$和$~[h, t + h]~$.类似地，如果将$~[0, t + h]~$拆分成$~[0, t]~$和$~[t, t + h]~$，就得到向前方程。用矩阵可以简写为$$ \frac{\mathrm d \mathbf P_t}{\mathrm d t } = \mathbf P_t \mathbf Q $$2.1. 极限行为如果对任意状态$~i~$和$~j~$，都有可能从$~i~$经过有限步转移到$~j~$，则称马氏链是不可约的，即$$ \exist k_0 = i, k_1, \cdots, k_n = j,\quad\text{s.t.}\quad q(k_{m - 1}, k_m) > 0, 1 \le m \le n $$若$~X_t~$不可约，且$~t > 0~$，则$~p_t(i, j) > 0~$.随机停留时间满足指数分布，这意味着我们不再需要考虑时间的非周期性，因此连续时间马氏链的极限行为比离散时间马氏链更简单。在离散时间下，平稳分布是$~\pi\mathbf P = \pi~$的一个解，而在连续时间下，对$~\forall t > 0~$，都有$~\pi\mathbf P_t = \pi~$成立。但在实际操作中，$~\mathbf P_t~$可能难以计算。分布$~\pi~$是一个平稳分布，当且仅当$~\pi\mathbf Q = 0~$.如果一个连续时间马氏链$~X_t~$不可约，且具有平稳分布$~\pi~$，则$$ \lim_{t\to\infty}p_t(i, j) = \pi(j) $$相应的细致平衡条件为对$~\forall j \neq k~$，满足$~\pi(k)q(k, j) = \pi(j)q(j, k)~$. 若时间连续马氏链满足细致平衡条件，则$~\pi~$是一个平稳分布。我们可以根据此定理，用细致平衡条件来求得平稳分布。