数量经济学的Python实现（六）|初探卡尔曼滤波器

译者：邱任翔，中国人民大学，qrx_math@ruc.edu.cn

「本部分的原文参考为：23. A First Look at the Kalman Filter^[1]」

本讲稿部分实例运行需要提前安装quantecon库

概览

本篇讲稿对于卡尔曼滤波器（Kalman Filter）进行了简单并且直观的介绍，适合以下两种人群学习：

听说过卡尔曼滤波器但是不知道原理的。
知道卡尔曼滤波器方程，但是不知道其如何得来的。
所需预备知识：熟悉矩阵运算，多元正态分布以及协方差矩阵等。

%matplotlib inline
import matplotlib.pyplot as plt
plt.rcParams['figure.figsize'] = (11, 5)  #set default figure size
from scipy import linalg
import numpy as np
import matplotlib.cm as cm
from quantecon import Kalman, LinearStateSpace
from scipy.stats import norm
from scipy.integrate import quad
from numpy.random import multivariate_normal
from scipy.linalg import eigvals

基本原理

卡尔曼滤波器在经济学中有很多应用，但现在首先假设我们是一个火箭科学家。

一枚导弹已经从Y国发射，我们的任务是追踪它。

让表示导弹当前的位置（代表经纬度坐标）。

但在当前时刻，x的准确位置我们并不知道，只能知道一些相关的信息（beliefs）。

去汇总我们所知信息的一种方法是预测

但如果总统想要知道这个导弹目前在日本海上方的概率该怎么办？
那最好用二元概率密度来总结我们的已知信息。

表示我们认为导弹在区域内的概率。密度被称为随机变量的「先验概率分布」。

为了简化我们的例子，我们假设先验概率呈高斯分布：

其中是整个分布的均值，是的协方差矩阵。在我们的模拟中，我们假设：

二元概率密度在下面显示为等高线图，红色椭圆的中心即代表。

# Set up the Gaussian prior density p
Σ = [[0.4, 0.3], [0.3, 0.45]]
Σ = np.matrix(Σ)
x_hat = np.matrix([0.2, -0.2]).T
# Define the matrices G and R from the equation y = G x + N(0, R)
G = [[1, 0], [0, 1]]
G = np.matrix(G)
R = 0.5 * Σ
# The matrices A and Q
A = [[1.2, 0], [0, -0.2]]
A = np.matrix(A)
Q = 0.3 * Σ
# The observed value of y
y = np.matrix([2.3, -1.9]).T

# Set up grid for plotting
x_grid = np.linspace(-1.5, 2.9, 100)
y_grid = np.linspace(-3.1, 1.7, 100)
X, Y = np.meshgrid(x_grid, y_grid)

def bivariate_normal(x, y, σ_x=1.0, σ_y=1.0, μ_x=0.0, μ_y=0.0, σ_xy=0.0):
    '''
    Compute and return the probability density function of bivariate normal
    distribution of normal random variables x and y

    Parameters
    ----------
    x : array_like(float)
        Random variable

    y : array_like(float)
        Random variable

    σ_x : array_like(float)
          Standard deviation of random variable x

    σ_y : array_like(float)
          Standard deviation of random variable y

    μ_x : scalar(float)
          Mean value of random variable x

    μ_y : scalar(float)
          Mean value of random variable y

    σ_xy : array_like(float)
           Covariance of random variables x and y

    '''

    x_μ = x - μ_x
    y_μ = y - μ_y

    ρ = σ_xy / (σ_x * σ_y)
    z = x_μ**2 / σ_x**2 + y_μ**2 / σ_y**2 - 2 * ρ * x_μ * y_μ / (σ_x * σ_y)
    denom = 2 * np.pi * σ_x * σ_y * np.sqrt(1 - ρ**2)
    return np.exp(-z / (2 * (1 - ρ**2))) / denom

def gen_gaussian_plot_vals(μ, C):
    'Z values for plotting the bivariate Gaussian N(μ, C)'
    m_x, m_y = float(μ[0]), float(μ[1])
    s_x, s_y = np.sqrt(C[0, 0]), np.sqrt(C[1, 1])
    s_xy = C[0, 1]
    return bivariate_normal(X, Y, s_x, s_y, m_x, m_y, s_xy)

# Plot the figure

fig, ax = plt.subplots(figsize=(10, 8))
ax.grid()

Z = gen_gaussian_plot_vals(x_hat, Σ)
ax.contourf(X, Y, Z, 6, alpha=0.6, cmap=cm.jet)
cs = ax.contour(X, Y, Z, 6, colors='black')
ax.clabel(cs, inline=1, fontsize=10)

