Jiming's blog - 随机过程
type: channel
link: https://www.jiming.site/tag/stochastic-process/
items: 过阻尼 Langevin 方程的路径积分形式
items: 随机过程的修正与不可分辨性
items: Python generator 实现 Gillespie 算法
items: 随机分析、Langevin 方程与 Fokker-Planck 方程
items: 关于随机热力学的简短介绍
items: 马尔可夫过程

过阻尼 Langevin 方程的路径积分形式
type: item
date: 2025-02-19T09:22:00+08:00
description: 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 \r…

随机过程的修正与不可分辨性
type: item
date: 2024-02-08T05:47:00+08:00
description: 考虑定义在同一个概率空间 $(\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$.如果我们进一步要求随机轨迹几乎处处右连续，那么“修正”可以推出“不可分辨”。

Python generator 实现 Gillespie 算法
type: item
date: 2023-10-03T11:16:00+08:00
description: 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 结构，它可以当作普通的可迭代对象使用，比如用在循…

随机分析、Langevin 方程与 Fokker-Planck 方程
type: item
date: 2023-08-21T11:00:00+08:00
description: 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_…

关于随机热力学的简短介绍
type: item
date: 2023-04-11T16:31:00+08:00
description: 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} …

马尔可夫过程
type: item
date: 2023-04-10T20:41:00+08:00
description: 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~$的最…