我编写了一个非常简单的卡尔曼滤波器,它对时间序列进行操作(包括数据间隙)。它工作得很好,但我碰巧有一个数据数据立方体(例如,形状Nt,Ny,Nx
的数组),并且我想应用我的时态卡尔曼对数据立方体中的每个像素进行过滤。我已经完成了显而易见的事情(循环最后两个维度),但这需要相当长的时间。
最终,我总是不得不从各个“像素”中提取数据,并构建相关的矩阵/向量,因此这个过程非常慢(请注意,每个单独的时间序列中的间隙是不同的,通常,所以是将状态与观测值联系起来的 H 矩阵)。我不熟悉 cython,它可能会有所帮助(只是我不熟悉它)。
我只是想知道是否巧妙地重新表述问题或巧妙的数据结构是否可以更有效地进行时间过滤。我宁愿只使用 numpy/scipy,而不是 OpenCV,否则依赖额外的包会很麻烦。
最佳答案
我制作了一个像这样的简单矢量化卡尔曼滤波器来处理电影帧。它非常快,但目前仅限于一维输入和输出,并且它不会对任何滤波器参数进行 EM 优化。
import numpy as np
def runkalman(y, RQratio=10., meanwindow=10):
"""
A simple vectorised 1D Kalman smoother
y Input array. Smoothing is always applied over the first
dimension
RQratio An estimate of the ratio of the variance of the output
to the variance of the state. A higher RQ ratio will
result in more smoothing.
meanwindow The initial mean and variance of the output are
estimated over this number of timepoints from the start
of the array
References:
Ghahramani, Z., & Hinton, G. E. (1996). Parameter Estimation for
Linear Dynamical Systems.
http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.55.5997
Yu, B., Shenoy, K., & Sahani, M. (2004). Derivation of Extented
Kalman Filtering and Smoothing Equations.
http://www-npl.stanford.edu/~byronyu/papers/derive_eks.pdf
"""
x,vpre,vpost = forwardfilter(y, RQratio=RQratio, meanwindow=meanwindow)
x,v = backpass(x, vpre, vpost)
x[np.isnan(x)] = 0.
return x
def forwardfilter(y, RQratio=10., meanwindow=10):
"""
the Kalman forward (filter) pass:
xpost,Vpre,Vpost = forwardfilter(y)
"""
y = np.array(y,copy=False,subok=True,dtype=np.float32)
xpre = np.empty_like(y)
xpost = np.empty_like(y)
Vpre = np.empty_like(y)
Vpost = np.empty_like(y)
K = np.empty_like(y)
# initial conditions
pi0 = y[:meanwindow].mean(0)
ystd = np.std(y,0)
R = ystd * ystd
Q = R / RQratio
V0 = Q
xpre[0] = xpost[0] = pi0
Vpre[0] = Vpost[0] = V0
K[0] = 0
# loop forwards through time
for tt in xrange(1, y.shape[0]):
xpre[tt] = xpost[tt-1]
Vpre[tt] = Vpost[tt-1] + Q
K[tt] = Vpre[tt] / (Vpre[tt] + R)
xpost[tt] = xpre[tt] + K[tt] * (y[tt] - xpre[tt])
Vpost[tt] = Vpre[tt] - K[tt] * (Vpre[tt])
return xpost,Vpre,Vpost
def backpass(x, Vpre, V):
"""
the Kalman backward (smoothing) pass:
xpost,Vpost = backpass(x,Vpre,V)
"""
xpost = np.empty_like(x)
Vpost = np.empty_like(x)
J = np.empty_like(x)
xpost[-1] = x[-1]
Vpost[-1] = V[-1]
# loop backwards through time
for tt in xrange(x.shape[0]-1, 0, -1):
J[tt-1] = V[tt-1] / Vpre[tt]
xpost[tt-1] = x[tt-1] + J[tt-1] * (xpost[tt] - x[tt-1])
Vpost[tt-1] = V[tt-1] + J[tt-1] * (Vpost[tt] - Vpre[tt]) * J[tt-1]
return xpost,Vpost
如果有人知道 Python 中支持多维输入/输出和 EM 参数优化的矢量化卡尔曼平滑器实现,我很想听听!我向 PyKalman 提出了功能请求维护者,但他们表示清晰度比速度更重要。
关于python - 网格数据上的高效卡尔曼滤波器实现,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/18855387/