plt.show()

滤波器的步骤

现在我们有一些好消息和一些坏消息。

好消息是导弹已经被我们的传感器定位，其报告的位置是在。

原始先验概率分布与的位置如下图所示：

fig, ax = plt.subplots(figsize=(10, 8))
ax.grid()

Z = gen_gaussian_plot_vals(x_hat, Σ)
ax.contourf(X, Y, Z, 6, alpha=0.6, cmap=cm.jet)
cs = ax.contour(X, Y, Z, 6, colors='black')
ax.clabel(cs, inline=1, fontsize=10)
ax.text(float(y[0]), float(y[1]), '$y$', fontsize=20, color='black')

plt.show()

然而坏消息是我们的传感器并不精确。

特别地，我们不应该把我们传感器的输出解释为；而应当是：

这里的和是的矩阵，并且是正定的，这两者都假定是已知的，并且噪声项与独立。

那么我们应该如何结合之前的和这个新信息明确导弹位置呢？

答案是使用贝叶斯公式，他告诉我们先验概率和后验概率之间的关系：

其中：

在求解时，我们有：

条件概率密度是
不依赖于，仅作为一个标准化的常数加入计算。

因为我们整个方程是线性的，并且涉及到的分布是高斯分布（「被称为线性高斯系统」），所以可以对总体进行线性回归来计算更新后的概率密度。

特别的，我们可以解得（具体解法参照C. M. Bishop. Pattern Recognition and Machine Learning. Springer, 2006.的93页，并且最后需要用到伍德伯里矩阵恒等式（Woodbury matrix identity））：

其中

这个新的概率密度如下图所示：

fig, ax = plt.subplots(figsize=(10, 8))
ax.grid()

Z = gen_gaussian_plot_vals(x_hat, Σ)
cs1 = ax.contour(X, Y, Z, 6, colors='black')
ax.clabel(cs1, inline=1, fontsize=10)
M = Σ * G.T * linalg.inv(G * Σ * G.T + R)
x_hat_F = x_hat + M * (y - G * x_hat)
Σ_F = Σ - M * G * Σ
new_Z = gen_gaussian_plot_vals(x_hat_F, Σ_F)
cs2 = ax.contour(X, Y, new_Z, 6, colors='black')
ax.clabel(cs2, inline=1, fontsize=10)
ax.contourf(X, Y, new_Z, 6, alpha=0.6, cmap=cm.jet)
ax.text(float(y[0]), float(y[1]), '$y$', fontsize=20, color='black')

plt.show()

上图中我们取。

预测

现在为止我们做了哪些工作呢？

我们根据先验概率密度和新的信息得到了导弹当前位置的概率。

这被称为“过滤”而不是“预测”，因为这只是过滤了噪音而并没有预测未来。

被称作过滤分布。

但是现在我们被赋予了一项新的使命——在一单位时间过去后去预测导弹的位置。

为此，我们需要建立一个状态更新的模型。

我们假设这样一个线性高斯系统：

我们的目标就是根据这个运动方程和当前的过滤分布计算出一个新的预测分布。

根据运动方程，我们所需要做的就是引入一个随机向量，并且计算出的分布情况，其中并且独立于。

由于高斯分布的线性组合仍然是高斯分布，所以也是高斯分布。

在经过一些基本计算后，我们有：

并且

矩阵通常被写作，称为「卡尔曼增益」。

下标告诉我们主要取决于，而不是或者

使用这个符号，我们将计算的结果总结如下：

更新后的概率分布是，其中：

计算过程中的参数设置如下：

fig, ax = plt.subplots(figsize=(10, 8))
ax.grid()

# Density 1
Z = gen_gaussian_plot_vals(x_hat, Σ)
cs1 = ax.contour(X, Y, Z, 6, colors='black')
ax.clabel(cs1, inline=1, fontsize=10)

# Density 2
M = Σ * G.T * linalg.inv(G * Σ * G.T + R)
x_hat_F = x_hat + M * (y - G * x_hat)
Σ_F = Σ - M * G * Σ
Z_F = gen_gaussian_plot_vals(x_hat_F, Σ_F)
cs2 = ax.contour(X, Y, Z_F, 6, colors='black')
ax.clabel(cs2, inline=1, fontsize=10)

