23、随机过程与伊藤随机微积分入门

随机过程与伊藤随机微积分入门

1. 化学反应用 Gillespie 方法模拟

在化学动力学模拟中，Gillespie 方法是一种有效的工具。它引入了几个关键参数，如 (a_1 = k_1ax)，(a_2 = k_2xy)，(a_3 = k_3y)，以及 (a_0 = a_1 + a_2 + a_3)。这些参数类似于概率 (p_1)，(p_2)，(p_3)。同时，还引入了一个从归一化均匀分布中选取的随机数 (r_2)。

根据 (r_2a_0) 的取值范围，会发生不同的反应：
- 若 (0 < r_2a_0 \leq a_1)，则第一个反应发生，(x(t + \Delta t) = x(t) + 1)，(y(t + \Delta t) = y(t))。
- 若 (a_1 < r_2a_0 \leq a_1 + a_2)，则第二个反应发生，(x(t + \Delta t) = x(t) - 1)，(y(t + \Delta t) = y(t) + 1)。
- 若 (a_1 + a_2 < r_2a_0 \leq a_0)，则第三个反应发生，(x(t + \Delta t) = x(t))，(y(t + \Delta t) = y(t) - 1)。

与其他方法的一个显著区别在于时间步长，它不再是常数，而是随时间变化，(\Delta t = \ln(1/r_1)/a_0)，其中 (r_1) 是从归一化均匀分布中选取的随机变量。

以下是使用 Gillespie 技术模拟化学反应的 MATLAB 代码示例：

% 这里需要根据具体的 k1, k2, k3, a, x, y 等参数进行初始化
% 然后按照 Gillespie 方法的规则进行模拟
% 最后可以绘制 x(t) vs y(t) 的图像

还可以计算特定时间下 (x) 和 (y) 分子存在的概率密度函数。

2. 泊松过程

泊松随机过程是一种计数过程，用于计算随着时间增加，特定事件的发生次数。对于每个时间 (t)，有一个离散随机变量 (N(t))，表示在区间 ([0, t]) 内发生的事件数量，其可能取值为 ({0, 1, 2, \ldots})。

泊松过程的数学表达式为 (N(t) = \sum_{n=0}^{\infty} H(t - T[n]))，其中 (T[n]) 是第 (n) 次到达的时间。确定 (T[n]) 的值需要三个物理假设：
1. (N(0) = 0)。
2. (N(t)) 具有独立和平稳增量。平稳意味着对于任意两个相等的时间间隔 (\Delta t_1) 和 (\Delta t_2)，在 (\Delta t_1) 内发生 (n) 个事件的概率等于在 (\Delta t_2) 内发生 (n) 个事件的概率；独立意味着对于任意时间间隔 ((t, t + \Delta t))，在该间隔内发生 (n) 个事件的概率与之前发生的事件数量和方式无关。
3. (P[N(t + \Delta t) - N(t) = k ] = \begin{cases} 1 - \lambda\Delta t, & k = 0 \ \lambda\Delta t, & k = 1 \ 0, & k > 1 \end{cases})，其中 (\lambda) 是单位时间内事件的期望发生次数。

通过一系列推导，可以得到泊松过程的概率质量函数 (P_k(t) = \exp(-\lambda t)(\lambda t)^k / k!)，(k = 0, 1, 2, \ldots)。

在泊松过程的实现中，到达时间 (t_n) 是一个重要的随机变量。相邻两次事件发生的时间间隔 (Z_i = t_i - t_{i - 1}) 是独立同分布的，且 (P(Z_n \leq x) = 1 - e^{-\lambda x})，(x \geq 0)，(n = 1, 2, 3, \ldots)。

示例：随机电报信号

可以利用泊松过程中时间间隔的独立性和同分布性来生成随机电报信号 (X(t) = (-1)^{N(t)})。以下是实现该信号的 MATLAB 代码：

clear
N = 100; % number of switches in realization
lambda = 0.15; % switching rate
X = [ ];
% generate N uniformly distributed random variables
S = rand(1,N);
% transform S into an exponential random variable
T = - log(S)/lambda;
V = cumsum(T); % compute switching times
t = [0.01:0.01:100]; % create time array
icount = 1; amplitude = -1; % initialize X(t)
for k = 1:10000
    if ( t(k) >= V(icount) ) % at each switching point
        icount = icount + 1;
        amplitude = - amplitude; % switch sign
    end
    X(k) = amplitude; % generate X(t)
