当时间步长不小于 Python 时，如何计算 Kramers-Moyal 系数？

Question

我正在 python 中编写一个简单的一维 Ornstein-Uhlenbeck 模拟来解释来自视频源的一些实验数据，即跟踪帧序列中的布朗粒子。相机的采样频率是 fps=30 帧/秒，这意味着我的时间分辨率大约是 dt=0.03 秒。另一方面，我能够测量的最小时间尺度约为 tau=0.5 秒，比 dt 大一个数量级（这是一个合理的差异，能够计算可观测值，例如过程）。

现在，对于模拟数据，我使用 kramersmoyal 计算前两个 Kramers-Moyal 系数。只有当 dt 不比 tau 小至少两个数量级时，第二个系数才会出现问题，在这种情况下，噪声强度的估计变为凸函数而不是平坦函数。

所以手头的问题是：

是否有一种简单的方法来编写一些更正代码来解决这个问题？
有没有其他软件包可以解决这个问题？
这是一个实现问题还是这些系数是在无限小的时间步长（即 dt->0）的限制下定义的结果？

为了更具体，下面的代码显示了两个不同时间步长的模拟，dt，分别与 tau 相差两个和一个数量级。然后为这两种情况计算第二个系数，并估计过程的噪声强度。第一种情况的噪声强度看起来很平坦，但第二种情况则不然。

# === Ornstein-Uhlenbeck simulation ===

import numpy as np
import kramersmoyal as km
import matplotlib.pyplot as plt

np.random.seed(0)

n = 100_000       # total number of time steps
dt = [.003, .03]  # time steps two/one order of magnitud smaller than tau
tau = .5          # relaxation or characteristic time
sigma = 5         # noise intensity

scale_dW = sigma * np.sqrt(dt)    # scale parameter of the Wiener process
scale_x = sigma / np.sqrt(2/tau)  # scale parameter of the process under study

dW = np.random.normal(scale=scale_dW, size=(n, 2))  # Wiener process array
x = np.empty((n, 2))                                # x array
x[0] = np.random.normal(scale=scale_x, size=2)      # initial condition

# Euler-Maruyama method to solve the SDE
for i in range(n - 1):
    x[i + 1] = x[i] - 1/tau * x[i] * dt + dW[i]

# === Second Kramers-Moyal coefficient ===

fig, axs = plt.subplots(1, 2, sharey=True, figsize=(9,3.5))
bins = np.array([100])         # number of bins
powers = np.array([[1], [2]])  # first and second power

for i in range(len(dt)):
    # computation of the coefficients using `kramersmoyal`
    kmc, edges = km.km(x[:,i], bins, powers)

    # sigma parameter estimation
    sigma_est = np.sqrt(2 * kmc[1] / dt[i])

    # plotting
    axs[i].plot(edges[0], sigma_est, '.-', label=f'dt = {dt[i]}')
    axs[i].legend()
    axs[i].set_title(['FLAT', 'NOT FLAT'][i])
    axs[i].set_xlabel('x')
    axs[i].set_ylabel('estimated sigma')
plt.show()

Answer 1

tl;dr：将校正添加到扩散中

D₂(x) = (M₂(x) - (M₁(x)²)/2)/dt

或在您的符号中使用 kmc

sigma_est = np.sqrt(2 * (kmc[1] - 0.5*kmc[0]**2) / dt[i])

这实际上是一个已知的“问题”，它源于在求解 Kramers-Moyal 方程（或类似的 Fokker-Planck 方程）时表达 Kramers-Moyal 运算符的数学近似。您可以在 On the definition and handling of different drift and diffusion estimates、New Journal of Physics 10、083034、2008 中找到此问题的 first-ish 示例。

为了直接得到答案，对

给出的第二个 Kramers-Moyal 系数进行了 finite-time 修正

D₂(x) = (M₂(x) - (M₁(x)²)/2)/dt

这应该会立即纠正您观察到的凸性（或二次效应）。真是finite-time效果。

您可以在 Arbitrary-Order Finite-Time Corrections for the Kramers–Moyal Operator、Entropy、23[=38= 中找到所有 Kramers-Moyal 系数的全套校正](5), 517, 2021.

披露：我是 kramersmoyal python 图书馆和最后提到的出版物的作者之一。

当时间步长不小于 Python 时，如何计算 Kramers-Moyal 系数？

How to compute Kramers-Moyal coefficients when the time step is not to small in Python?

python

simulation

numpy

stochastic