高斯分布中均值，方差，协方差的计算及matlab实现

原帖请见：http://blog.sina.com.cn/s/blog_4aa4593d01012am3.html

今天看论文的时候又看到了协方差矩阵这个破东西，以前看模式分类的时候就特困扰，没想到现在还是搞不清楚，索性开始查协方差矩阵的资料，恶补之后决定马上记录下来，嘿嘿~本文我将用自认为循序渐进的方式谈谈协方差矩阵。

统计学的基本概念

学过概率统计的孩子都知道，统计里最基本的概念就是样本的均值，方差，或者再加个标准差。首先我们给你一个含有n个样本的集合，依次给出这些概念的公式描述，这些高中学过数学的孩子都应该知道吧，一带而过。

均值：


标准差：


方差：

很 显然，均值描述的是样本集合的中间点，它告诉我们的信息是很有限的，而标准差给我们描述的则是样本集合的各个样本点到均值的距离之平均。以这两个集合为 例，[0，8，12，20]和[8，9，11，12]，两个集合的均值都是10，但显然两个集合差别是很大的，计算两者的标准差，前者是8.3，后者是 1.8，显然后者较为集中，故其标准差小一些，标准差描述的就是这种“散布度”。之所以除以n-1而不是除以n，是因为这样能使我们以较小的样本集更好的 逼近总体的标准差，即统计上所谓的“无偏估计”。而方差则仅仅是标准差的平方。

为什么需要协方差？

上 面几个统计量看似已经描述的差不多了，但我们应该注意到，标准差和方差一般是用来描述一维数据的，但现实生活我们常常遇到含有多维数据的数据集，最简单的 大家上学时免不了要统计多个学科的考试成绩。面对这样的数据集，我们当然可以按照每一维独立的计算其方差，但是通常我们还想了解更多，比如，一个男孩子的 猥琐程度跟他受女孩子欢迎程度是否存在一些联系啊，嘿嘿~协方差就是这样一种用来度量两个随机变量关系的统计量，我们可以仿照方差的定义：

来度量各个维度偏离其均值的程度，标准差可以这么来定义：

协方差的结果有什么意义呢？如果结果为正值，则说明两者是正相关的(从协方差可以引出“相关系数”的定义)，也就是说一个人越猥琐就越受女孩子欢迎，嘿嘿，那必须的~结果为负值就说明负相关的，越猥琐女孩子越讨厌，可能吗？如果为0，也是就是统计上说的“相互独立”。

从协方差的定义上我们也可以看出一些显而易见的性质，如：

协方差多了就是协方差矩阵

上一节提到的猥琐和受欢迎的问题是典型二维问题，而协方差也只能处理二维问题，那维数多了自然就需要计算

多个协方差，比如n维的数据集就需要计算个协方差，那自然而然的我们会想到使用矩阵来组织这些数据。给出协方差矩阵的定义：

这个定义还是很容易理解的，我们可以举一个简单的三维的例子，假设数据集有三个维度，则协方差矩阵为

可见，协方差矩阵是一个对称的矩阵，而且对角线是各个维度上的方差。

Matlab协方差实战

上面涉及的内容都比较容易，协方差矩阵似乎也很简单，但实战起来就很容易让人迷茫了。必须要明确一点，协方差矩阵计算的是不同维度之间的协方差，而不是不同样本之间的。这个我将结合下面的例子说明，以下的演示将使用Matlab，为了说明计算原理，不直接调用Matlab的cov函数(蓝色部分为Matlab代码)。

首先，随机产生一个10*3维的整数矩阵作为样本集，10为样本的个数，3为样本的维数。
MySample = fix(rand(10,3)*50)

根据公式，计算协方差需要计算均值，那是按行计算均值还是按列呢，我一开始就老是困扰这个问题。前面我们也特别强调了，协方差矩阵是计算不同维度间的协方差，要时刻牢记这一点。样本矩阵的每行是一个样本，每列为一个维度，所以我们要按列计算均值。为了描述方便，我们先将三个维度的数据分别赋值：

