【滤波】多元高斯

本文主要翻译自rlabbe/Kalman-and-Bayesian-Filters-in-Python的第5章节05-Multivariate-Gaussians（多元高斯）。

%matplotlib inline

#format the book
import book_format
book_format.set_style()

简介

上一篇文章中的技术非常强大，但它们只适用于一个变量或维度。它们无法表示多元数据，例如狗在野外的位置和速度。位置和速度是相互关联的，而我们永远不应该丢弃信息。在本文中，我们将学习如何从概率的角度来描述这种关系，并获得明显更好的滤波器性能。

多元正态分布

我们用高斯函数来表示一个标量随机变量，表示为$N(\mu ,{\sigma }^2)$。一个更正式的术语是：一元正态分布。

多元正态分布会是什么样子呢？多元是指多个变量。我们的目标是能够表示多元正态分布。我指的不仅仅只是空间维度——比如在$(x,y,z)$坐标系下中跟踪飞机的位置、速度和加速度，这就给出了一个九维问题。考虑一个二维的例子：它可能是机器人的$x$和$y$坐标，可能是狗在$x$轴上的位置和速度，也可能是奶牛场的产奶量和喂奶量。这些都不重要，我们可以看到，对于$N$维问题，我们需要$N$个值来描述，我们可以将其排列在列矩阵（向量）中，如下所示：

$$\mu =\left[\begin{array}{c}{\mu }_1\\ {\mu }_2\\ ...\\ {\mu }_N\end{array}\right]$$

假设一个二维位置为$x=2$，$y=17$。我们会有：

$$\mu =\left[\begin{array}{c}2\\ 17\end{array}\right]$$

下一步是表示方差。乍一看，我们可能认为我们也需要为$N$个维度设置$N$个方差。我们可以说$x$的方差是10，$y$的方差是4，就像这样：

$${\sigma }^2=\left[\begin{array}{c}10\\ 4\end{array}\right]$$

但这是不完整的，因为它没有考虑更一般的情况。在高斯一章中，我们计算了学生身高的方差。这是一个衡量高度如何相对变化的指标。如果所有学生的身高都一样，那么方差为0，如果他们的身高相差很大，那么方差会很大。

我们知道，身高和体重是存在关系的。一般来说，高个子比矮个子重，身高和体重是相关的。我们需要一种方法来表达我们所认为的身高和体重的差异，以及它们之间的关联程度。换言之，我们想知道体重与身高相比是如何变化的。我们称之为协方差。

在我们理解多元正态分布之前，我们需要理解相关性和协方差背后的数学。

相关性与相关变异

协方差描述了两个变量一起变化的程度，是相关方差的缩写。换言之，方差是衡量一个变量本身是如何变化的指标，协方差是衡量两个变量之间相互变化关系的指标。例如，随着身高的增加，体重也普遍增加。这些变量是相关的，并且是正相关的，因为一个变量越大，另一个变量也越大。随着室外温度的降低，家庭取暖费用也随之增加。这些变量呈负相关的，因为当一个变量变大时，另一个变量变小。茶叶的价格和我的狗摇尾巴的次数彼此没有关系，我们说它们是不相关或独立的——每一个都可以独立地改变。

相关性会对预测产生影响：如果你明显比我高，我可以预测你也比我重；冬天来了，我预计我会花更多的钱来取暖；如果我的狗再摇尾巴，我不认为茶叶价格会改变。

例如，这里有一个关于学校田径队学生身高和体重的曲线图。如果一个学生是68英寸高，我可以预测他们大约160磅重。因为相关性并不完美，我的预测也不完美。

from kf_book.gaussian_internal import plot_correlated_dataheight = [60, 62, 63, 65, 65.1, 68, 69, 70, 72, 74]
weight = [95, 120, 127, 119, 151, 143, 173, 171, 180, 210]
plot_correlated_data(height, weight, 'Height (in)', 'Weight (lbs)', False)

在本文中，我们只考虑线性相关。我们假设变量之间的关系是线性的。我已把上表中的数据拟合成一条直线。非线性相关的概念是存在的，但我们不会使用它。

$X$和$Y$之间的协方差方程为：

$$COV(X,Y)={\sigma }_{xy}=\mathbb{E}[(X-{\mu }_x)(Y-{\mu }_y)]$$

其中，$\mathbb{E}[X]$是$X$的期望值，定义为：

