协方差和它的c/c++算法

今天写个协方差的c++算法， 在初始步骤又犯了2，彻底总结一下协方差的数学概念和它的算法。
参考文章
http://zipperary.com/2014/01/12/covariance/
http://grunt1223.iteye.com/blog/961063
数学理论方面的都是copy第一篇的，我认为他/她已经写得非常好了。
在机器学习中经常需要计算协方差矩阵，本科时没学过这个概念，一直对此非常头疼。现在试图通过实例的计算、图形化的表示来梳理一下什么是协方差矩阵。
A numerical example
问题：
有一组数据（如下），分别为二维向量，这四个数据对应的协方差矩阵是多少？

解答：
由于数据是二维的，所以协方差矩阵是一个2*2的矩阵，矩阵的每个元素为：
元素(i,j) = (第 i 维所有元素 - 第 i 维的均值) * (第 j 维所有元素 - 第 j 维的均值) 。
其中「*」代表向量内积符号，即两个向量求内积，对应元素相乘之后再累加。
我们首先列出第一维：
D1: (1,3,4,5) 均值：3.25
D2: (2,6,2,2) 均值：3
下面计算协方差矩阵第(1,2)个元素：
元素(1,2)=(1-3.25,3-3.25,4-3.25,5-3.25)*(2-3,6-3,2-3,2-3)=-1
类似的，我们可以把2*2个元素都计算出来：

这个题目的最终结果就是：

An explanation
我们来分析一下上面的例子。首先看一下元素(1,1)的计算过程：
把所有数据的第一个维度拿出来，求出均值，之后的求解过程完全是我们熟悉的「方差」的求法。也就是说，这完完全全就是在求所有数据第一维元素（共4个）的方差（8.75）嘛。类似地，元素(2,2)求的是第二维(共4个)元素的方差（12）。
再来看元素(1,2)，这分明就是我们高数里面学的求 x 和 y 的协方差，不再单独计算某一维度的分散程度，而是把两个维度的分散值结合起来，这里才真正体现了「协方差矩阵」中「协方差」的意味。从计算过程和计算结果都能看出，元素(2,1)与元素(1,2)是一样的。也就是说，所有协方差矩阵都是一个对称阵。
总结一下协方差矩阵的特点：
对角线元素(i,i)为数据第 i 维的方差。
非对角线元素(i,j)为第 i 维和第 j 维的协方差。
协方差矩阵是对称阵。
现在只需要了解这些就够了。
//——————————-理论分割线—————————
那么如何编写程序呢，首先推导一下
由于协方差矩阵中的每个元素的计算方法都是一样，我们只举一例就可以了。 例如计算第一行第一列的元素，基本公式如下

X11=∑1N(xi1−x¯1)2

这里，上标i代表第i个数据， 下标1代表第1维， 总的意思就是第i个数据的第一维。后面的代表第一维数据的平均值。 这里也可以看出 对角元素就是对应维度分量的方差。
继续对公式变换，为了简便易看， 由于只讨论第一维数据， 所以以下的公式就不再特意写明， 下标i代表的就是第i个数据， 平均值代表的就是第一维数据的平均值，这点请注意。
对上面的公式展开

X11=∑1Nx2i−2∑1nxix¯+∑1nx¯2

中间一项拿出来推导

∑1nxix¯=2x¯∑1nxi=2x¯×n×x¯=2∑1nx¯2

然后带回到第二个公式，看看结果

X11=∑1Nx2i−∑1nx¯2

如果按照最基本的公式计算， 那么需要O(n)去计算平均值，然后再一个O(n)计算协方差。
推导到这里就可以发现，只需要一个循环就可以了， 计算均值的同时可以得到∑n1x2,最后减掉n倍的x¯2就可以了。

如果是二维数据的话，代码大概是这样子的

 covMat11 = covMat12 = covMat21 = covMat22 = 0;     for(i=0; i<dataSize; i++) {          meanX +=data[i]->x;          meanY +=data[i]->y;          covMat11     +=data[i]->x * data[i]->x;          covMat12     +=data[i]->x * data[i]->y;          covMat22     +=data[i]->y * data[i]->y;     }     meanX/=dataSize;     meanY/=dataSize;     covMat11 -= dataSize*meanX*meanX;        covMat12 -= dataSize*meanX*meanY;     covMat22 -= dataSize*meanY*meanY;     covMat21 = covMat12;

其他的元素计算方法类似。