2016-10-28 14:37:44 shenziheng1 阅读数 35399

1.前言

看论文的时候又看到了协方差矩阵这个破东西,以前看图像处理的书籍的时候就特困扰,没想到现在还是搞不清楚,索性开始查协方差矩阵的资料,恶补之后决定马上记录下来。

2.拼出身—统计学的定义

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

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

那么问题来了,上面介绍的参量难道不能描述统计学所有特性吗?为啥还要搞出一个协方差出来

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

度量各个维度偏离其均值的程度,协方差可以这么来定义:

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

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



3.协方差矩阵的由来

好几十年前,鲁迅爷爷就说过,世界上本没有路,走的人多了也就有了路。协方差矩阵也是这样,好多协方差凑合到一起就形成了协方差矩阵。当然,数学的定义,不能如我这样随意。对于一个二维矩阵,每一个因子都可以视为两个不同随机变量的关系,这正好和协方差矩阵多少有点牵连,因此数学家们就把协方差矩阵引入到了二维矩阵中,衡量各个变量之间的紧密程度(就是关系度啦)。根据协方差的性质,我们可以类似的推出协方差矩阵的性质:


1.协方差矩阵一定是个对称的方阵

2.协方差矩阵对角线上的因子其实就是变量的方差:cov(X,X)=var(X)


这个定义还是很容易理解的,我们可以举一个简单的三变量的例子,假设数据集有{x,y,z}{x,y,z}三个维度,则协方差矩阵为:

再一次可以看出,协方差矩阵是一个对称的矩阵,而且对角线是各个变量上的方差。

3.MATLAB实战练习

上面涉及的内容都比较容易,协方差矩阵似乎也很简单,但实战起来就很容易让人迷茫了。必须要明确一点,### 协方差矩阵计算的是不同维度之间的协方差,而不是不同样本之间的
为了说明计算原理,不直接调用Matlab的cov函数:
首先,随机产生一个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 ) % 得到  -106.4000
搞清楚了这个后面就容易多了,协方差矩阵的对角线就是各个维度上的方差,下面我们依次计算:
std(dim1)^2 % 得到   108.3222
std(dim2)^2 % 得到   260.6222
std(dim3)^2 % 得到   94.1778
这样,我们就得到了计算协方差矩阵所需要的所有数据:

C11=108.3222 C12=74.5333 C13=-10.0889
C21=74.5333 C22=260.6222 C23=-106.4000
C31=-10.0889 C32=-106.4000 C33=94.1778


调用Matlab自带的cov函数进行验证:


4.心得感悟

理解协方差矩阵的关键就在于牢记它计算的是不同维度之间的协方差,而不是不同样本之间,拿到一个样本矩阵,我们最先要明确的就是一行是一个样本还是一个维度,心中明确这个整个计算过程就会顺流而下,这么一来就不会迷茫了。
其实还有一个更简单的容易记还不容易出错的方法:协方差矩阵一定是一个对称的方阵,一定是一个对称的方阵一定是一个对称的方阵!!!记住就好啦~
2016-06-19 16:30:21 cwcww1314 阅读数 2986

http://pinkyjie.com/2011/02/24/covariance-pca/


自从上次谈了协方差矩阵之后,感觉写这种科普性文章还不错,那我就再谈一把协方差矩阵吧。上次那篇文章在理论层次介绍了下协方差矩阵,没准很多人觉得这东西用处不大,其实协方差矩阵在好多学科里都有很重要的作用,比如多维的正态分布,再比如今天我们今天的主角——主成分分析(Principal Component Analysis,简称PCA)。结合PCA相信能对协方差矩阵有个更深入的认识~

PCA的缘起

PCA大概是198x年提出来的吧,简单的说,它是一种通用的降维工具。在我们处理高维数据的时候,为了能降低后续计算的复杂度,在“预处理”阶段通常要先对原始数据进行降维,而PCA就是干这个事的。本质上讲,PCA就是将高维的数据通过线性变换投影到低维空间上去,但这个投影可不是随便投投,要遵循一个指导思想,那就是:找出最能够代表原始数据的投影方法。这里怎么理解这个思想呢?“最能代表原始数据”希望降维后的数据不能失真,也就是说,被PCA降掉的那些维度只能是那些噪声或是冗余的数据。这里的噪声和冗余我认为可以这样认识:

  • 噪声:我们常说“噪音污染”,意思就是“噪声”干扰我们想听到的真正声音。同样,假设样本中某个主要的维度A,它能代表原始数据,是“我们真正想听到的东西”,它本身含有的“能量”(即该维度的方差,为啥?别急,后文该解释的时候就有啦~)本来应该是很大的,但由于它与其他维度有那么一些千丝万缕的相关性,受到这些个相关维度的干扰,它的能量被削弱了,我们就希望通过PCA处理后,使维度A与其他维度的相关性尽可能减弱,进而恢复维度A应有的能量,让我们“听的更清楚”!
  • 冗余:冗余也就是多余的意思,就是有它没它都一样,放着就是占地方。同样,假如样本中有些个维度,在所有的样本上变化不明显(极端情况:在所有的样本中该维度都等于同一个数),也就是说该维度上的方差接近于零,那么显然它对区分不同的样本丝毫起不到任何作用,这个维度即是冗余的,有它没它一个样,所以PCA应该去掉这些维度。

