运行(一次通过)计算协方差

Question

我有一组3d矢量(x,y,z),我想计算协方差矩阵而不存储矢量.

我将在C#中完成它,但最终我将在C上在微控制器上实现它,所以我需要算法本身,而不是库.

伪代码也很棒.

Answer 1

该公式很简单,如果你有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可以具有固定的尺寸,并且可以为此目的轻松编码.您甚至可以将代码重写为离散浮点计算,并避免计算协方差矩阵的对称部分.

请注意,第二行代码中有一个向量外积.并非所有矢量库都能正确解释它.

Answer 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);
    }
}