# Density 3
new_x_hat = A * x_hat_F
new_Σ = A * Σ_F * A.T + Q
new_Z = gen_gaussian_plot_vals(new_x_hat, new_Σ)
cs3 = ax.contour(X, Y, new_Z, 6, colors='black')
ax.clabel(cs3, inline=1, fontsize=10)
ax.contourf(X, Y, new_Z, 6, alpha=0.6, cmap=cm.jet)
ax.text(float(y[0]), float(y[1]), '$y$', fontsize=20, color='black')

plt.show()

递归过程

让我们回顾一下我们的过程。

我们从当前位置的先验概率分布开始。

然后使用了当前位置的测量值去获得。

最后根据的运动方程得到了。

如果我们进行下一个周期的循环，就把作为一个新的先验概率分布。

我们用来代替表示；用来代替表示，那么整个递归的过程如下所示：

从当期的先验概率分布开始。
观察得到当期的测量值 .
应用贝叶斯公式和条件分布，根据和，计算得到过滤分布。
根据状态更新方程和过滤分布计算得到
增加一个单位，返回步骤1。

根据之前计算的结果，和的动态路径如下：

这就是标准的卡尔曼过滤器的动态方程。

收敛性

矩阵衡量了我们对的预测值的不确定性。

除了一些特殊情况，这种不确定性都会永远存在，不会被完全解决。

原因之一是我们的预测值是根据期的有效信息计算的，而不是期。

即使我们知道的准确值（事实上我们并不知道），转移方程（运动方程）意味着。

而冲击在期也并不能被观测，因此任何期对的预测都会出现一些偏误。（除非是退化的（degenerate））

然而，当时收敛到一个恒定的矩阵是完全有可能的。

我们将关于的等式扩展一下：

这是一个关于的非线性差分方程。

其不动点可以通过如下式子来求得：

上面两个式子分别被称为离散时间的「Riccati方程」和离散时间的「代数Riccati方程」。

不动点存在以及收敛的条件在D. B. O. Anderson and J. B. Moore. Optimal Filtering. Dover Publications, 2005.第四章中有详细说明。

其中一个充分（但不必要）的条件是矩阵的所有特征值满足(上书中的77页)

这个加强条件使得非条件分布当时收敛。

在这种情况下，对于任意非负对称的初值都会收敛到同样非负对称的矩阵。

实现

Quantecon库中的Kalman类可以实现卡尔曼滤波器。

需要的实例数据有：

先验概率分布的矩。
库中的另一个类LinearStateSpace实例。后者是用来表示如下的一种线性状态空间：

其中冲击和是满足IID的标准正态分布。

为了和本文中使用的符号建立练习，我们有：

Kalman类中有大量的实现函数，其中的一些我们将在系列之后的文章中学到更高级的应用时使用到。
与本文相关的函数有：

prior_to_filtered，可以根据先验分布计算出过滤分布
filtered_to_forecast，可以根据过滤分布计算出预测分布
update，将上述两个函数结合
stationary_values，计算出不动点以及相应的卡尔曼增益。

练习

我们看一下卡尔曼滤波器的简单应用：

假设

所有变量都为标量。
未知的实际位置为常数等于

因此状态更新方程中。

测量方程为,其中且满足IID。

练习的任务是使用Kalman类绘制前五个预测概率密度分布

在模拟中取。

答案

# Parameters
θ = 10  # Constant value of state x_t
A, C, G, H = 1, 0, 1, 1
ss = LinearStateSpace(A, C, G, H, mu_0=θ)

# Set prior, initialize kalman filter
x_hat_0, Σ_0 = 8, 1
kalman = Kalman(ss, x_hat_0, Σ_0)

# Draw observations of y from state space model
N = 5
x, y = ss.simulate(N)
y = y.flatten()

# Set up plot
fig, ax = plt.subplots(figsize=(10,8))
xgrid = np.linspace(θ - 5, θ + 2, 200)

for i in range(N):
    # Record the current predicted mean and variance
    m, v = [float(z) for z in (kalman.x_hat, kalman.Sigma)]
    # Plot, update filter
    ax.plot(xgrid, norm.pdf(xgrid, loc=m, scale=np.sqrt(v)), label=f'$t={i}$')
    kalman.update(y[i])

ax.set_title(f'First {N} densities when $\\theta = {θ:.1f}$')
ax.legend(loc='upper left')
plt.show()

Reference

[1]

23. A First Look at the Kalman Filter: https://python.quantecon.org/kalman.html

本站仅提供存储服务，所有内容均由用户发布，如发现有害或侵权内容，请点击举报。

目录

概览

基本原理

滤波器的步骤

预测

递归过程

收敛性

实现

练习

答案

Reference