这么一分析,那么PCA的最终目的就是“降噪”和消灭这些“冗余”的维度,以使降低维度的同时保存数据原有的特征不失真。后面我们将结合例子继续讨论。

协方差矩阵——PCA实现的关键

前面我们说了,PCA的目的就是“降噪”和“去冗余”。“降噪”的目的就是使保留下来的维度间的相关性尽可能小,而“去冗余”的目的就是使保留下来的维度含有的“能量”即方差尽可能大。那首先的首先,我们得需要知道各维度间的相关性以及个维度上的方差啊!那有什么数据结构能同时表现不同维度间的相关性以及各个维度上的方差呢?自然是非协方差矩阵莫属。回忆下 浅谈协方差矩阵的内容,协方差矩阵度量的是维度与维度之间的关系,而非样本与样本之间。协方差矩阵的主对角线上的元素是各个维度上的方差(即能量),其他元素是两两维度间的协方差(即相关性)。我们要的东西协方差矩阵都有了,先来看“降噪”,让保留下的不同维度间的相关性尽可能小,也就是说让协方差矩阵中非对角线元素都基本为零。达到这个目的的方式自然不用说,线代中讲的很明确——矩阵对角化。而对角化后得到的矩阵,其对角线上是协方差矩阵的特征值,它还有两个身份:首先,它还是各个维度上的新方差;其次,它是各个维度本身应该拥有的能量(能量的概念伴随特征值而来)。这也就是我们为何在前面称“方差”为“能量”的原因。也许第二点可能存在疑问,但我们应该注意到这个事实,通过对角化后,剩余维度间的相关性已经减到最弱,已经不会再受“噪声”的影响了,故此时拥有的能量应该比先前大了。看完了“降噪”,我们的“去冗余”还没完呢。对角化后的协方差矩阵,对角线上较小的新方差对应的就是那些该去掉的维度。所以我们只取那些含有较大能量(特征值)的维度,其余的就舍掉即可。PCA的本质其实就是对角化协方差矩阵。

下面就让我们跟着上面的感觉来推推公式吧。假设我们有一个样本集X,里面有N个样本,每个样本的维度为d。即:

X={X1,,XN}Xi=(xi1,,xid)Rd,i=1,,NX={X1,…,XN}Xi=(xi1,…,xid)∈Rd,i=1,…,N

将这些样本组织成样本矩阵的形式,即每行为一个样本,每一列为一个维度,得到样本矩阵S:SRN×dS∈RN×d。我们先将样本进行中心化,即保证每个维度的均值为零,只需让矩阵的每一列除以减去对应的均值即可。很多算法都会先将样本中心化,以保证所有维度上的偏移都是以零为基点的。然后,对样本矩阵计算其协方差矩阵,按照《浅谈协方差矩阵》里末尾的update,我们知道,协方差矩阵可以简单的按下式计算得到:

C=STSN1CRd×dC=STSN−1C∈Rd×d

下面,根据我们上文的推理,将协方差矩阵C对角化。注意到,这里的矩阵C是是对称矩阵,对称矩阵对角化就是找到一个正交矩阵P,满足:PTCP=ΛPTCP=Λ。具体操作是:先对C进行特征值分解,得到特征值矩阵(对角阵)即为ΛΛ,得到特征向量矩阵并正交化即为PP。显然,P,ΛRd×dP,Λ∈Rd×d。假如我们取最大的前p(p<d)个特征值对应的维度,那么这个p个特征值组成了新的对角阵Λ1Rp×pΛ1∈Rp×p,对应的p个特征向量组成了新的特征向量矩阵P1Rd×pP1∈Rd×p