end
plot(t,X) % plot results
xlabel('\it t','FontSize',25);
ylabel('\it X(t)/a','FontSize',25);
axis([0 max(t) -1.1 1.1])

示例：工作任务完成概率

假设一个工人平均每小时组装一个小部件，将该任务建模为泊松过程。那么在八小时轮班内组装 12 个或更多小部件的概率可以通过泊松概率质量函数计算：
(P[N(t) \geq 12] = e^{-8} \sum_{n = 12}^{\infty} \frac{8^n}{n!} = 0.1119)

也可以使用 MATLAB 代码来模拟：

t_uniform = rand(1,12);
T = - log(1 - t_uniform);
total_time = sum(T);
% 多次执行上述代码，统计 total_time <= 8 的次数 icount
% 概率等于 icount / N

3. 随机微分方程

考虑一阶常微分方程 (\frac{dx}{dt} = a(t, x))，(x(0) = x_0)，其解为 (x(t) = x(0) + \int_{0}^{t} a[\eta, x(\eta)] d\eta)。

类似地，随机微分方程 (dX(t) = a[t, X(t)] dt)，(X(0) = X_0)，其中 (dX(t) = X(t + dt) - X(t))。当引入随机强迫项时，方程变为 (dX(t) = a[t, X(t)] dt + b[t, X(t)] dB(t))，(X(0) = X_0)，这里 (dB(t) = B(t + dt) - B(t))，(B(t)) 表示布朗运动，(a[t, X(t)]) 和 (b[t, X(t)]) 是确定性函数。这种由该方程支配的随机过程称为伊藤过程。

伊藤过程的解可以形式上写为 (X(t) = X_0 + \int_{0}^{t} a[\eta, X(\eta)] d\eta + \int_{0}^{t} b[\eta, X(\eta)] dB(\eta))，其中第一个积分是常规的黎曼积分，第二个积分是伊藤随机积分。

4. 随机微分方程的应用示例

LR 电路

LR 电路的数学模型为 (L\frac{dI}{dt} + RI = E(t))，其中 (I(t)) 是电流，(L) 是电感，(R) 是电阻，(E(t)) 是平均电动势。使用积分因子求解该一阶常微分方程，其解为 (I(t) = I(0) \exp(-\frac{Rt}{L}) + \frac{1}{L} \exp(-\frac{Rt}{L}) \int_{0}^{t} F(\tau) \exp(\frac{R\tau}{L}) d\tau)。

如果电动势是随机的，那么电流也是随机的。可以计算电流的均值和方差：
- 均值：若 (E[F(t)] = 0)，则 (E[I(t)] = I(0) \exp(-\frac{Rt}{L}))。
- 方差：(\sigma_{X(t)}^2 = E[I^2(t)] - {E[I(t)]}^2)，通过一系列推导和积分计算，最终得到不同情况下的方差表达式。

例如，当采用特定的自相关函数 (E[F(\tau)F(\tau’)] = 2D\delta(\tau - \tau’)) 时，方差为 (\sigma_{X(t)}^2 = \frac{DL}{R} (1 - e^{-2Rt/L}))。

阻尼谐波运动

阻尼谐波振荡器的方程为 (y’’ + 2\xi\omega_0y’ + \omega_0^2y = F(t))，其中 (0 \leq \xi < 1)，(y) 是位移，(t) 是时间，(\omega_0^2 = \frac{k}{m})，(2\xi\omega_0 = \frac{\beta}{m})，(m) 是振荡器的质量，(k) 是线性弹簧常数，(\beta) 是粘性阻尼器的常数。

其解为 (y(t) = y(0)e^{-\xi\omega_0t} [\cos(\omega_1t) + \frac{\xi\omega_0}{\omega_1} \sin(\omega_1t)] + \frac{y’(0)}{\omega_1} e^{-\xi\omega t} \sin(\omega_1t) + \int_{0}^{t} h(t - \tau)F(\tau) d\tau)，其中 (\omega_1 = \omega_0 \sqrt{1 - \xi^2})，(h(t) = \frac{e^{-\xi\omega_0t}}{\omega_1} \sin(\omega_1t)H(t))。

同样可以计算均值和方差，在 (E[F(t)] = 0) 的情况下，均值只取决于初始条件。方差的计算较为复杂，通过一系列积分和推导得到不同情况下的表达式。

5. 项目实践

