协方差矩阵的定义及其计算公式 （转）

协方差矩阵的定义及其计算公式

 (2014-03-20 18:36:56)

转载▼

 分类： matlab

协方差矩阵在机器学习中经常用到，查看wiki：http://zh.wikipedia.org/wiki/协方差矩阵 可知协方差矩阵的具体计算公式如下：

在统计学与概率论中，协方差矩阵是一个矩阵，其每个元素是各个向量元素之间的协方差。这是从标量随机变量到高维度随机向量的自然推广。

假设是以个标量随机变量组成的列向量，

并且是其第i个元素的期望值，即, 。协方差矩阵被定义的第i，j项是如下：

即：

矩阵中的第个元素是与的协方差。这个概念是对于标量随机变量方差的一般化推广。

 

二、理解的关键

　　1、理解的关键是两个随机变量x₁,x₂的协方差如何计算，有cov(x1,x2) = E{[(x₁-E(x₁)][x₂-E(x₂)]}，对于离散的随机变量x₁,x₂,协方差矩阵描述的x₁,x₂相互联系的偏差，所以两个随机变量是一一对应的，即假设有m个样本值，则分别为(x₁₁,x₂₁),(x₁₂,x₂₂),(x₁₃,x₂₃),...(x_1m,x_2m)，这便可以写成一个2*m的矩阵的形式。则x₁,x₂协方差表示的是两个随机变量对应的样本值分别都减去各自均值后的乘积的均值（因为无偏性估计的缘故，除以的不是m而是m-1)；

　　2、所以对于一个n*m的样本矩阵，得出的协方差矩阵C是n*n的矩阵，协方差矩阵每个元素C_ij表示的随机变量x_i,x_j的协方差。所以协方差矩阵是一个对称矩阵，且对角线上元素为每个随机变量的方差（如果是信号的话可以看成是能量）；如果各个变量相关性很小的话，互相的协方差接近0，即协方差矩阵基本上为对角阵；

　　3、可以证明，协方差矩阵是非负定矩阵，这可以有非负定矩阵的定义得到；

　　4、同样地，为了表示各个随机变量相关性到底有多大，可以引入相关性矩阵。

统计学的基本概念

学过概率统计的孩子都知道，统计里最基本的概念就是样本的均值，方差，或者再加个标准差。首先我们给你一个含有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为样本的维数。

1MySample = fix(rand(10,3)*50)

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

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

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

 1sum( (dim1-mean(dim1)) .* (dim2-mean(dim2)) ) / ( size(MySample,1)-1 ) %(10行；（10-1）) 得到 74.53332；sum( (dim1-mean(dim1)) .* (dim3-mean(dim3)) ) / ( size(MySample,1)-1 ) % 得到  -10.08893；sum( (dim2-mean(dim2)) .* (dim3-mean(dim3)) ) / ( size(MySample,1)-1 ) % 得到  -106.4000

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

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

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

1cov(MySample)

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

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

1X = MySample - repmat(mean(MySample),10,1);    % 中心化样本矩阵，使各维度均值为02C = (X'*X)./(size(X,1)-1);%B = repmat(A,m,n)将矩阵 A 复制 m×n 块，即把 A 作为 B 的元素，B 由 m×n 个 A 平铺而成。B 的维数是 [size(A,1)*m, (size(A,2)*n]

>> A = [1,2;3,4]

A =

1 2

3 4

>> B = repmat(A,2,3)

B =

1 2 1 2 1 2

3 4 3 4 3 4

1 2 1 2 1 2

3 4 3 4 3 4

 

总结

理解协方差矩阵的关键就在于牢记它计算的是不同维度之间的协方差，而不是不同样本之间，拿到一个样本矩阵，我们最先要明确的就是一行是一个样本还是一个维度。

matlab：

1、均值
数学定义：
Matlab 函数：mean
>>X=[1,2,3]
>>mean(X)=2
如果X 是一个矩阵，则其均值是一个向量组。mean(X,1)为列向量的均值，mean(X,2)为行向量的均
值。
>>X=[1 2 3
         4 5 6]
>>mean(X,1)=[2.5, 3.5, 4.5]
>>mean(X,2)=[2
                       5]
若要求整个矩阵元素的均值，则为mean(mean(X))。
>>mean(mean(X))=3.5
也可使用mean2函数：
>>mean2(X)=3.5
median，求一组数据的中值，用法与mean 相同。
>>X=[1,2,9]
>>mean(X)=4
>>median(X)=2
2、方差
数学定义：
均方差：
Matlab 函数：var
要注意的是var 函数所采用公式中，分母不是，而是。这是因为var 函数实际上求的并不是方差，
而是误差理论中“有限次测量数据的标准偏差的估计值”。
>>X=[1,2,3,4]
>>var(X)=1.6667
>> sum((X(1,:)-mean(X)).^2)/length(X)=1.2500
>> sum((X(1,:)-mean(X)).^2)/(length(X)-1)=1.6667
var 没有求矩阵的方差功能，可使用std 先求均方差，再平方得到方差。
std，均方差，std(X,0,1)求列向量方差，std(X,0,2)求行向量方差。
>>X=[1 2
         3 4]
>>std(X,0,1)=1.4142 1.4142
>>std(X,0,2)=0.7071
                     0.7071
若要求整个矩阵所有元素的均方差，则要使用std2函数：
>>std2(X)=1.2910
4、协方差矩阵
A=[61.45,55.9,61.95,59,58.14,53.61,55.48,54.21,61.52,54.92];
B=[40.36,39.8,49.2,48,51.5,49.39,51.13,58.06,61,62.35];
C=[8.61,8.91,10.43,13.32,13.48,15.75,18.14,19.95,21.95,23.53];
D=[14.31,14.72,15.28,15.91,14.67,15,15.86,15.16,13.72,12.94];
E=[7.67,7.75,8.15,9.24,10.68,10.58,10.31,10,8.91,8.51];
>> q=[A',B',C',D',E'];
>> w=cov(q)
w =
10.3710 -4.7446 -6.6023 -0.1873 -1.8881
-4.7446 59.1503 38.7606 -3.0743 3.0982
-6.6023 38.7606 28.6966 -2.0199 2.4166
-0.1873 -3.0743 -2.0199 0.8474 0.3936

-1.8881 3.0982 2.4166 0.3936 1.3412

例：

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

for i=1:size(a,2) %size(a,1)求矩阵的行数；size(a,2)求矩阵的列数，相当于length(a)；size(a)同时求矩阵的行和列数
    for j=1:size(a,2) %A(1,:)表示这个矩阵的第一行；A(:,1)表示这个矩阵的第一列

        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，减少了一行的样本值个数。

clear all;
a =[-1,1,2;-2,3,1;4,0,3]
b=cov(a)
for i=1:size(a,2) %size(a,1)求矩阵的行数；size(a,2)求矩阵的列数
    for j=1:size(a,2) %A(1,:)表示这个矩阵的第一行；A(:,1)表示这个矩阵的第一列
    c(i,j)=sum((a(:,i)-mean(a(:,i))).*(a(:,j)-mean(a(:,j))))/(size(a,1)-1)
    end
end