模拟SDE的解

以最简单的随机微分方程为例，考虑随机游走
dxt=ξtdtdx_t = \xi_t dtdxt=ξtdt
其中ξt\xi_tξt是白噪声，假设⟨ξt1,ξt2⟩=2Dδ(t1−t2)\langle \xi_{t_1},\xi_{t_2} \rangle = 2D\delta(t_1-t_2)⟨ξt1,ξt2⟩=2Dδ(t1−t2)，其中D>0D>0D>0，δ\deltaδ是Dirac函数。假设随机过程xtx_txt的概率密度函数为ρ(t,x)\rho(t,x)ρ(t,x)，它满足Fokker-Planck方程：
∂ρ∂t=D∂2ρ∂x2\frac{\partial \rho}{\partial t} = D \frac{\partial^2 \rho}{\partial x^2}∂t∂ρ=D∂x2∂2ρ
这个是一维的扩散方程，DDD叫做扩散系数。这个方程的一个解是
ρ(t,x)=14πDtexp⁡(−x24Dt)\rho(t,x) = \frac{1}{\sqrt{4\pi Dt}}\exp \left(- \frac{x^2}{4Dt} \right)ρ(t,x)=4πDt1exp(−4Dtx2)
即xtx_txt是高斯过程，分布可以记为N(0,2Dt)N(0,2Dt)N(0,2Dt)。

现在对时间做离散化，记time step为τ\tauτ，简单回顾一下解数值ODE的Runge-Kutta方法 (RK)：对于动力系统x˙=f(t,x)∈Rm\dot{x}=f(t,x)∈\mathbb{R}^mx˙=f(t,x)∈Rm，定义RK的系数为A∈Rs×s,b∈Rs,c∈RsA∈\mathbb{R}^{s×s}, b∈\mathbb{R}^s,c∈\mathbb{R}^sA∈Rs×s,b∈Rs,c∈Rs
ki=hf(t+cih,x(t)+∑j=1sAijkj)x(t+h)=x(t)+∑j=1sbjkjk_i=hf(t+c_i h,x(t)+\sum_{j=1}^s A_{ij} k_j) \\ x(t+h)=x(t)+\sum_{j=1}^s b_j k_jki=hf(t+cih,x(t)+j=1∑sAijkj)x(t+h)=x(t)+j=1∑sbjkj
其中sss是RK方法的阶数，hhh是步长，根据上述递推关系，可以计算出x(t)x(t)x(t)的序列。

对于随机微分方程dxt=ξtdtdx_t = \xi_t dtdxt=ξtdt，考虑RK2：
A=[00120],b=[0,1]T,c=[0,12]T,h=τA = \left[ \begin{matrix} 0 & 0 \\ \frac{1}{2} & 0 \end{matrix} \right],b=[0,1]^T,c=[0,\frac{1}{2}]^T,h=\tauA=[02100],b=[0,1]T,c=[0,21]T,h=τ
则
k1=τξi,k2=τξi+12xi+1=xi+τξi+12k_1 = \tau \xi_i,k_2 = \tau \xi_{i+\frac{1}{2}} \\ x_{i+1} = x_i + \tau \xi_{i+\frac{1}{2}}k1=τξi,k2=τξi+21xi+1=xi+τξi+21
xi+1=xi+τξi+12x_{i+1} = x_{i} + \tau \xi_{i+\frac{1}{2}}xi+1=xi+τξi+21
其中⟨ξi,ξj⟩=2Dδij/τ\langle \xi_{i},\xi_{j} \rangle = 2D\delta_{ij}/\tau⟨ξi,ξj⟩=2Dδij/τ，其中δij\delta_{ij}δij为Kronecker符号，所有ξi\xi_iξi都是独立同分布的。因此可以将上式写成：
xi+1=xi+2Dτζζ∼N(0,1)x_{i+1} = x_{i} + \sqrt{2D\tau} \zeta\\\zeta \sim N(0,1)xi+1=xi+2Dτζζ∼N(0,1)

要解这个随机过程可以用随机模拟的方法：

第一步，从标准正态分布中抽一个样本ζ\zetaζ；
第二步，根据xi+1=xi+2Dτζx_{i+1} = x_{i} + \sqrt{2D\tau} \zetaxi+1=xi+2Dτζ更新随机过程的值；

重复这两步，直到达到设置的时间上限。

由此我们可以总结出模拟一个随机微分方程的一般性做法：

用数值ODE的方法导出SDE解的递推关系
生成正态分布的随机数，根据递推关系模拟这个SDE

例1下面的SDE表示几何随机游走
dxt=xtξtdtdx_t =x_t \xi_t dtdxt=xtξtdt
其中ξt\xi_tξt是白噪声，假设⟨ξt1,ξt2⟩=2Dδ(t1−t2)\langle \xi_{t_1},\xi_{t_2} \rangle = 2D\delta(t_1-t_2)⟨ξt1,ξt2⟩=2Dδ(t1−t2)，D>0D>0D>0。用这个例子介绍两种离散化的规则：