低通滤波器与随机输入

考虑初始值问题 (y’ + y = f(t))，(y(0) = y_0)，其解为 (y(t) = y_0e^{-t} + e^{-t} \int_{0}^{t} e^{\tau}f(\tau) d\tau)。该微分方程与 RC 电路的方程相同，具有过滤高频干扰的特性。

操作步骤如下：
1. 使用 MATLAB 内置函数 randn 生成长度为 (N) 的平稳白噪声激励。
2. 根据生成的高斯随机强迫，编写 MATLAB 代码计算 (y(t))。
3. 生成多个 (y(t)) 的实现，使用 MATLAB 的 mean 和 var 函数计算均值和方差。
4. 计算 (y(t)) 在不同时间的概率密度函数。
5. 修改代码计算自协方差。

随机振动的首次通过问题

对于一个简单的轻微阻尼谐波振荡器 (y’’ + 2\zeta\omega_0y’ + \omega_0^2y = f(t))，(0 < \zeta \ll 1)，当受到白噪声强迫时，需要计算其振幅超过特定幅度的概率。

操作步骤如下：
1. 使用 MATLAB 命令 randn 生成长度为 (N) 的平稳白噪声激励。
2. 根据给定的高斯随机强迫，编写 MATLAB 代码计算 (y(t))。
3. 计算多个实现的 (y(t))，并存储在 (y(n,m)) 中，使用 MATLAB 命令计算均值和方差。
4. 存储 (y(n)) 首次超过特定幅度 (b > 0) 的时间 (T(m))，并使用 histc 估计概率密度函数。

随机强迫产生的波动

考虑一维波动方程 (\frac{\partial^2 u}{\partial t^2} - \frac{\partial^2 u}{\partial x^2} = \cos(\omega t)\delta[x - X(t)])，其解为 (u(x, t) = \frac{1}{2} \int_{0}^{t} H[t - \tau - |X(\tau) - x|] \cos(\omega\tau) d\tau)。

总结

本文介绍了随机过程和伊藤随机微积分的相关知识，包括化学反应用 Gillespie 方法模拟、泊松过程、随机微分方程及其应用示例，还给出了几个项目实践的操作步骤和代码示例。通过这些内容，可以更好地理解和处理随机现象在工程和科学中的应用。

流程图

graph TD;
    A[开始] --> B[化学反应用 Gillespie 方法模拟];
    B --> C[泊松过程];
    C --> D[随机微分方程];
    D --> E[随机微分方程的应用示例];
    E --> F[项目实践];
    F --> G[结束];

表格

项目	描述
化学反应用 Gillespie 方法模拟	引入参数，时间步长可变，可模拟化学反应
泊松过程	计数过程，有三个物理假设，可计算概率质量函数
随机微分方程	包括常规和引入随机强迫项的情况，解涉及伊藤随机积分
随机微分方程的应用示例	LR 电路和阻尼谐波运动，可计算均值和方差
项目实践	低通滤波器、随机振动首次通过问题、随机强迫产生的波动

随机过程与伊藤随机微积分入门（续）

6. 伊藤随机微积分的深入理解

伊藤随机微积分是处理随机微分方程的重要工具。在前面提到的随机微分方程 (dX(t) = a[t, X(t)] dt + b[t, X(t)] dB(t)) 中，伊藤随机积分 (\int_{0}^{t} b[\eta, X(\eta)] dB(\eta)) 与常规的黎曼积分有很大不同。

由于布朗运动 (B(t)) 处处不可微，伊藤随机积分不能用传统的微积分方法来定义。它的定义基于极限和期望的概念，是一种均方收敛意义下的积分。

伊藤引理是伊藤随机微积分中的一个重要定理。假设 (Y(t)=f(t,X(t)))，其中 (X(t)) 是满足上述随机微分方程的伊藤过程，那么 (Y(t)) 的随机微分 (dY(t)) 可以表示为：

(dY(t)=\left(\frac{\partial f}{\partial t}+a[t,X(t)]\frac{\partial f}{\partial X}+\frac{1}{2}b^{2}[t,X(t)]\frac{\partial^{2} f}{\partial X^{2}}\right)dt + b[t,X(t)]\frac{\partial f}{\partial X}dB(t))

伊藤引理为我们处理复合随机过程提供了有力的工具，使得我们能够计算更复杂的随机过程的变化率。

7. 随机过程的统计特性

