如何使用 R 中的 DSE 包模拟卡尔曼滤波器的后验滤波估计
How to simulate the posterior filtered estimates of a Kalman Filter using the DSE package in R
如何使用 DSE 包从 R 中的卡尔曼滤波器模拟调用后验(精细)状态估计?
我在下面添加了一个示例。假设我创建了一个简单的随机游走状态 space,误差为标准正态分布。该模型是使用 SS 函数创建的,初始化状态和协方差估计为零。因此理论模型形式为:
X(t) = X(t-1) + e(t)~N(0,1) 用于状态演化
Y(t) = X(t) + w(t)~N(0,1)
我们现在按照《统计软件杂志》"Kalman Filtering in R" 文章第 6 页和第 7 页的说明在 R 中实现这一点。首先,我们使用 SS() 函数创建状态 space 模型并将其存储在名为 kalman.filter:
的变量中
kalman.filter=dse::SS(F = matrix(1,1,1),
Q = matrix(1,1,1),
H = matrix(1,1,1),
R = matrix(1,1,1),
z0 = matrix(0,1,1),
P0 = matrix(0,1,1)
)
然后我们使用 simulate() 从模型中模拟 100 个观测值并将它们放入名为 simulate 的变量中。kalman.filter:
simulate.kalman.filter=simulate(kalman.filter, start = 1, freq = 1, sampleT = 100)
然后我们 运行 使用 l() 对测量值进行卡尔曼滤波并将其存储在名为 test 的变量下:
test=l(kalman.filter, simulate.kalman.filter)
根据输出结果,哪些是我过滤后的估计值?
我找到了这个问题的答案。
首先,模型的过滤估计值没有在 l() 函数中给出。此功能仅给出领先一步的预测。我的问题的上述框架编码为:
kalman.filter=dse::SS(F = matrix(1,1,1),
Q = matrix(1,1,1),
H = matrix(1,1,1),
R = matrix(1,1,1),
z0 = matrix(0,1,1),
P0 = matrix(0,1,1)
)
simulate.kalman.filter=simulate(kalman.filter, start = 1, freq = 1, sampleT = 100)
test=l(kalman.filter, simulate.kalman.filter)
领先一步的预测是:
predictions = test$estimates$pred
以下给出了一种快速可视化的方法:
tfplot(test)
这使您可以根据实际数据快速绘制出您提前一步的预测。要获得过滤后的估计值,您需要在同一个 dse 包中使用 smoother() 函数。它输入状态模型和数据,在本例中它分别是 kalman.filter 和 simulate.kalman.filter。输出是所有时间点的平滑估计。但请注意,它是在考虑完整数据集之后执行此操作的,因此它不会在每次观察时都执行此操作。请参见下面的代码。代码的第一行为您提供平滑估计,以下几行绘制示例:
smooth = smoother(test, simulate.kalman.filter)
plot(test$estimates$pred, ylim=c(max(test$estimates$pred,smooth$filter$track,simulate.kalman.filter$outpu), min(test$estimates$pred,smooth$filter$track,simulate.kalman.filter$output)))
points(smooth$smooth$state, col = 3)
points(simulate.kalman.filter$output, col = 4)
上图绘制了所有实际数据、模型估计值和平滑估计值之间的对比。
如何使用 DSE 包从 R 中的卡尔曼滤波器模拟调用后验(精细)状态估计?
我在下面添加了一个示例。假设我创建了一个简单的随机游走状态 space,误差为标准正态分布。该模型是使用 SS 函数创建的,初始化状态和协方差估计为零。因此理论模型形式为: X(t) = X(t-1) + e(t)~N(0,1) 用于状态演化 Y(t) = X(t) + w(t)~N(0,1)
我们现在按照《统计软件杂志》"Kalman Filtering in R" 文章第 6 页和第 7 页的说明在 R 中实现这一点。首先,我们使用 SS() 函数创建状态 space 模型并将其存储在名为 kalman.filter:
的变量中kalman.filter=dse::SS(F = matrix(1,1,1),
Q = matrix(1,1,1),
H = matrix(1,1,1),
R = matrix(1,1,1),
z0 = matrix(0,1,1),
P0 = matrix(0,1,1)
)
然后我们使用 simulate() 从模型中模拟 100 个观测值并将它们放入名为 simulate 的变量中。kalman.filter:
simulate.kalman.filter=simulate(kalman.filter, start = 1, freq = 1, sampleT = 100)
然后我们 运行 使用 l() 对测量值进行卡尔曼滤波并将其存储在名为 test 的变量下:
test=l(kalman.filter, simulate.kalman.filter)
根据输出结果,哪些是我过滤后的估计值?
我找到了这个问题的答案。
首先,模型的过滤估计值没有在 l() 函数中给出。此功能仅给出领先一步的预测。我的问题的上述框架编码为:
kalman.filter=dse::SS(F = matrix(1,1,1),
Q = matrix(1,1,1),
H = matrix(1,1,1),
R = matrix(1,1,1),
z0 = matrix(0,1,1),
P0 = matrix(0,1,1)
)
simulate.kalman.filter=simulate(kalman.filter, start = 1, freq = 1, sampleT = 100)
test=l(kalman.filter, simulate.kalman.filter)
领先一步的预测是:
predictions = test$estimates$pred
以下给出了一种快速可视化的方法:
tfplot(test)
这使您可以根据实际数据快速绘制出您提前一步的预测。要获得过滤后的估计值,您需要在同一个 dse 包中使用 smoother() 函数。它输入状态模型和数据,在本例中它分别是 kalman.filter 和 simulate.kalman.filter。输出是所有时间点的平滑估计。但请注意,它是在考虑完整数据集之后执行此操作的,因此它不会在每次观察时都执行此操作。请参见下面的代码。代码的第一行为您提供平滑估计,以下几行绘制示例:
smooth = smoother(test, simulate.kalman.filter)
plot(test$estimates$pred, ylim=c(max(test$estimates$pred,smooth$filter$track,simulate.kalman.filter$outpu), min(test$estimates$pred,smooth$filter$track,simulate.kalman.filter$output)))
points(smooth$smooth$state, col = 3)
points(simulate.kalman.filter$output, col = 4)
上图绘制了所有实际数据、模型估计值和平滑估计值之间的对比。