Ito规则：xi+1=xi+2Dτxiζx_{i+1}=x_i + \sqrt{2D\tau}x_i \zetaxi+1=xi+2Dτxiζ
Stratonovich规则：xi+1=xi+2Dτxi+xi+12ζx_{i+1}=x_i + \sqrt{2D\tau}\frac{x_i + x_{i+1}}{2} \zetaxi+1=xi+2Dτ2xi+xi+1ζ

在这个例子中，Ito规则导出的结果是Exi+1=Exi=x0Ex_{i+1}=Ex_i = x_0Exi+1=Exi=x0，这个比较明显。但Stratonovich规则导出的是
xi+1=xi1+2Dτζ21−2Dτζ2=xi(1+2Dτζ2)(1+2Dτζ2+2Dτζ24+o(ζ2))=xi(1+2Dτζ+Dτζ2+o(ζ2))→xi(1+Dτ+2Dτζ)x_{i+1} = x_i \frac{1+ \sqrt{2D\tau} \frac{\zeta}{2}}{1- \sqrt{2D\tau} \frac{\zeta}{2}} = x_i(1+ \sqrt{2D\tau} \frac{\zeta}{2})(1+\sqrt{2D\tau} \frac{\zeta}{2}+2D\tau \frac{\zeta^2}{4}+o(\zeta^2)) \\ = x_i(1+\sqrt{2D\tau}\zeta + D\tau \zeta^2 + o(\zeta^2)) \to x_i(1+D\tau+\sqrt{2D\tau}\zeta )xi+1=xi1−2Dτ2ζ1+2Dτ2ζ=xi(1+2Dτ2ζ)(1+2Dτ2ζ+2Dτ4ζ2+o(ζ2))=xi(1+2Dτζ+Dτζ2+o(ζ2))→xi(1+Dτ+2Dτζ)
因此
Exi+1=Exi(1+Dτ)Ex_{i+1} = Ex_i (1+D \tau)Exi+1=Exi(1+Dτ)
因此这两个规则模拟出来的结果会有明显差别。

用蒙特卡罗方法估计SDE解的自相关函数

随着时间流逝，xtx_txt信息减少的规律可以用自相关函数来表示，
KXX(τ)=⟨(xt+τ−xt)2⟩=∫tt+τdt1∫tt+τdt2⟨ξt1ξt2⟩=2DτK_{XX}(\tau) = \langle (x_{t+\tau}-x_t)^2 \rangle = \int_t^{t+\tau} dt_1 \int_t^{t+\tau} dt_2 \langle \xi_{t_1} \xi_{t_2}\rangle =2D \tauKXX(τ)=⟨(xt+τ−xt)2⟩=∫tt+τdt1∫tt+τdt2⟨ξt1ξt2⟩=2Dτ
随机游走的自相关函数显然只与两个时点的差有关，这种随机过程叫平稳随机过程。如果SDE的解比较复杂，甚至根本写不出解析解，很难直接计算自相关函数，我们可以用Monte Carlo方法估计自相关函数：
KXX(τ)=Extxt+τ≈1T∑i=1Txixi+τK_{XX}(\tau) = Ex_tx_{t+\tau} \approx \frac{1}{T} \sum_{i=1}^T x_i x_{i+\tau}KXX(τ)=Extxt+τ≈T1i=1∑Txixi+τ
其中xi,xi+τx_i,x_{i+\tau}xi,xi+τ可以用随机模拟得到。

例2 一个简单的单位根过程，假设扩散系数D=12D=\frac{1}{2}D=21
dxt=(−xt+ξt)dtdx_t = (-x_t + \xi_t )dtdxt=(−xt+ξt)dt
参考SDE的第一讲，这个方程定义了一个Ornstein-Uhlenbeck过程，它的解是
xt=∫−∞te−(t−s)ξsdsx_t = \int_{-\infty}^t e^{-(t-s)}\xi_s dsxt=∫−∞te−(t−s)ξsds
根据这个解计算自相关函数，第一行到第二行用到的性质是Ito Isometry：
KXX(τ)=Extxt+τ=E∫−∞te−(t−s)ξsds∫−∞t+τe−(t+τ−r)ξrdr=∫−∞t∫−∞t+τe−(t−s)e−(t+τ−r)Eξsξrdsdr=∫−∞t∫−∞t+τe−(t−s)e−(t+τ−r)δ(r−s)dsdr=12exp⁡(−τ),τ≥0K_{XX}(\tau) = Ex_{t}x_{t+\tau} = E \int_{-\infty}^t e^{-(t-s)}\xi_s ds \int_{-\infty}^{t+\tau} e^{-(t+\tau-r)}\xi_r dr \\ =\int_{-\infty}^t \int_{-\infty}^{t+\tau} e^{-(t-s)}e^{-(t+\tau-r)} E\xi_s\xi_r dsdr \\ =\int_{-\infty}^t \int_{-\infty}^{t+\tau} e^{-(t-s)}e^{-(t+\tau-r)} \delta(r-s) dsdr = \frac{1}{2} \exp(-\tau),\tau\ge0 KXX(τ)=Extxt+τ=E∫−∞te−(t−s)ξsds∫−∞t+τe−(t+τ−r)ξrdr=∫−∞t∫−∞t+τe−(t−s)e−(t+τ−r)Eξsξrdsdr=∫−∞t∫−∞t+τe−(t−s)e−(t+τ−r)δ(r−s)dsdr=21exp(−τ),τ≥0