随机过程的统计特性对于理解和分析随机现象非常重要。除了前面提到的均值和方差，还有协方差、相关函数等。

协方差

对于两个随机过程 (X(t)) 和 (Y(t))，它们的协方差定义为：

(Cov[X(t),Y(s)] = E[(X(t)-E[X(t)])(Y(s)-E[Y(s)])])

协方差衡量了两个随机过程在不同时间点上的线性相关性。如果 (Cov[X(t),Y(s)] = 0)，则称 (X(t)) 和 (Y(s)) 是不相关的。

相关函数

相关函数描述了随机过程在不同时间点上的相关性。自相关函数 (R_X(t,s)=E[X(t)X(s)]) 描述了同一个随机过程在不同时间点上的相关性，而互相关函数 (R_{XY}(t,s)=E[X(t)Y(s)]) 描述了两个不同随机过程在不同时间点上的相关性。

8. 项目实践拓展

低通滤波器与随机输入的拓展

在前面的低通滤波器项目中，我们可以进一步研究不同参数对滤波器性能的影响。例如，改变时间间隔 (\Delta t) 或噪声的强度，观察输出信号的均值、方差和概率密度函数的变化。

操作步骤如下：
1. 修改 MATLAB 代码中的参数，如时间间隔 (\Delta t) 和噪声强度。
2. 重新运行代码，计算不同参数下 (y(t)) 的均值、方差和概率密度函数。
3. 绘制不同参数下的均值、方差和概率密度函数曲线，进行对比分析。

随机振动的首次通过问题的拓展

在随机振动首次通过问题中，我们可以考虑不同的阻尼比 (\zeta) 和幅度 (b) 的组合，观察首次通过时间的概率密度函数的变化。

操作步骤如下：
1. 在 MATLAB 代码中设置不同的阻尼比 (\zeta) 和幅度 (b)。
2. 多次运行代码，计算不同组合下 (y(n)) 首次超过幅度 (b) 的时间 (T(m))。
3. 使用 histc 函数估计不同组合下的概率密度函数，并绘制对比图。

随机强迫产生的波动的拓展

对于随机强迫产生的波动问题，我们可以研究不同的随机过程 (X(t)) 对波动传播的影响。

操作步骤如下：
1. 定义不同形式的随机过程 (X(t))，如高斯随机过程、泊松随机过程等。
2. 修改 MATLAB 代码，将不同的随机过程 (X(t)) 代入波动方程的解 (u(x, t) = \frac{1}{2} \int_{0}^{t} H[t - \tau - |X(\tau) - x|] \cos(\omega\tau) d\tau) 中。
3. 计算不同随机过程下的波动解 (u(x, t))，并绘制波形图进行对比分析。

9. 随机过程在实际中的应用

随机过程在许多领域都有广泛的应用，如通信工程、金融、生物医学等。

通信工程

在通信工程中，随机过程用于描述信号的噪声和干扰。例如，高斯白噪声是一种常见的随机过程，用于模拟通信信道中的噪声。通过对随机过程的分析，可以设计更有效的通信系统，提高信号的传输质量。

金融

在金融领域，随机过程用于描述股票价格、利率等金融变量的变化。例如，几何布朗运动常用于模拟股票价格的变化，通过对随机过程的建模和分析，可以进行风险评估和投资决策。

生物医学

在生物医学中，随机过程用于描述生物系统中的随机现象，如神经元的放电、生物分子的扩散等。通过对随机过程的研究，可以更好地理解生物系统的行为和机制。

流程图

graph TD;
    A[伊藤随机微积分深入理解] --> B[随机过程的统计特性];
    B --> C[项目实践拓展];
    C --> D[随机过程在实际中的应用];

表格

项目	描述
伊藤随机微积分深入理解	介绍伊藤随机积分和伊藤引理，用于处理复合随机过程
随机过程的统计特性	包括协方差和相关函数，描述随机过程的相关性
项目实践拓展	对低通滤波器、随机振动首次通过问题和随机强迫产生的波动项目进行拓展研究
随机过程在实际中的应用	涉及通信工程、金融和生物医学等领域

总结

本文进一步深入介绍了随机过程和伊藤随机微积分的相关知识，包括伊藤随机微积分的深入理解、随机过程的统计特性、项目实践的拓展以及随机过程在实际中的应用。通过这些内容，我们可以更全面地掌握随机过程和伊藤随机微积分的理论和应用，为解决实际问题提供有力的工具。