$$\mathbb{E}[X]=\{\begin{array}{cc}{\sum }_{i=1}^n{p}_i{x}_i& 离散\\ {\int }_{-\infty }^{+\infty }f(x)xdx& 连续\end{array}$$

我们假设每个数据点的概率相等，因此每个数据点的概率为$\frac{1}{N}$。对于离散情况，我们有

$$\mathbb{E}[X]=\frac{1}{N}\sum _{i=1}^n{x}_i$$

比较协方差方程和方差方程，如你所见，它们非常相似：

$$VAR(X)={\sigma }_x^2=\mathbb{E}[(X-\mu {)}^2]$$

$$COV(X,Y)={\sigma }_{xy}=\mathbb{E}[(X-{\mu }_x)(Y-{\mu }_y)]$$

特别是，如果你计算$COV(X,X)$，你会得到$VAR(X)$的方程。这支持我的说法，即方差计算随机变量如何在自身之间变化。

我们使用协方差矩阵来表示多元正态分布的协方差，它如下所示：

$$\sum =\left[\begin{array}{cccc}{\sigma }_1^2& {\sigma }_{12}& ...& {\sigma }_{1n}\\ {\sigma }_{21}& {\sigma }_2^2& ...& {\sigma }_{2n}\\ ⋮& ⋮& \ddots & ⋮\\ {\sigma }_{n1}& {\sigma }_{n2}& ...& {\sigma }_n^2\end{array}\right]$$

对角线包含每个变量的方差，非对角线元素包含第$i$个和第$j$个变量之间的协方差。所以${\sigma }_3^2$是第三个变量的方差，${\sigma }_{13}$是第一个和第三个变量之间的协方差。

协方差为0表示没有相关性。如果$x$的方差是10，$y$的方差是4，并且$x$和$y$之间没有线性相关，那么我们将写

$$\sum =\left[\begin{array}{cc}10& 0\\ 0& 4\end{array}\right]$$

如果$x$和$y$之间有少量的正相关，我们可能会

$$\sum =\left[\begin{array}{cc}10& 1.2\\ 1.2& 4\end{array}\right]$$

其中，1.2是$x$和$y$之间的协方差。我说相关性很小，因为1.2的协方差相对于10的方差很小。如果$x$和$y$之间有大量的负相关性，我们可能有

$$\sum =\left[\begin{array}{cc}10& -9.7\\ -9.7& 4\end{array}\right]$$

协方差矩阵是对称的。毕竟，$x$和$y$之间的协方差总是等于$y$和$x$之间的协方差。也就是说，对于任何$x$和$y$，${\sigma }_{xy}={\sigma }_{yx}$。

我担心你可能会困扰，所以让我们举个例子。在高斯一章中，我们有一个学生的身高$H=[1.8,2.0,1.7,1.9,1.6]$米。我们计算：

$$VAR(H)=E[(H-{\mu }_H{)}^2]=\frac{1}{N}\sum _{i=1}^n(H_i-{\mu }_H{)}^2$$

$$=\frac{1}{5}[(1.8-1.8{)}^2+(2-1.8{)}^2+(1.7-1.8{)}^2+(1.9-1.8{)}^2+(1.6-1.8{)}^2]=0.02$$

很简单，对吧？如果我们给学生称重，我们可能会发现他们的体重是$W=[70.1,91.2,59.5,93.2,53.5]$。我们能用协方差方程来建立协方差矩阵吗？当然。它看起来像：

$$\sum =\left[\begin{array}{cc}{\sigma }_H^2& {\sigma }_{HW}\\ {\sigma }_{WH}& {\sigma }_W^2\end{array}\right]$$

我们刚刚计算了高度的方差，它会出现在矩阵的左上角，右下角应该是体重的方差。用同样的方程式我们得到：

$${\mu }_W=\frac{1}{5}(70.1+91.2+59.5+93.2+53.5)=73.5$$

$${\sigma }_W^2=\frac{1}{5}[(70.1-73.5{)}^2+(91.2-73.5{)}^2+(59.5-73.5{)}^2+(93.2-73.5{)}^2+(53.5-73.5{)}^2]=261.8$$

现在是协方差。使用上述公式，我们计算：

$${\sigma }_{HW}=\mathbb{E}[(H-{\mu }_H)(W-{\mu }_W)]=\frac{1}{N}\sum _{i=1}^n(H_i-{\mu }_H)(W_i-{\mu }_W)$$

$$=\frac{1}{5}[(1.8-1.8)(70.1-73.5)+(2-1.8)(91.2-73.5)+(1.7-1.8)(59.5-73.5)+...+(1.9-1.8)(93.2-73.5)+(1.6-1.8)(53.5-73.5)]=2.18$$

这尽管很容易，但更乏味。我们再也不会这样做了，因为NumPy会帮你计算。

import numpy as npW = [70.1, 91.2, 59.5, 93.2, 53.5]
H = [1.8, 2.0, 1.7, 1.9, 1.6]
np.cov(H, W)

array([[  0.025,   2.727],[  2.727, 327.235]])

这与我们的计算不符！是出了什么问题么？其实NumPy会对小样本数据进行校正：它使用$\frac{1}{N-1}$作为标准化项，而不是$\frac{1}{N}$。

这超出了这本书的范围。简单地说，假设实际的班级规模是200名学生，我们选取了5名学生作为样本来进行计算，因为我们无法对所有200名学生进行测高和称重。但几乎可以肯定的是，我们的估计量会有一些误差，因为样本不可能完全代表总体。当我们的样本量接近200时，误差将接近0。我们说后者没有偏差，我们有一个无偏估计。相比之下，当我们采取小样本有偏差（误差是非零），我们有一个有偏差的估计。

如果误差为零，则除以$N$是有意义的。但对于有偏估计量，我们使用$\frac{1}{N-1}$对小样本进行校正。NumPy默认是这样做的，因为在实践中，我们几乎总是从更大的集合中提取数据样本。如果你需要使用Numpy计算无偏估计量，请在调用中使用bias=1。

np.cov(H, W, bias=1)

array([[  0.02 ,   2.182],[  2.182, 261.788]])

这与我们的计算一致。在本文中，我们将不再使用bias=1。因为我们使用的是随机变量，而这些随机变量是从我们跟踪的对象的无限多个位置集合中采样的。

这个矩阵告诉我们什么？它告诉我们身高的方差是$0.02m^2$，体重的方差是$261.788k{g}^2$。此外，它告诉我们体重和身高是正相关的——随着身高的增加，体重也会增加。

让我们创建完全相关的数据。我的意思是，数据完全符合一条直线——这条直线没有变化。

X = np.linspace(1, 10, 100)
Y = np.linspace(1, 10, 100)
np.cov(X, Y)

array([[6.956, 6.956],[6.956, 6.956]])

从协方差矩阵可以看出，协方差等于$x$和$y$的方差。

现在首先使$Y$和$X$呈现负相关性，再给它添加一些噪声，使它们不再完全相关。

X = np.linspace(1, 10, 100)
Y = -(np.linspace(1, 5, 100) + np.sin(X)*.2)
plot_correlated_data(X, Y)
print(np.cov(X, Y))

[[ 6.956 -3.084][-3.084  1.387]]

此时，数据不再形成一条直线。协方差为${\sigma }_{xy}=-3.08$，与${\sigma }_x^2$和${\sigma }_y^2$的量级相比，它并不接近于零，因此我们仍然觉得存在一定的相关性。我们可以通过查看图表来验证这一点，这些数据几乎形成一条直线。

现在我将随机噪声添加到一条直线上：

from numpy.random import randnX = np.linspace(1, 10, 1000) + randn(1000)*2
Y = np.linspace(1, 5, 1000) + randn(1000)
plot_correlated_data(X, Y)
print(np.cov(X, Y))

[[10.433  2.934][ 2.934  2.317]]

我们看到协方差相对于方差较小，反映了$X$和$Y$之间的相关性较低。尽管我们仍然可以通过这些数据拟合出一条直线，但数据中的变化还是较大的。

最后，如果是完全随机数据之间的协方差，又会变成什么样呢？

X = randn(100000)
Y = randn(100000)
plot_correlated_data(X, Y)
print(np.cov(X, Y))

[[1.005 0.005][0.005 1.005]]

此时的协方差，非常接近于零。从图中可以看出，没有明确的方法来绘制适合数据的直线。一条垂直线和我展示的水平线一样，都没有说服力。

多元正态分布方程

回想一下一元正态分布方程：

$$f(x,\mu ,{\sigma }^2)=\frac{1}{\sqrt[]{2\pi {\sigma }^2}}exp[-\frac{1}{2}(x-\mu {)}^2/{\sigma }^2]$$

这是$n$维的多元正态分布：

$$f(\mathbf{x},\mu ,\sum )=\frac{1}{\sqrt[]{(2\pi {)}^n|\sum |}}exp[-\frac{1}{2}(\mathbf{x}-\mu {)}^T{\sum }_{}^{-1}(\mathbf{x}-\mu )]$$

多元版本只是用矩阵代替一元版本的标量。如果你相当精通线性代数，这个方程看起来应该很容易处理。如果你并不怎么了解线性代数，也不用担心，FilterPy和SciPy都提供函数来计算它。让我们暂时忽略如何计算，先绘制它，看看它是什么样子。

import kf_book.mkf_internal as mkf_internalmean = [2., 17.]
cov = [[10., 0.], [0., 4.]]mkf_internal.plot_3d_covariance(mean, cov)

这是一个多元高斯图，平均值$\mu =\left[\begin{array}{c}2\\ 17\end{array}\right]$，协方差$\sum =\left[\begin{array}{cc}10& 0\\ 0& 4\end{array}\right]$。三维形状的z轴坐标表示$(X,Y)$值对应的概率密度，我可以把$x$和$y$的值投影到图表的墙上，你可以看到它们呈现高斯钟形曲线的形状。$X$的曲线比$Y$的曲线宽，可以用${\sigma }_x^2=10$和${\sigma }_y^2=4$来解释。三维曲面的最高点位于$X$和$Y$的均值处。

所有的多元高斯函数都是这个形状。如果我们认为这是狗的二维位置的高斯分布，每个$(X,Y)$的点对应的$z$值就是狗在那个位置的概率密度。严格地说，这是联合概率密度函数，我将很快定义它。因此，狗接近$(2,17)$的概率最高，接近$(5,14)$的概率适中，接近$(10,10)$的概率很低。和一元的情况一样，这是一个概率密度，而不是一个概率。连续分布的范围是无限的，因此精确到$(2,17)$，或任何其他点，它们的概率都是$0\%$。我们可以通过用积分计算表面下的体积来计算在给定范围内的概率。

FilterPy的filterpy.stats模块使用multivariate_gaussian()实现多元高斯分布。SciPy的stats模块用multivariate_normal()来实现。在使用它时，你只需要设置一次均值和协方差，然后就可以计算任意数量的$x$值的概率密度，并可以任意调用次数。我将函数命名为multivariate_gaussian()，以确保它不会与SciPy的版本混淆。

from filterpy.stats import gaussian, multivariate_gaussian

我将演示如何使用它，然后继续讨论更有趣的事情。

首先，让我们计算一下：如果我们认为狗在$(2,7)$，而事实上它在$(2.5,7.3)$的概率密度。其中，$x$的方差为8，$y$的方差为3。

首先，将$x$设置为$(2.5,7.3)$，可以使用元组、列表或NumPy数组。

x = [2.5, 7.3]

接下来，我们设定均值：

mu = [2.0, 7.0]

最后，我们必须定义协方差矩阵。在问题陈述中，我们没有提到$x$和$y$之间的任何相关性，我们将假设没有相关性。这是有道理的：狗可以选择在$x$方向或$y$方向上任意移动而不影响另一个方向。我将使用变量名$\mathbf{P}$。卡尔曼滤波器使用$\mathbf{P}$表示协方差矩阵，我们需要熟悉这些约定。

P = [[8., 0.], [0., 3.]]

现在调用函数：

%precision 4
multivariate_gaussian(x, mu, P)

0.0315

我们可以从scipy.stats模块得到同样的结果：

import scipy
from scipy.stats import multivariate_normalprint(f'{multivariate_normal(mu, P).pdf(x):.4f}')

0.0315

是时候定义一些术语了：

联合概率$P(x,y)$，指$x$和$y$都发生的概率。例如，如果你滚动两个骰子，$P(2,5)$表示第一个骰子的结果是2、第二个骰子的结果是5的概率。假设骰子是六面且公平的，概率$P(2,5)=\frac{1}{6}\times \frac{1}{6}=\frac{1}{36}$。上面的3D图表，显示的就是联合概率的密度函数。

边际概率指一个事件发生的概率，不考虑任何其他事件。在上面的图表中，投影画在左边的高斯曲线是$Y$的边缘，这是狗忽略$X$的值后，在$Y$的任何位置的概率。早些时候我写过，把$x$和$y$的值投影到图表的墙上，这些是$x$和$y$的边际概率。高斯的另一个计算优势是，多元高斯的边缘概率也是高斯！

让我们用一种稍微不同的方式来看待这个问题。我将生成均值为$\left[\begin{array}{c}2\\ 7\end{array}\right]$，协方差为$\left[\begin{array}{cc}8& 0\\ 0& 3\end{array}\right]$的二元高斯分布的1000个点，而不是绘制显示概率分布的曲面。

mkf_internal.plot_3d_sampled_covariance(mu, P)

我们可以认为：这些由给定均值和协方差而生成的采样点，都是我们的狗可能的位置。侧面的等高线显示了$X$和$Y$的边际概率。我们可以看到，它更可能在$(2,7)$，因为有很多点，密度很大；而不是在$(-5,5)$，因为仅有几个点，很稀疏。

尽管这些图很美，但很难从中得到有用的信息。例如，很难判断$X$和$Y$的方差是否相同，以及它们之间的相关性有多大。因此，在这本书的大部分内容中，我都会将高斯函数显示为等高线图。

等高线图显示了多元高斯函数在特定标准差下的取值范围，看起来就像从3D绘图中取出的水平切片。

下面的图显示了3个不同协方差矩阵的标准差的切片形状。

mkf_internal.plot_3_covariances()

对于在线或在电脑上的Juptyer笔记本中查看此内容的人，这里有一个在保持方差不变的同时，改变协方差的动画：

这些图看起来像圆和椭圆。事实上，通过多元高斯函数的任何切片都是椭圆。因此，在统计学中，我们不称之为等高线图，而是称之为误差椭圆或置信椭圆。

此代码使用filterpy.stats模块的plot_covariance_ellipse()函数。默认情况下，函数显示一个标准差，但你可以使用variance或std参数来控制显示的内容。例如，variance=3**2或std=3将显示第3个标准偏差，variance=[1, 4, 9]或std=[1, 2, 3]将显示第1、第2和第3个标准偏差。

from filterpy.stats import plot_covariance_ellipse
import matplotlib.pyplot as pltP = [[2, 0], [0, 6]]
plot_covariance_ellipse((2, 7), P, fc='g', alpha=0.2, std=[1, 2, 3],title='|2 0|\n|0 6|')
plt.gca().grid(b=False)

纯色可能会误导你：同一标准差之间的概率分布是均匀的，其实这是不正确的。正如你可以从高斯曲线的3D图中看出的。下图是协方差为$\left[\begin{array}{cc}2& 1.2\\ 1.2& 1.3\end{array}\right]$的概率分布的2D着色表示，较深的灰色对应较高的概率密度。

from kf_book.nonlinear_plots import plot_cov_ellipse_colormapplot_cov_ellipse_colormap(cov=[[2, 1.2], [1.2, 1.3]])

思考这些图的物理解释，就可以阐明它们的意义。第一个图的均值和协方差为：

$$\mu =\left[\begin{array}{c}2\\ 7\end{array}\right]$$

$$\sum =\left[\begin{array}{cc}2& 0\\ 0& 2\end{array}\right]$$

x = [2, 7]
P = [[2, 0], [0, 2]]
plot_covariance_ellipse(x, P, fc='g', alpha=0.2, title='|2 0|\n|0 2|')
plt.gca().grid(b=False)

一种贝叶斯的思考方式是，椭圆向我们展示了我们概率中的误差程度。一个小圆圈表示我们有一个非常小的误差，一个非常大的圆圈表示我们的概率有很多误差。椭圆的形状向我们展示了$X$和$Y$中误差的几何关系。这里我们是一个圆，所以$X$和$Y$中的误差是同样可能的。

第二个图的均值和协方差是：

$$\mu =\left[\begin{array}{c}2\\ 7\end{array}\right]$$

$$\sum =\left[\begin{array}{cc}2& 0\\ 0& 6\end{array}\right]$$

x = [2, 7]
P = [[2, 0], [0, 6]]
plot_covariance_ellipse(x, P, fc='g', alpha=0.2, title='|2 0|\n|0 6|')
plt.gca().grid(b=False)

这次我们对$X$（${\sigma }_x^2=2$）和$Y$（${\sigma }_y^2=6$），使用不同的方差，结果是一个又高又窄的椭圆。我们可以看到$Y$和$X$之间有更多的不确定性。在这两种情况下，我们都相信狗是在$(2,7)$，但误差程度（不确定性）是不同的。

第三个图的均值和协方差是：

$$\mu =\left[\begin{array}{c}2\\ 7\end{array}\right]$$

$$\sum =\left[\begin{array}{cc}2& 1.2\\ 1.2& 6\end{array}\right]$$

x = [2, 7]
P = [[2, 1.2], [1.2, 2]]
plot_covariance_ellipse(x, P, fc='g', alpha=0.2, title='|2 1.2|\n|1.2 2|')

这是协方差矩阵的非对角元素中有值的等高线，此时等高线图的形状是倾斜的椭圆。这不是巧合。一个倾斜的椭圆告诉我们，$x$和$y$的值在某种程度上是相关的。协方差矩阵中的非对角元素不为零，表示存在相关性。

回想一下，上文中的身高和体重的图，它形成了一组倾斜的点。我们可以计算两个或更多变量的协方差，通过把它们放入一个二维数组，然后调用NumPy的cov()函数。让我们这样做，然后在数据顶部绘制$2\sigma$协方差椭圆。我们需要使用bias=1，因为数据代表了整个总体，它不是一个样本。

cov_hw = np.cov(np.vstack((height, weight)), bias=1)
cov_hw

array([[  18.5249,  135.701 ],[ 135.701 , 1092.29  ]])

plt.scatter(height, weight, s=120, marker='s')
plt.title('Track Team Height vs. Weight')
plt.xlabel('Height (in)'); plt.ylabel('Weight (lbs)')
plot_covariance_ellipse((np.mean(height), np.mean(weight)), cov_hw, fc='g', alpha=0.2, axis_equal=False, std=2)

这将有助于你对协方差矩阵的含义和用法形成强烈的直觉。协方差椭圆显示了数据是如何相互分散的，像这样的窄的斜椭圆告诉你：数据是非常相关的。对于任何给定的高度，重量的范围都很窄。椭圆向右倾斜，告诉我们正相关——随着$x$的增加，$y$也随之增加；如果椭圆向左倾斜，那么相关性将为负——随着$x$的增加，$y$也随之减少。我们可以在下面的图中看到这一点：

max_temp = [200, 250, 300, 400, 450, 500]
lifespan = [10, 9.7, 5, 5.4, 4.3, 0.3]plt.scatter(max_temp, lifespan, s=80)
cov = np.cov(np.vstack((max_temp, lifespan)))
plot_covariance_ellipse((np.mean(max_temp), np.mean(lifespan)), cov, fc='g', alpha=0.2, axis_equal=False, std=2)
plt.title('Engine Temperature vs Lifespan')
plt.xlabel('Temperature (C)'); plt.ylabel('Years')

方差和协方差之间的关系很难通过死记硬背来弄清楚，所以这里有一个交互图：

from ipywidgets import interact
from kf_book.book_plots import figsize, FloatSliderfig = None
def plot_covariance(var_x, var_y, cov_xy):global figif fig: plt.close(fig)fig = plt.figure(figsize=(4,4))P1 = [[var_x, cov_xy], [cov_xy, var_y]]plot_covariance_ellipse((10, 10), P1, axis_equal=False,show_semiaxis=True)plt.xlim(4, 16)plt.gca().set_aspect('equal')plt.ylim(4, 16)with figsize(y=6):interact (plot_covariance,           var_x=FloatSlider(5, min=0, max=20), var_y=FloatSlider(5, min=0, max=20), cov_xy=FloatSlider(1.5, min=0, max=50, step=.2))

皮尔逊相关系数

我们在本文中不使用这个系数，但你可能在别处看到它。如果你不感兴趣，可以跳过这一部分。

两个变量之间的相关性可以用皮尔逊相关系数给出一个数值，它定义为：

$${\rho }_{xy}=\frac{COV(X,Y)}{{\sigma }_x{\sigma }_y}$$

此值的范围为-1到1。如果协方差为0，则$\rho =0$。如果$\rho$大于0的值表示该关系为正相关，负值表示存在负相关；接近-1或1的值表示非常强的相关性，接近0的值表示非常弱的相关性。

相关性和协方差是密切相关的。协方差有与其相关的单位，而相关性是一个无单位的比率。例如，我们的狗${\sigma }_{xy}$的单位是$m^2$。

我们可以用scipy.stats.pearsonr函数来计算皮尔逊相关系数。它返回一个元组，元组的第一个值就是皮尔逊系数。这里我们计算学生运动员身高与体重的$\rho$：

from scipy.stats import pearsonrpearsonr(height, weight)[0]

0.9539731096080194

同样，我们也计算一下发动机温度和寿命之间的关系。

pearsonr(max_temp, lifespan)[0]

-0.9178223453527254

利用相关性改进估计

假设我们相信狗所在的位置是$(5,10)$，并且有一定的协方差。如果$x$和$y$的标准差各为2，且它们是强相关的，协方差等值线可能如下所示：

P = [[4, 3.9], [3.9, 4]]plot_covariance_ellipse((5, 10), P, ec='k', std=[1, 2, 3])
plt.xlabel('X')
plt.ylabel('Y')

现在我告诉你$x=7.5$，那关于$y$的值，我们能推断出什么？该位置极有可能位于$3\sigma$协方差椭圆内。我们可以根据协方差矩阵推断$y$中的位置，因为$x$和$y$之间存在相关性。我将$y$的可能值范围用蓝色填充的圆表示。

mkf_internal.plot_correlation_covariance()

这个圆在数学上并不正确，但它能让人理解这个意思。目前，我们可以预测$y$可能接近12，$y=-10$的值极不可能。

一些关于相关性和独立性的辨析。如果变量是独立的，意味着它们可以单独变化。如果你走在开阔的田野上，你可以沿着$x$方向（东西方向）、$y$方向（南北方向）或它们的任意组合移动。也就是说，如果变量是独立的，那么肯定是不相关的。如果变量是不相关的，但是却有可能是独立，有可能是不独立的。例如，考虑$y=x^2$。相关性是一个线性度量，因此$x$和$y$是不相关的。但是，$y$依赖于$x$，并不独立。

多元高斯乘

在上一章中，我们将不确定的观测值与不确定的先验值相结合，将它们的高斯相乘，结果是另一个方差较小的高斯分布。即：如果两个都不太确定的信息相互印证，我们就可以得到更确定的信息。图表如下所示：

mkf_internal.plot_gaussian_multiply()

观测值1和观测值2的组合产生了更大的确定性，因此新的高斯分布更高、更窄——方差变得更小。同样的情况也发生在多元高斯分布中。

下面是多元高斯数相乘的方程。$\sum$表示这些是矩阵，而不是标量。具体来说，它们是协方差矩阵：

$$\mu ={\sum }_2^{}({\sum }_1^{}+{\sum }_2^{}{)}^{-1}{\mu }_1+{\sum }_1^{}({\sum }_1^{}+{\sum }_2^{}{)}^{-1}{\mu }_2$$

$${\sum }_{}^{}={\sum }_1^{}({\sum }_1^{}+{\sum }_2^{}{)}^{-1}{\sum }_2^{}$$

它们是通过将先验值和观测值插入贝叶斯定理而产生的多元高斯函数。

当然，你不需要记住这些方程。因为它们将在卡尔曼滤波方程中计算，稍后将介绍这些方程。FilterPy中还提供了使用multivariate_multiply()来进行此计算，你也可以从filterpy.stats中导入。

为了给你一些关于这方面的直觉，回想一下一元高斯乘的方程：

$$\mu =\frac{{\sigma }_2^2{\mu }_1+{\sigma }_1^2{\mu }_2}{{\sigma }_1^2+{\sigma }_2^2}$$

$${\sigma }^2=\frac{{\sigma }_1^2{\sigma }_2^2}{{\sigma }_1^2+{\sigma }_2^2}$$

多元高斯积的结果和这个，从结构上看很类似。如果你认识到用-1次方来表示矩阵的逆，就像是一个倒数，因为$A{A}^{-1}=I$。我将把矩阵的逆重写为除法——这在数学上是不正确的，因为矩阵的除法没有定义，但它确实帮助我们比较方程。

$$\mu \approx \frac{{\sum }_2^{}{\mu }_1+{\sum }_1^{}{\mu }_2}{{\sum }_1^{}+{\sum }_2^{}}$$

$${\sum }_{}^{}\approx \frac{{\sum }_1^{}{\sum }_2^{}}{{\sum }_1^{}+{\sum }_2^{}}$$

在这种形式下，一元和多元方程之间的关系是很清晰的。

现在让我们用一个具体的例子来探讨多元高斯函数。假设我们用两个雷达系统跟踪一架飞机，忽略高度，这样我就可以使用二维坐标来表示该飞机的位置。雷达能提供的是：目标的距离和方向。由于我们刚开始并不确定飞机的位置，所以协方差（对位置的不确定性）可能如下图所示。在贝叶斯统计语言中，即先验值。

P0 = [[6, 0], [0, 6]]
plot_covariance_ellipse((10, 10), P0, fc='y', alpha=0.6)

现在假设飞机左下角有一个雷达，进一步假设雷达的方向观测是准确的，但是距离观测却不那么准确。观测误差的协方差可能如下所示（以绿色绘制在黄色之前的顶部）：

P1 = [[2, 1.9], [1.9, 2]]
plot_covariance_ellipse((10, 10), P0, fc='y', alpha=0.6)
plot_covariance_ellipse((10, 10), P1, fc='g', alpha=0.9)

仔细分析一下图像的含义：椭圆指向雷达。它是非常长的，因为距离观测是不准确的，飞机可能有一个相当大的距离观测范围；它是非常窄的，因为方向观测非常准确，因此飞机必须非常接近方位估计。

我们想找到后验值，也就是把观测值和先验值结合起来，得到新的均值和协方差。就像其他章节一样，我们通过将两者相乘来组合。

from filterpy.stats import multivariate_multiplyP2 = multivariate_multiply((10, 10), P0, (10, 10), P1)[1]
plot_covariance_ellipse((10, 10), P0, ec='k', fc='y', alpha=0.2)
plot_covariance_ellipse((10, 10), P1, ec='k', fc='g', alpha=0.9)
plot_covariance_ellipse((10, 10), P2, ec='k', fc='b')

我用非常透明的黄色绘制了原始值（先验值），用绿色绘制了雷达读数（观测值），用蓝色绘制了最终值（后验值）。

后验保留了与雷达观测相同的形状和位置，但大小更小。我们已经在一维高斯中看到了这一点：两个高斯相乘使得方差变小，因为我们包含了更多的信息，因此我们的不确定性降低。另一点是，协方差的形状反映了飞机和雷达的物理布局。这一点的重要性将在下一步变得更加明显。

现在让我们从第二个雷达得到一个观测值，这个雷达位于飞机的右下角。上一步的后验值变成了新的先验值，我用黄色绘制，而新的观测值用绿色表示。

P3 = [[2, -1.9], [-1.9, 2.2]]
plot_covariance_ellipse((10, 10), P2, ec='k', fc='y', alpha=0.6)
plot_covariance_ellipse((10, 10), P3, ec='k', fc='g', alpha=0.6)

我们通过高斯相乘来合并这些信息：

P4 = multivariate_multiply((10, 10), P2, (10, 10), P3)[1]
plot_covariance_ellipse((10, 10), P2, ec='k', fc='y', alpha=0.6)
plot_covariance_ellipse((10, 10), P3, ec='k', fc='g', alpha=0.6)
plot_covariance_ellipse((10, 10), P4, ec='k', fc='b')

飞机唯一可能的位置就是，两个椭圆相交的地方。由先验值和观测值相乘形成的交集是一种新的高斯分布。这使我们能够对飞机进行三角测量，从而得出非常准确的估计。我们没有显式地编写任何代码来执行三角测量，这是将每个观测值的高斯数直接相乘的自然结果。

让我们考虑一个不同的雷达布局。假设第一个雷达正对着飞机的左边。我可以用：

$${\sum }_{}^{}=\left[\begin{array}{cc}2& 0\\ 0& 0.2\end{array}\right]$$

这里我们看到的是先验值与观测值相乘的结果：

P1 = [[2, 0], [0, .2]]
P2 = multivariate_multiply((10, 10), P0, (10, 10), P1)[1]
plot_covariance_ellipse((10, 10), P0, ec='k', fc='y', alpha=0.2)
plot_covariance_ellipse((10, 10), P1, ec='k', fc='g', alpha=0.6)
plot_covariance_ellipse((10, 10), P2, ec='k', fc='b')

现在我们可以合并第二个雷达系统的观测结果，我们将把它放在与以前相同的位置：

P3 = [[2, -1.9], [-1.9, 2.2]]
P4 = multivariate_multiply((10, 10), P2, (10, 10), P3)[1]
plot_covariance_ellipse((10, 10), P2, ec='k', fc='y', alpha=0.2)
plot_covariance_ellipse((10, 10), P3, ec='k', fc='g', alpha=0.6)
plot_covariance_ellipse((10, 10), P4, ec='k', fc='b')

我们的估计不如前一个例子准确：两个雷达系统相对于飞机的位置不再相互正交，因此这种三角测量方式不是最优的。

最后一个例子是，想象一下最糟糕的情况，就是两个雷达系统重合。这和在短时间内从同一个雷达上进行两次观测，效果应该是差不多的。

P5 = multivariate_multiply((10,10), P2, (10.1, 9.97), P2)
plot_covariance_ellipse((10, 10), P2, ec='k', fc='y', alpha=0.2)
plot_covariance_ellipse((10.1, 9.97), P2, ec='k', fc='g', alpha=0.6)
plot_covariance_ellipse(P5[0], P5[1], ec='k', fc='b')
plt.xlim(6, 14)

隐藏量

你可能已经能够明白了为什么多元卡尔曼滤波器比一元卡尔曼滤波器性能更好：变量之间的相关性可以显著改善我们的估计。我们还可以更进一步。

假设我们正在跟踪一架飞机，我们得到了时间$t=$ 1、2和3秒时，$x$和$y$坐标。坐标如图，你的直觉告诉你在$t=4$秒时，$x$的值是多少？

mkf_internal.show_position_chart()

看起来飞机是直线飞行的，我们知道飞机一秒钟之内转不成弯。因此，最合理的猜测是，在$t=4$时，飞机处于$(4,4)$。我将用绿色箭头来描绘它。

mkf_internal.show_position_prediction_chart()

这样推断是因为，假定飞机的速度是恒定的（速度的大小和方向）。因此，飞机在每个时间步的$x$和$y$方向各移动相同的距离。

如果我们考虑速度，在这种情况下，我们跟踪一架飞机的位置和速度。当然，我不可能绘制4D图，包括$x$、$y$，以及它们各自方向上的速度，所以让我们只在$x$维度上进行绘制。

在时间0，我们可能会确定位置（$x=0$），却不知道速度。我们可以用这样的协方差矩阵来描绘：较窄的宽度表示我们对位置的相对确定，而较高的高度表示我们对速度缺乏了解。

mkf_internal.show_x_error_chart(1)

然而，位置和速度是相关的。如果速度是$5m/s$，那么在1秒内位置将是$5m$。如果速度是$-10m/s$，那么在1秒内位置将是$-10m$。让我们用绘制在对角线上的协方差来可视化1秒后的预测：

mkf_internal.show_x_error_chart(2)

因为不知道速度是多少，我们并不能根据该协方差来预测一个新的位置。但是一秒钟后我们得到了$x=5$的位置观测的更新。

mkf_internal.show_x_error_chart(3)

对于新的观测，由于我们对距离还是相对比较确定的，因此$x$方向比较窄。当然，我们对速度的不确定性意味着在$y$轴上非常分散。但正如我之前所说，位置与速度有关。

因此，这两个协方差的叠加的地方，就是需要关注的地方。时间$t=1$（$x=5$）的唯一合理估计大致是，预测协方差和观测协方差之间的交集！更确切地说，我们可以使用上文的数学原理，将两个协方差相乘。从贝叶斯的观点来看，我们将先验值与观测值的概率（可能性）相乘，得到后验概率。如果我们使用贝叶斯方程将两个协方差相乘，我们得到以下结果：

mkf_internal.show_x_error_chart(4)

新的协方差（后验）位于原本两个协方差的交点处。它略微倾斜，表明位置和速度之间有一定的相关性。更重要的是，它比原本的两个协方差都小得多。

在上一章中，每次执行update()时，我们的方差都会变小，同样的情况也发生在这里。然而，这里的改善明显更好。这是因为我们使用了两种不同的信息，但它们是相互关联的。由于速度和位置是相互关联的，那么大致知道位置，就足够使我们作出非常准确的估计。我们知道这一点是因为新的协方差在$x$轴上并没有达到$t=1$时的观测协方差那么远。所以我们不仅对速度更确定，而且对位置也更确定，如果我们只使用位置观测而不考虑速度的话！

这是一个关键点，所以仔细阅读！雷达能够探测飞机的位置，这称为观测量。根据位置估计我们可以计算速度，我们称速度为隐藏量。隐藏意味着：没有传感器对速度进行观测，因此它的价值对我们是隐藏的。我们却能够利用位置和速度之间的相关性非常准确地推断出它的值。

当然，还有一些不可观察的变量。例如，飞机的状态包括航向、发动机转速、重量、颜色、飞行员的名字等。我们无法使用位置传感器直接感应到这些，因此无法观察到它们；同样，也无法从传感器观测值和相关性中推断出它们（比如：红色飞机的速度不会比白色飞机快），因此它们没有被隐藏。相反，它们是不可观察的。如果在滤波状态中包含一个不可观察的变量，那么对该变量的估计将是无意义的。

究竟是什么让这成为可能？想象一下，飞机的速度告诉我们一些非常重要的东西——飞行的方向和速度。只要飞机不改变速度，速度就允许我们预测下一个位置。想想看，如果你突然改变方向，你的位置也会有很大的改变。如果位置的观测不在速度变化的方向上，则不太可能是真的。这两者是相关的，所以如果速度改变，位置也必须改变，并且以一种可预测的方式。

重要的是要明白，我们是在利用速度和位置是相关的这一事实。我们从两次观测之间的距离和时间得到了速度的粗略估计，并且使用贝叶斯定理只需几次观测就可以得到非常精确的估计。如果你有任何疑问，请重新阅读。如果你不明白这一点，你会很快发现你在接下去的章节中无法推理。

更高维度

到目前为止，我已经向你们展示了二维高斯，但数学并没有把你们限制在二维。在后面的章节中，我们将研究9维，甚至12维。如果你在诸如天气预报之类的领域工作，你最终会得到上千个维度。

这些更高的维度看起来像什么？二维高斯可以用误差椭圆来表示，所以三维高斯可以用三维误差椭圆来表示。FilterPy提供了一个函数来绘制这个椭球。

首先，让我们用一个给定的协方差生成一些有噪声的数据，这样我们就可以在椭球体内部绘制它。

from filterpy.stats import plot_3d_covariancemu = [0.3, 5., 10.]
C = np.array([[1.0, .03, .2],[.03,  4.0, .0],[.2,  .0, 16.1]])sample = np.random.multivariate_normal(mu, C, size=1000)

现在我们用FilterPy的plot_3d_covariance()函数绘制椭球体，再对样本进行散点绘制：

ax = plot_3d_covariance(mu, C, alpha=.4, std=3, limit_xyz=True)
ax.scatter(sample[:, 0], sample[:, 1], zs=sample[:, 2],)

理论上说，大约99%的分布会在3个标准差内，这似乎是事实。

九个维度？我还没有办法在2D屏幕上绘制9D椭球体，所以不会有图形。但概念是一样的，分布的标准差误差可以用一个9维椭球来描述。

总结

我们利用了系统的相关性，得出了非常精确的估计。数学并不关心我们是在处理两个位置，还是一个位置和一个相关的速度，或者其他的空间维度。如果建筑面积和房价相关，你可以写一个卡尔曼滤波器来跟踪房价；如果年龄与疾病发生率相关，你可以编写一个卡尔曼滤波器来跟踪疾病；如果僵尸数量与猎枪数量成反比，那么你可以编写一个卡尔曼滤波器来跟踪僵尸数量。我向你们展示了三角定位，那只是为了建立你的直觉。你当然也可以，为无法用几何表示的状态量编写一个卡尔曼滤波器。如果我们可以将不确定性表示为多维高斯分布，我们就可以将先验值与似然值相乘，得到更精确的结果。

相关阅读

• Kalman-and-Bayesian-Filters-in-Python/05-Multivariate-Gaussians