如果我们不知道Ito Isometry，要计算自相关函数的解析解就会比较麻烦，这时就可以用Monte Carlo方法来估计自相关函数。先用RK2得到SDE的数值算法：
xi+1=e−τxi+τζ≈(1−τ)xi+τζx_{i+1}=e^{-\tau}x_i + \sqrt{\tau} \zeta \approx (1-\tau)x_i + \sqrt{\tau} \zeta xi+1=e−τxi+τζ≈(1−τ)xi+τζ
根据这个递推关系模拟SDE的解，然后根据
KXX(τ)=Extxt+τ≈1T∑i=1Txixi+τK_{XX}(\tau) = Ex_tx_{t+\tau} \approx \frac{1}{T} \sum_{i=1}^T x_i x_{i+\tau}KXX(τ)=Extxt+τ≈T1i=1∑Txixi+τ
估计自相关函数。

UA MATH575B 数值分析下VI 统计物理的随机模拟方法1相关推荐

UA MATH575B 数值分析下VI 统计物理的随机模拟方法2
UA MATH575B 数值分析下VI 统计物理的随机模拟方法2 扩散过程的停时问题理论解模拟解扩散过程的停时问题简单一点,考虑一维的随机扩散问题,假设扩散系数D=12D=\frac{1}{2 ...
UA MATH575B 数值分析下计算统计物理例题2
UA MATH575B 数值分析下计算统计物理例题2 理论解法 C-K方程法特征值法(近似解) 模拟解法 Rejection Sampling Importance Sampling 一个位于原点 ...
UA MATH575B 数值分析下计算统计物理例题1
UA MATH575B 数值分析下计算统计物理例题1 统计物理方法的解析解 Markov链理论解数值解 Monte Carlo模拟. 一道有趣的统计物理的题目.下面这个简单的迷宫中,一只老鼠一开 ...
UA MATH575B 数值分析下统计物理的随机模拟方法5
UA MATH575B 数值分析下统计物理的随机模拟方法5 Ising Model Gibbs Sampling Glauber Dynamics 这一讲介绍Ising Model,它是MCMC与G ...
UA MATH575B 数值分析下统计物理的随机模拟方法4
UA MATH575B 数值分析下统计物理的随机模拟方法4 这一讲介绍MCMC方法,这个方法最早出现在Metropolis在1953年发在J Chem Phys上的Equation of state ...
UA MATH575B 数值分析下IV 带约束的优化
UA MATH575B 数值分析下IV 带约束的优化问题带等式约束的优化问题带不等式约束的优化问题同时带等式约束与不等式约束的优化问题今天不想敲公式,就不写理论了,反正方法也就是前面的Newt ...
UA MATH575B 数值分析下III 图像恢复
UA MATH575B 数值分析下III 图像恢复图像去噪的优化模型算法实现图像去噪的优化模型假设VVV代表原图,观测到的图像是WWW,WWW是原图与噪声的叠加W=V+ξW=V+\xiW=V+ ...
UA MATH575B 数值分析下I 梯度下降
UA MATH575B 数值分析下I 梯度下降梯度下降最速下降下降算法实现举例梯度下降的应用对于凸优化问题min⁡xf(x)\min_x f(x)minxf(x),最主流的数值计算方法是下 ...
UA MATH575B 数值分析下II 牛顿算法
UA MATH575B 数值分析下II 牛顿算法 Pure Newton算法 Damped Newton算法 Levenberg-Marquardt算法 Quasi-Newton算法割线法 Broy ...

UA MATH575B 数值分析下VI 统计物理的随机模拟方法1

UA MATH575B 数值分析下VI 统计物理的随机模拟方法1

模拟SDE的解

用蒙特卡罗方法估计SDE解的自相关函数

UA MATH575B 数值分析下VI 统计物理的随机模拟方法1相关推荐

最新文章

热门文章