实际上,这个新的特征向量矩阵P1P1就是投影矩阵,为什么这么说呢?假设PCA降维后的样本矩阵为S1S1,显然,根据PCA的目的,S1S1中的各个维度间的协方差基本为零,也就是说,S1S1的协方差矩阵应该为Λ1Λ1。即满足:

ST1S1N1=Λ1S1TS1N−1=Λ1

而我们又有公式:

PTCP=ΛPT1CP1=Λ1PTCP=Λ⇒P1TCP1=Λ1

代入可得:

ST1S1N1=Λ1=PT1CP1=PT1(STSN1)P1=(SP1)T(SP1)N1S1TS1N−1=Λ1=P1TCP1=P1T(STSN−1)P1=(SP1)T(SP1)N−1

S1=SP1S1RN×p⇒S1=SP1S1∈RN×p

由于样本矩阵SN×dSN×d的每一行是一个样本,特征向量矩阵P1(d×p)P1(d×p)的每一列是一个特征向量。右乘P1P1相当于每个样本以P1P1的特征向量为基进行线性变换,得到的新样本矩阵S1RN×pS1∈RN×p中每个样本的维数变为了p,完成了降维操作。实际上,P1P1中的特征向量就是低维空间新的坐标系,称之为“主成分”。这就是“主成分分析”的名称由来。同时,S1S1的协方差矩阵Λ1Λ1为近对角阵,说明不同维度间已经基本独立,噪声和冗余的数据已经不见了。至此,整个PCA的过程已经结束,小小总结一下:

  1. 形成样本矩阵,样本中心化
  1. 计算样本矩阵的协方差矩阵
  1. 对协方差矩阵进行特征值分解,选取最大的p个特征值对应的特征向量组成投影矩阵
  1. 对原始样本矩阵进行投影,得到降维后的新样本矩阵

Matlab中PCA实战

首先,随机产生一个10*3维的整数矩阵作为样本集,10为样本的个数,3为样本的维数。

1
S = fix(rand(10,3)*50);

计算协方差矩阵:

