用Python进行随机积分

Question

我希望能够对包含白噪声的积分进行数值求解。

在数学上，白噪声可以用变量X(t)来描述，它是一个具有时间平均值Avg[X(t)] = 0和相关函数Avg[X(t), X(t')] = delta_distribution(t-t')的随机变量。

一个简单的例子是计算从t=0到t=1上X(t)的积分。虽然平均值为零，但我需要的是这个积分的不同实现方式。

问题在于，使用numpy.integrate.quad()并不能解决这个问题。

是否有任何适用于随机积分的python包？

- physicsGuy

---

如果这篇文章对您没有帮助，那么在stats.stackexchange上提问可能会更好。 - rocksportrocker

你能详细解释一下你的问题吗？为什么quad()不够用？我知道它通常假设你要积分的函数具有某种平滑性，但似乎结果的准确性会严重依赖于函数相对于噪声的幅度。为什么你不能使用np.trapz()计算一个有限的积分呢？ - user545424

我相信我当时想要问的问题在rocksportrocker发布的链接中得到了澄清。然而，我现在无法删除或编辑我的问题... - physicsGuy

Answer 1

这是一个关于数值SDE方法的良好起点：http://math.gmu.edu/~tsauer/pre/sde.pdf。

下面是我去年为课程项目编写的一种简单的numpy求解器，用于随机微分方程dX_t = a(t,X_t)dt + b(t,X_t)dW_t。它基于常规微分方程的向前欧拉方法，在实践中解决SDEs时被广泛使用。

def euler_maruyama(a,b,x0,t):
    N = len(t)
    x = np.zeros((N,len(x0)))
    x[0] = x0
    for i in range(N-1):
        dt = t[i+1]-t[i]
        dWt = np.random.normal(0,dt)
        x[i+1] = x[i] + a(t[i],x[i])*dt + b(t[i],x[i])*dWt
    return x

基本上，在每个时间步长，函数的确定性部分使用向前欧拉法进行积分，并且随机部分通过生成一个均值为0，方差为$dt$的正态随机变量$dWt$来进行积分。

我们之所以这样生成$dWt$，是基于布朗运动的定义。特别地，如果$W$是布朗运动，则$(W_t-W_s)$服从均值为0，方差为$t-s$的正态分布。因此，$dWt$是对小时间间隔内$W$的变化的离散化。

这是上述函数的文档字符串：

Parameters
----------
a : callable a(t,X_t),
    t is scalar time and X_t is vector position
b : callable b(t,X_t),
    where t is scalar time and X_t is vector position
x0 : ndarray
    the initial position
t : ndarray
    list of times at which to evaluate trajectory

Returns
-------
x : ndarray
    positions of trajectory at each time in t

- tch