dim1 = MySample(:,1);
dim2 = MySample(:,2);
dim3 = MySample(:,3);

计算dim1与dim2，dim1与dim3，dim2与dim3的协方差：

sum( (dim1-mean(dim1)) .* (dim2-mean(dim2)) ) / ( size(MySample,1)-1 ) % 得到  74.5333
sum( (dim1-mean(dim1)) .* (dim3-mean(dim3)) ) / ( size(MySample,1)-1 ) % 得到  -10.0889
sum( (dim2-mean(dim2)) .* (dim3-mean(dim3)) ) / ( size(MySample,1)-1 ) % 得到  -10***000

搞清楚了这个后面就容易多了，协方差矩阵的对角线就是各个维度上的方差，下面我们依次计算：

std(dim1)^2 % 得到   108.3222
std(dim2)^2 % 得到   260.6222
std(dim3)^2 % 得到  94.1778

这样，我们就得到了计算协方差矩阵所需要的所有数据，调用Matlab自带的cov函数进行验证：

cov(MySample)

把我们计算的数据对号入座，是不是一摸一样？

Update： 今天突然发现，原来协方差矩阵还可以这样计算，先让样本矩阵中心化，即每一维度减去该维度的均值，使每一维度上的均值为0，然后直接用新的到的样本矩阵乘 上它的转置，然后除以(N-1)即可。其实这种方法也是由前面的公式通道而来，只不过理解起来不是很直观，但在抽象的公式推导时还是很常用的！同样给出 Matlab代码实现：

X = MySample – repmat(mean(MySample),10,1);    % 中心化样本矩阵，使各维度均值为0
C = (X’*X)./(size(X,1)-1)

总结

理解协方差矩阵的关键就在于牢记它计算的是不同维度之间的协方差，而不是不同样本之间，拿到一个样本矩阵，我们最先要明确的就是一行是一个样本还是一个维度，心中明确这个整个计算过程就会顺流而下，这么一来就不会迷茫了~

P.S.写论文要选Latex，在wordpress里编辑公式还得用Latex，用Latex还真对得起咱学计算机这张脸~

例2：

a =
    -1     1     2
    -2     3     1
     4     0     3

for i=1:size(a,2) 
    for j=1:size(a,2) 
        c(i,j)=sum((a(:,i)-mean(a(:,i))).*(a(:,j)-mean(a(:,j))))/(size(a,1)-1);
    end 
end

c =

   10.3333   -4.1667    3.0000

   -4.1667    2.3333   -1.5000

    3.0000   -1.5000    1.0000 

 c为求得的协方差矩阵，在matlab以矩阵a的每一列为变量，对应的每一行为样本。这样在矩阵a中就有3个列变量分别为a（:,1), a(:,2), a(:,3)。

 在协方差矩阵c中，每一个元素c(i,j)为对第i列与第j列的协方差，例如c(1,2) = -4.1667为第一列与第二列的协方差。

 

 拿c(1,2)的求解过程来说

 c(1,2)=sum((a(:,1)-mean(a(:,1))).*(a(:,2)-mean(a(:,2))))/(size(a,1)-1);

 1. a(:,1)-mean(a(:,1))，第一列的元素减去该列的均值得到

   -1.3333

   -2.3333

    3.6667

2,  a(:,2)-mean(a(:,2))，第二列的元素减去该列的均值得到

   -0.3333

    1.6667

   -1.3333

3, 再将第一步与第二部的结果相乘

   -1.3333        -0.3333           0.4444

   -2.3333  .*     1.6667  =     -3.8889

    3.6667         -1.3333          -4.8889

 

4, 再将结果求和/size(a,1)-1 得 -4.1667,该值即为c（1,2）的值。

 

再细看一下是不是与协方差公式：Cov(X,Y) = E{ [ (X-E(X) ] [ (Y-E(Y) ] } 过程基本一致呢，只是在第4步的时候matlab做了稍微的调整，自由度为n-1，减少了一行的样本值个数。