1
2
3
4
S = S - repmat(mean(S),10,1);
C = (S'*S)./(size(S,1)-1);
or
C = cov(S);

对协方差矩阵进行特征值分解:

1
[P,Lambda] = eig(C);

这里由于三个方差没有明显特别小的,所以我们都保留下来,虽然维度没有降,但观察Lambda(即PCA后的样本协方差矩阵)和C(即原始的协方差矩阵),可以发现,3个维度上的方差都有增大,也就是能量都比原来增大了,3个维度上的方差有所变化,但对角线之和没有变,能量重新得到了分配,这就是“降噪”的功劳。最后我们得到降维后的样本矩阵:

1
S1 = S*P;

为了验证,我们调用matlab自带的主成分分析函数princomp

1
[COEFF,SCORE] = princomp(S) % COEFF表示投影矩阵,SCORE表示投影后新样本矩阵

对比,可以发现,SCORE和S1S1在不考虑维度顺序和正负的情况下是完全吻合的,之所以我们计算的S1S1的维度顺序不同,是因为通常都是将投影矩阵P按能量(特征值)的降序排列的,而刚才我们用eig函数得到的结果是升序。另外,在通常的应用中,我们一般是不使用matlab的princomp函数的,因为它不能真正的降维(不提供相关参数,还是我没发现?)。一般情况下,我们都是按照协方差矩阵分解后特征值所包含的能量来算的,比如取90%的能量,那就从最大的特征值开始加,一直到部分和占特征值总和的90%为止,此时部分和含有的特征值个数即为p。

经过了一番推公式加敲代码的过程,相信大家对主成分分析应该不陌生了吧,同时对协方差矩阵也有了更深层次的认识了吧,它可不只是花花枪啊。我个人觉得PCA在数学上的理论还是很完备的,想必这也是它能在多种应用中博得鳌头的原因吧。


2014-08-31 12:30:00 weixin_34375054 阅读数 192

协方差的定义

 

 

对于一般的分布,直接代入E(X)之类的就可以计算出来了,但真给你一个具体数值的分布,要计算协方差矩阵,根据这个公式来计算,还真不容易反应过来。网上值得参考的资料也不多,这里用一个例子说明协方差矩阵是怎么计算出来的吧。

记住,X、Y是一个列向量,它表示了每种情况下每个样本可能出现的数。比如给定

 

则X表示x轴可能出现的数,Y表示y轴可能出现的。注意这里是关键,给定了4个样本,每个样本都是二维的,所以只可能有X和Y两种维度。所以

 

 

 

 

用中文来描述,就是:

协方差(i,j)=(第i列的所有元素-第i列的均值)*(第j列的所有元素-第j列的均值)

这里只有X,Y两列,所以得到的协方差矩阵是2x2的矩阵,下面分别求出每一个元素:

 

       所以,按照定义,给定的4个二维样本的协方差矩阵为:

 

 

    

用matlab计算这个例子

z=[1,2;3,6;4,2;5,2]

cov(z)

ans =

    2.9167   -0.3333

   -0.3333    4.0000

可以看出,matlab计算协方差过程中还将元素统一缩小了3倍。所以,协方差的matlab计算公式为:

    协方差(i,j)=(第i列所有元素-第i列均值)*(第j列所有元素-第j列均值)/(样本数-1

       下面在给出一个4维3样本的实例,注意4维样本与符号X,Y就没有关系了,X,Y表示两维的,4维就直接套用计算公式,不用X,Y那么具有迷惑性的表达了。

 

 

 

 

 

 

 

 

    

                

        (3)与matlab计算验证

                     Z=[1 2 3 4;3 4 1 2;2 3 1 4]

                     cov(Z)

                     ans =

                          1.0000    1.0000   -1.0000   -1.0000

                          1.0000    1.0000   -1.0000   -1.0000

                         -1.0000   -1.0000    1.3333    0.6667

                          -1.0000   -1.0000    0.6667    1.3333

       可知该计算方法是正确的。我们还可以看出,协方差矩阵都是方阵,它的维度与样本维度有关(相等)。参考2中还给出了计算协方差矩阵的源代码,非常简洁易懂,在此感谢一下!

 

参考:

[1] http://en.wikipedia.org/wiki/Covariance_matrix

[2] http://www.cnblogs.com/cvlabs/archive/2010/05/08/1730319.html

 

http://blog.csdn.net/ybdesire/article/details/6270328

2018-10-22 21:29:22 u010420283 阅读数 474

0. 前言

MATLAB或者OpenCV里有很多封装好的函数,我们可以使用一行代码直接调用并得到处理结果。然而当问到具体是怎么实现的时候,却总是一脸懵逼,答不上来。前两天参加一个算法工程师的笔试题,其中就考到了这几点,感到非常汗颜!赶紧补习!

1. 双线性插值

在图像处理中,我们有时需要改变图像的尺寸,放大或者缩小。线性插值则是这类操作的关键算法。不管是放大还是缩小操作,其实都是一个像素映射的处理。如下图从小图到大图的映射,以及从大图到小图的映射。

图像来源:https://www.cnblogs.com/sdxk/p/4056223.html

然而,这两种操作都有一定的缺点。对于把小图放大的操作,因为小图中的像素点到大图中的像素点不是满射,因此大图中的点不能完全有像素值;对于将大图缩小的操作,大图中的点逆映射为小图中的点时,得到的像素坐标值可能不是整数。一种解决办法是采用最近邻方法,即将得到的坐标值与相邻的原图像中的像素坐标值比较,取离得最近的坐标值对应的像素值作为缩放后的图像对应的坐标值的像素值,但是这种办法可能导致图像失真,因此可以采用双线性差值的办法来进行计算相应的像素值。

 

对于图中红色的四个点(Q11,Q12,Q21,Q22)为源图像中存在的点,需要求在目标图像的插值(绿色点P)的坐标对应的像素值。

首先在X轴进行插值,R1,R2是两个插值过程中过渡的点.

然后在 y 方向进行线性插值,得到:

 

这样就得到所要的结果 f \left( x, y \right),

讲一个具体的例子:

比如源图像是尺寸是(100,150),现在要缩小尺寸0.6倍,即目标图像的尺寸是(60,90),则求目标图像在坐标为P[10,4]的点的像素值怎么求呢?

假设源图像是ori_im,目标图像是tar_im,tra_im[10,4]表示在行列分别是10和4时候图像的像素值。

此时,x=10/0.6=16.67,  y=4/0.6=6.67,而x1=16,x2=17,y1=6,y2=7, (x1,y1), (x1,y2), (x2,y1), (x2,y2)是在源图像中最接近tra_im[10,4]的4个点。

tra_im[10,4]=ori_im[x1,y1]*(17-16.67)*(7-6.67)+ori_im[x2,y1]*(16.67-16)*(6.67-6)+ori_im[x2,y1]*(17-16.67)*(6.67-6)+ori_im[x2,y2]*(16.67-16)*(6.67-6)

带入4个点在源图像中对应的像素值即可得到缩小后图像的像素值。

下面我用python实现了这个双线性插值,并和python中自带的skimage函数里封装的resize进行了效果对比,感觉效果差不多。(代码可能有点冗余)

我的代码:

# -*- coding: utf-8 -*-
# Author: lmh
# Time: 2018.10.22
from skimage import transform
import matplotlib.pyplot as plt
import matplotlib.image as mping
import numpy as np
def chazhi(x,y,im):
#x,y分别是缩放或者放大后对应源图像的浮点坐标位置,im是源图像,返回目标图像根据插值计算得到的像素值
    x1,y1=int(x),int(y) #x1,x2,x3,x4分别是插值坐标对应在源图像上下左右最近的点的坐标
    x2,y2=x1+1,y1+1
    pixel11,pixel21,pixel12,pixel22=im[x1-1,y1-1],im[x2-1,y1-1],im[x2-1,y1-1],im[x2-1,y2-1]
	#以下是根据双线性插值的公式求得的目标图像的该位置的像素值
    new_pixel=(x2-x)*(y2-y)*pixel11+(x-x1)*(y2-y)*pixel21+(x2-x)*(y-y1)*pixel12+(x-x1)*(y-y1)*pixel22
    return new_pixel


im=mping.imread('C:\\Users\\shou\\Desktop\\photo.png')
im11=im
scale=0.4 #缩小程度
row,col=int(im.shape[0]*scale),int(im.shape[1]*scale)
im_sml=np.zeros([row,col,3])
for k in range(3):
    im1 = im[:, :, k]
    for i in range(row):
        for j in range(col):
            value=chazhi(i/scale,j/scale,im1)#3通道图像逐像素计算缩小或者放大后的新像素值
            im_sml[i][j][k]=value

im_narrow=im_sml

scale=1.7  #扩大程度
row, col = int(im.shape[0] * scale), int(im.shape[1] * scale)
im_sml = np.zeros([row, col, 3])
for k in range(3):
    im1 = im[:, :, k]
    for i in range(row):
        for j in range(col):
            value = chazhi(i / scale, j / scale, im1)
            im_sml[i][j][k] = value

im_enlarge=im_sml
python_narrow=transform.resize(im, (175, 145)) #使用skimage自带函数resize图像,也可以直接写像上面一样写比例(0.4,1.7)
python_enlarge=transform.resize(im, (746, 617))#为了对比,两种方法特意放大和缩小一样大小

plt.figure()
plt.subplot(151)
plt.imshow(im11,plt.cm.gray)
plt.title('Original')
# plt.axis('off')

plt.subplot(152)
plt.imshow(im_narrow,plt.cm.gray)
plt.title('my_narrow')
# plt.axis('off')

plt.subplot(153)
plt.imshow(python_narrow,plt.cm.gray)
plt.title('skimage_narrow')
# plt.axis('off')

plt.subplot(154)
plt.imshow(im_enlarge,plt.cm.gray)
plt.title('my_enlarge')
# plt.axis('off')

plt.subplot(1,5,5)
plt.imshow(python_enlarge,plt.cm.gray)
plt.title('skimage_enlarge')
# plt.axis('off')
plt.savefig('image_comp.png')

skimage的源码:

 

 效果图对比(原图[439,363,3];缩小0.4倍后为[175,145,3];放大1.7倍后为[746,617,3])

在博客:https://www.cnblogs.com/sdxk/p/4056223.html中,博主讲根据双线性插值的定义,我们自己写的函数图像处理的结果会因为坐标系的原因,而和MATLAB,OpenCV结果的完全不同。最好的解决方法就是,两个图像的几何中心重合,并且目标图像的每个像素之间都是等间隔的,并且都和两边有一定的边距,这也是matlab和openCV的做法。并给出了以下解决方法,其中m,n 是源图像尺寸;a,b是目标函数尺寸。然而我没有这些修改,与python封装的方法也没发现太大区别。

int x=(i+0.5)*m/a-0.5

int y=(j+0.5)*n/b-0.5

代替

int x=i*m/a

int y=j*n/b

补:经过查看Skimge的resize源码,发现确实里面添加了上面所说的策略(下面的代码),那么为什么我写的双线性代码没有这样做,图像的结果确和这样做的几乎一样?更加疑惑了

# take into account that 0th pixel is at position (0.5, 0.5)
dst_corners[:, 0] = col_scale * (src_corners[:, 0] + 0.5) - 0.5
dst_corners[:, 1] = row_scale * (src_corners[:, 1] + 0.5) - 0.5

2. 协方差矩阵

 

 

 

3. 求矩阵的特征值,特征向量,以及主成分

 

 

 

2018-06-28 14:21:01 qq_40213457 阅读数 60

1.前言

看论文的时候又看到了协方差矩阵这个破东西,以前看图像处理的书籍的时候就特困扰,没想到现在还是搞不清楚,索性开始查协方差矩阵的资料,恶补之后决定马上记录下来。

2.拼出身—统计学的定义

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

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

那么问题来了,上面介绍的参量难道不能描述统计学所有特性吗?为啥还要搞出一个协方差出来

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

度量各个维度偏离其均值的程度,协方差可以这么来定义:

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

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



3.协方差矩阵的由来

好几十年前,鲁迅爷爷就说过,世界上本没有路,走的人多了也就有了路。协方差矩阵也是这样,好多协方差凑合到一起就形成了协方差矩阵。当然,数学的定义,不能如我这样随意。对于一个二维矩阵,每一个因子都可以视为两个不同随机变量的关系,这正好和协方差矩阵多少有点牵连,因此数学家们就把协方差矩阵引入到了二维矩阵中,衡量各个变量之间的紧密程度(就是关系度啦)。根据协方差的性质,我们可以类似的推出协方差矩阵的性质:


1.协方差矩阵一定是个对称的方阵

2.协方差矩阵对角线上的因子其实就是变量的方差:cov(X,X)=var(X)


这个定义还是很容易理解的,我们可以举一个简单的三变量的例子,假设数据集有{x,y,z}{x,y,z}三个维度,则协方差矩阵为:

再一次可以看出,协方差矩阵是一个对称的矩阵,而且对角线是各个变量上的方差。

3.MATLAB实战练习

上面涉及的内容都比较容易,协方差矩阵似乎也很简单,但实战起来就很容易让人迷茫了。必须要明确一点,### 协方差矩阵计算的是不同维度之间的协方差,而不是不同样本之间的
为了说明计算原理,不直接调用Matlab的cov函数:
首先,随机产生一个10*3维的整数矩阵作为样本集,10为样本的个数,3为样本的维数(3列)
MySample = fix(rand(10,3)*50)


根据公式,计算协方差需要计算均值,那是按行计算均值还是按列呢,我一开始就老是困扰这个问题。前面我们也特别强调了,协方差矩阵是计算不同维度间的协方差,要时刻牢记这一点。样本矩阵的每行是一个样本,每列为一个维度,所以我们要按列计算均值。为了描述方便,我们先将三个维度的数据分别赋值:
  1. dim1 = MySample(:,1);1列
  2. dim2 = MySample(:,2);2列
  3. dim3 = MySample(:,3);3列
计算dim1与dim2,dim1与dim3,dim2与dim3的协方差:
  1. sum( (dim1-mean(dim1)) .* (dim2-mean(dim2)) ) / ( size(MySample,1)-1 ) % 得到 74.5333
  2. sum( (dim1-mean(dim1)) .* (dim3-mean(dim3)) ) / ( size(MySample,1)-1 ) % 得到 -10.0889
  3. sum( (dim2-mean(dim2)) .* (dim3-mean(dim3)) ) / ( size(MySample,1)-1 ) % 得到 -106.4000
搞清楚了这个后面就容易多了,协方差矩阵的对角线就是各个维度上的方差,下面我们依次计算:
  1. std(dim1)^2 % 得到 108.3222
  2. std(dim2)^2 % 得到 260.6222
  3. std(dim3)^2 % 得到 94.1778
这样,我们就得到了计算协方差矩阵所需要的所有数据:

C11=108.3222C12=74.5333C13=-10.0889
C21=74.5333C22=260.6222C23=-106.4000
C31=-10.0889C32=-106.4000C33=94.1778


调用Matlab自带的cov函数进行验证:


4.心得感悟

理解协方差矩阵的关键就在于牢记它计算的是不同维度之间的协方差,而不是不同样本之间,拿到一个样本矩阵,我们最先要明确的就是一行是一个样本还是一个维度,心中明确这个整个计算过程就会顺流而下,这么一来就不会迷茫了。
其实还有一个更简单的容易记还不容易出错的方法:协方差矩阵一定是一个对称的方阵,一定是一个对称的方阵一定是一个对称的方阵!!!记住就好啦~

浅谈协方差矩阵

阅读数 273

没有更多推荐了,返回首页