我正在撰寫一個簡單的一維 Ornstein-Uhlenbeck 模擬python來解釋來自視頻源的一些實驗資料，即跟蹤幀序列中的布朗粒子。相機的采樣頻率是fps=30每秒幀數，這意味著我的時間解析度dt=0.03大約是秒。另一方面，我能夠測量的最小時間尺度大約是tau=0.5秒，比dt（一個合理的差異能夠計算可觀察量，例如程序的自相關函式）大一個數量級。
現在，對于模擬資料，我正在計算前兩個 Kramers-Moyal 系數，使用kramersmoyal. 僅當第二個系數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()
uj5u.com熱心網友回復：
tl; dr：將校正添加到擴散中
D?(x) = (M?(x) - (M?(x)2)/2)/dt
或kmc在您的符號中使用
sigma_est = np.sqrt(2 * (kmc[1] - 0.5*kmc[0]**2) / dt[i])
這實際上是一個已知的“問題”，它源于在求解 Kramers-Moyal 方程（或類似的 Fokker-Planck 方程）時表達 Kramers-Moyal 算子的數學近似。您可以在關于不同漂移和擴散估計的定義和處理，New Journal of Physics 10 , 083034, 2008中找到這個問題的第一個示例。
為了直接得到答案，對第二個 Kramers-Moyal 系數進行有限時間校正，由下式給出
D?(x) = (M?(x) - (M?(x)2)/2)/dt
這應該立即糾正您觀察到的凸度（或二次效應）。這確實是一種有限時間的效果。
您可以在Arbitrary-Order Finite-Time Corrections for the Kramers–Moyal Operator , Entropy , 23 (5), 517, 2021中找到所有 Kramers-Moyal 系數的全套校正。
披露：我是kramersmoyalpython 庫和最后提到的出版物的作者之一。
轉載請註明出處，本文鏈接：https://www.uj5u.com/qukuanlian/479493.html