如何在 MATLAB 中没有任何内置函数或循环的情况下计算协方差矩阵？

Question

是否可以在不使用 MATLAB 中的任何内置函数或循环的情况下找到矩阵的协方差？我完全不知道解决这个问题的想法。

我在想这样的事情：

cov(x,y) = 1/(n-1) .* (x*y)

但是，我认为这行不通。有什么想法吗？

Answer 1

这里有一个很好的例子来说明如何数值计算协方差矩阵。 http://www.itl.nist.gov/div898/handbook/pmc/section5/pmc541.htm。但是，为了完整起见，让我们把它放在这篇文章中。 我对你所说的“内置”函数有点困惑，因为协方差要求你对矩阵的列求和。如果你不能使用任何内置函数来总结这些元素，那么如果不使用for 循环，我看不到你怎么能做到这一点。 编辑：我想出了如何在不使用内置函数或循环的情况下做到这一点，但您需要使用size 来确定矩阵中有多少行......除非您在函数中将其指定为常量。

在数值上，您可以像这样计算协方差矩阵：

基本上，协方差矩阵的第 i^th 行和第 j^th 列是这样的，即您将 i 列的乘积之和减去平均值列i 与列j 减去列j 的平均值。现在，把这些加起来，然后除以n - 1。这被称为unbiased estimator。您还会注意到该矩阵是对称的，因为即使您翻转顺序（即查看j 列，然后查看i 列），答案应该仍然是相同的。我假设你也不能使用 MATLAB 中的mean，所以让我们从第一原则开始吧。

首先，计算一个行向量，计算每列的平均值。您可以在不使用sum 的情况下计算所有列的总和，因为它也是一个内置函数，将这个 1 的行向量与您的矩阵 A 相乘，输出将是一个行向量包含所有列的总和。因此，请执行以下操作：

one_vector(1:size(A,1)) = 1;
mu = (one_vector * A) / size(A,1);

第一行代码的诀窍在于，我们动态创建了一个与矩阵A 中的行数长度相同的数组。我们把这个完全填满 1。请注意，您可以使用ones，但您说您不能使用任何内置函数。 mu 将包含我们所有列的向量。

现在，让我们通过用平均值减去每一列来预处理数据，因为这就是定义所说的我们要做的。要在没有任何内置函数的情况下执行此操作，您可以做的是用各自的方法减去所有列，重复 mu 的次数与 one_vector 中的 1 一样多。因此：

A_mean_subtract = A - mu(one_vector, :);

这里有点棘手（而且很酷）。如果我们转置矩阵A，你会看到行变成了列，列变成了行。如果我们采用这个转置并乘以原始矩阵，我们实际上会得到矩阵A 的列i 和列j 之间的乘积之和。这是我们协方差计算的第一部分。然后我们除以n - 1。因此，我们的协方差很简单：

covA = (A_mean_subtract.' * A_mean_subtract) / (size(A,1) - 1);

这是一个简单的示例，以及我在上面向您展示的那个网站上看到的内容。假设A 是这样的：

A = [4 2 0.5; 4.2 2.1 0.59; 3.9 2.0 0.58; 4.3 2.1 0.62; 4.1 2.2 0.63]

A =

    4.0000    2.0000    0.5000
    4.2000    2.1000    0.5900
    3.9000    2.0000    0.5800
    4.3000    2.1000    0.6200
    4.1000    2.2000    0.6300

运行上面的代码，我们得到的是这样的：

covA =

    0.0250    0.0075    0.0042
    0.0075    0.0070    0.0034
    0.0042    0.0034    0.0026

您会看到这也与 MATLAB 中的 cov 函数匹配：

>> cov(A)

ans =

0.0250    0.0075    0.0042
0.0075    0.0070    0.0034
0.0042    0.0034    0.0026

一点提示

如果您在 MATLAB 命令提示符中输入 edit cov，您实际上可以看到它们如何在没有任何 for 循环的情况下计算协方差矩阵...。这与我给您的答案基本相同 :)

如果您想更有效地做到这一点

假设您可以使用sum 和bsxfun，我们可以用更少（并且更有效..）的代码行来做到这一点。首先，像我们上面使用sum 一样计算你的平均向量：

mu = sum(A) / size(A,1);

现在，要将矩阵A 与每列对应的平均值相减，您可以使用bsxfun 来帮助您进行减法：

A_mean_subtract = bsxfun(@minus, A, mu);

现在，像以前一样计算协方差矩阵：

covA = (A_mean_subtract.' * A_mean_subtract) / (size(A,1) - 1);

您应该得到与我们之前看到的完全相同的结果。

关于稳定性的小提示

我们使用直接定义来计算两列之间的协方差。但是，已经表明，如果您提供某些类型的数据，则使用直接定义可能会导致数值不稳定。请参阅this Wikipedia page，它通过各种算法计算两个更稳定的n 长度向量之间的协方差。