运行（一次）协方差计算

Question

运行（一次）协方差计算

c#statisticscovariance

7

我有一组三维向量 (x,y,z)，想要计算协方差矩阵，但不想存储这些向量。

我将在 C# 中实现该功能，但最终需要在微控制器上的 C 语言中实现，因此我需要算法本身而不是库。

提供伪代码将会很好。

- Terje Kolderup

4个回答

2

< p > emu 的代码很优雅，但需要一个额外的步骤才能正确：< /p>

Vector mean;
Matrix covariance;
for (int i = 0; i < points.size(); ++i) {
  Vector diff = points[i] - mean;
  mean += diff / (i + 1);
  covariance += diff * diff.transpose() * i / (i + 1);
}

covariance = covariance/(points.size()-1);

注意协方差归一化的最后一步。

- Menos

我的公式中的协方差矩阵已经正确地针对一个精确的数据集进行了归一化处理。你所说的是“样本协方差”，它是从有限数量的测量值中进行统计估计得出的。然后，你最后一行应该写成 covariance = covariance * (points.size() / (points.size() - 1))。 - emu

2

我认为我已经找到了解决方案。它基于这篇关于如何手动计算协方差的文章和这篇关于计算运行方差的文章。接着，我根据我从第一篇文章中的理解，修改了后者中的算法来计算协方差而不是方差。

public class CovarianceMatrix
{
    private int _n;
    private Vector _oldMean, _newMean,
                    _oldVarianceSum, _newVarianceSum,
                    _oldCovarianceSum, _newCovarianceSum;

    public void Push(Vector x)
    {
        _n++;
        if (_n == 1)
        {
            _oldMean = _newMean = x;
            _oldVarianceSum = new Vector(0, 0, 0);
            _oldCovarianceSum = new Vector(0, 0, 0);
        }
        else
        {
            //_newM = _oldM + (x - _oldM) / _n;
            _newMean = new Vector(
                _oldMean.X + (x.X - _oldMean.X) / _n,
                _oldMean.Y + (x.Y - _oldMean.Y) / _n,
                _oldMean.Z + (x.Z - _oldMean.Z) / _n);

            //_newS = _oldS + (x - _oldM) * (x - _newM);
            _newVarianceSum = new Vector(
                _oldVarianceSum.X + (x.X - _oldMean.X) * (x.X - _newMean.X),
                _oldVarianceSum.Y + (x.Y - _oldMean.Y) * (x.Y - _newMean.Y),
                _oldVarianceSum.Z + (x.Z - _oldMean.Z) * (x.Z - _newMean.Z));

            /* .X is X vs Y
             * .Y is Y vs Z
             * .Z is Z vs X
             */
            _newCovarianceSum = new Vector(
                _oldCovarianceSum.X + (x.X - _oldMean.X) * (x.Y - _newMean.Y),
                _oldCovarianceSum.Y + (x.Y - _oldMean.Y) * (x.Z - _newMean.Z),
                _oldCovarianceSum.Z + (x.Z - _oldMean.Z) * (x.X - _newMean.X));

            // set up for next iteration
            _oldMean = _newMean;
            _oldVarianceSum = _newVarianceSum;
        }
    }
    public int NumDataValues()
    {
        return _n;
    }

    public Vector Mean()
    {
        return (_n > 0) ? _newMean : new Vector(0, 0, 0);
    }

    public Vector Variance()
    {
        return _n <= 1 ? new Vector(0, 0, 0) : _newVarianceSum.DivideBy(_n - 1);
    }
}

- Terje Kolderup

是的，这个方法可行！如果你想移动窗口并删除当前窗口中的值，你只需要从旧的协方差总和中减去相应的值。例如，对于x，你首先计算出“oldMean”作为tempMean.X = ((oldMean.x * n) - x) / (n-1)，然后使用oldCov.x -= (y - oldMean.y) * (x - tempMean.x)更新协方差总和，最后将你的oldMean更新为tempMean即可。 - optional

0

以下是一个在 R 中演示原理的简单示例：

a <- matrix(rnorm(22), ncol = 2)
a1 <- a[1:10, ]
a2 <- a[2:11, ]
cov(a1)
cov(a2)
m <- 10

# initial step
m1.1 <- mean(a1[, 1]) 
m1.2 <- mean(a1[, 2]) 

c1.11 <- cov(a1)[1, 1]
c1.22 <- cov(a1)[2, 2]
c1.12 <- cov(a1)[1, 2]


#step 1->2
m2.1 <- m1.1 + (a[11, 1] - a[1, 1])/m
m2.2 <- m1.2 + (a[11, 2] - a[1, 2])/m

c2.11 <- c1.11 + (a[11, 1]^2 - a[1, 1]^2)/(m - 1) + (m1.1^2 - m2.1^2) * m/(m - 1)
c2.22 <- c1.22 + (a[11, 2]^2 - a[1, 2]^2)/(m - 1) + (m1.2^2 - m2.2^2) * m/(m - 1)
c2.12 <- c1.12 + (a[11, 1] * a[11, 2] - a[1, 1]*a[1, 2])/(m - 1) + 
   (m1.1 * m1.2 - m2.1 * m2.2) * m/(m - 1)

cov(a2) - matrix(c(c2.11, c2.12, c2.12, c2.22), ncol=2)

- maxatSOflow

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- emu · Accepted Answer

如果你手头有Matrix和Vector类，那么这个公式就很简单：

Vector mean;
Matrix covariance;
for (int i = 0; i < points.size(); ++i) {
  Vector diff = points[i] - mean;
  mean += diff / (i + 1);
  covariance += diff * diff.transpose() * i / (i + 1);
}
covariance *= 1 / points.size()

相较于双重计算，我个人更喜欢这种方式。代码简短，结果完美无缺。

Matrix和Vector可以有固定尺寸并且容易为此目的编写代码。你甚至可以将代码重写成离散浮点计算，避免计算协方差矩阵的对称部分。

请注意，在代码的倒数第二行有一个向量的外积。并非所有向量库都能正确解释它。