Matlab 求矩阵均值和协方差矩阵 ∈ Matlab 使用笔记
今天想用 Matlab 来算一个矩阵数据的均值和方差给我整迷了，好气！！所以做点笔记，以免下次在跳进坑里。
主博客：https://blog.csdn.net/Gou_Hailong/article/details/106092705
1.预备知识
关于算方差等，matlab提供了好几个函数：mean、var、std、cov，首先我们先看看这几个函数咋用的。
以下用到的矩阵都为：
X = [1 2 3; 3 3 6; 4 6 8; 4 7 7; 8 5 9]
=============================
X =
1 2 3
3 3 6
4 6 8
4 7 7
8 5 9
1.mean
/*这个函数还好理解，算平均值
如果参数是向量就直接是平均值
如果参数是矩阵，就会有俩种模式
1 返回一个行向量，计算各列的平均值（默认）
2 返回一个列向量，计算各行的平均值*/
>> mean(X(:,1))
ans =
4
>> mean(X(1,:))
ans =
2
>> mean(X)
ans =
4.0000 4.6000 6.6000
>> mean(X,2)
ans =
2.0000
4.0000
6.0000
6.0000
7.3333
如果一个矩阵中有nan的话，如何来求这个矩阵的均值呢？
AAt = reshape(AA, size(AA,1)*size(AA,2),1);
mean(AAt,'omitnan');
2.var
这个是用来算方差的，但是分母为（n-1），好像是为了搞成无偏估计，概率论的知识，忘了（😓）
var(X,W,DIM)
/*这个函数求方差，公式参上，默认：var(X,0,1)
参数X为向量时：w=0:n-1 w=1:n
参数X为矩阵时：w=0:n-1 w=1:n
dim=1 返回一个行向量，计算各列的方差（默认）
dim=2 返回一个列向量，计算各行的方差
*/
>> var(X(1,:))
ans =
1
>> var(X(:,1))
ans =
6.5000
>> var(X,1,1)
ans =
5.2000 3.4400 4.2400
>> var(X,1,2)
ans =
0.6667
2.0000
2.6667
2.0000
2.8889
>> var(X,0,1)
ans =
6.5000 4.3000 5.3000
>> var(X,0,2)
ans =
1.0000
3.0000
4.0000
3.0000
4.3333
3.std
这个函数用来求标准差，就是var开方，调用形式和var类似都是var(X,W,DIM) 这里就不在赘述。
4.cov
我们先来回顾一下协方差矩阵是咋算的：
设X为一个向量，其均值表示为E(X)，那么它的方差为：D(X)=E((E(X)-Xi)^2)
这时，又有另一个向量Y(与X长度相同) ,其均值表示为E(Y)
那么它俩的协方差为：cov(X,Y)=E((E(X)-Xi)*(E(Y)-Yi)) = E(XY)－E(X)*E(Y)
协方差与方差有如下关系：
D(X+Y)=D(X)+D(Y)+2Cov(X，Y)
D(X-Y)=D(X)+D(Y)-2Cov(X，Y)
ok，接下来我们一起用matlab 来算算。
cov(X,Y,w)
/*用来计算向量X和Y 的协方差矩阵
没有Y时，用来算X的方差。
w=0:n-1 w=1:n 0 默认
当X是矩阵时：
算的是将每一列当做一个X(或Y),
X若为mxn 的矩阵，得到的是nxn的协方差矩阵*/
>> cov(X,1)
ans =
5.2000 2.2000 4.2000
2.2000 3.4400 2.8400
4.2000 2.8400 4.2400
>> cov(X,0)
ans =
6.5000 2.7500 5.2500
2.7500 4.3000 3.5500
5.2500 3.5500 5.3000
>> cov(X(:,1),X(:,2))
ans =
6.5000 2.7500
2.7500 4.3000
ok，至此已经介绍完了，下面的内容是我自己的思想斗争，大家可以忽略，忽略…
5.一个奇怪的东西
我现在突然意识到是我的指导思想错了，我写这篇博客的目的是为了得到一组遥感影像数据（某一类）的协方差矩阵，可现在发现一个cov 完全可以解决，之前想的是：
其中n是n个波段，k是k个像素点，也就是说，每个像素点有n个灰度；也即有一组数据，有n个维度，这组数据共有k个。之后想求其协方差矩阵。
均值肯定是这k个样本点的均值，为一个n维的向量
从一维向量（行向量1xk）求方差 合理外延就会想到上述公式，但是这个并不是我所需要的，这个是一个kxk的矩阵，我想要的是一个nxn的矩阵，这两个东西到底是何含义，想到爆也想不出来（看来需要时间的积淀了！）
还是看个实例比较容易理解一点。
下面的实例是用来解决上面那个公式引申出来的。
2.实例
1.我们先建个矩阵
X = [1 2 3; 3 3 6; 4 6 8; 4 7 7; 8 5 9]
=============================
X =
1 2 3
3 3 6
4 6 8
4 7 7
8 5 9
如上，这是一个5行乘3列的矩阵，我们假设这个矩阵每一列代表遥感影像一个位置处的5个波段的灰度值，也就是说，有个遥感影像，他有三个点，5波段。我们想算他的均值和协方差矩阵，期望得到的均值是一个5行1列的列向量，希望得到的协方差矩阵是一个3行3列 （实际上是5x5）的矩阵。
我们还想得到一种协方差矩阵，它是5行5列的矩阵。
2.为了验证，matlab 计算的正确性，我们先手动算一遍：
均值向量为：
mean=
2.0000
4.0000
6.0000
6.0000
7.3333
协方差矩阵1为：
cov=
10.4444 -2.5556 -7.8889
-2.5556 7.4444 -4.8889
-7.8889 -4.8889 12.7778
协方差矩阵2为：
cov=
10.4444 -2.5556 -7.8889
-2.5556 7.4444 -4.8889
-7.8889 -4.8889 12.7778
3.我们来看看用 matlab 自带的函数mean、cov算出来的都是啥：
>> mean(X,1)
ans =
4.0000 4.6000 6.6000
>> mean(X,2)%很遗憾，只有这一个达到了我们的期望
ans =
2.0000
4.0000
6.0000
6.0000
7.3333
>> cov(X)
ans =
6.5000 2.7500 5.2500
2.7500 4.3000 3.5500
5.2500 3.5500 5.3000
那我们将X转置一下在计算会不会好点呢？
>> X=X'
X =
1 3 4 4 8
2 3 6 7 5
3 6 8 7 9
>> mean(X,1)
ans =
2.0000 4.0000 6.0000 6.0000 7.3333
>> mean(X,2)
ans =
4.0000
4.6000
6.6000
>> cov(X)
ans =
1.0000 1.5000 2.0000 1.5000 0.5000
1.5000 3.0000 3.0000 1.5000 2.5000
2.0000 3.0000 4.0000 3.0000 1.0000
1.5000 1.5000 3.0000 3.0000 -1.0000
0.5000 2.5000 1.0000 -1.0000 4.3333
结果还是不符合我们的期望。这是为什么呢？原来这个Matlab计算协方差矩阵的时候，是用每一列作为一个元素（这与我们的期望相符）对角线的方差是计算这一列中的方差（就是其他的元素不考虑，算出这一列中的几个数的均值，之后在算这几个数的方差），非对角线咋算呢？非对角线列是相应列的协方差。
实在想不清楚，可以看看这篇文章：
https://blog.csdn.net/clam_clam/article/details/7335588
3.解决方案
OK，到目前为止，我们已经知道指望不上Matlab 自带的函数了（均值可以求，方差不行），那么我们就自食其力，自己写了个函数，用来计算上面想得到的方差。
%mycov.m by@GHL
%此函数用来计算方差，
%m是一个n行，k列的矩阵
%每一个元素是m中的列向量，有n维
%m中共有k个元素
%最后得到的协方差矩阵是一个k行k列的矩阵
function sigma = mycov(m)%计算
X=mean(m,2);%均值
n=size(m,1);%n维
k=size(m,2);%k个元素
sigma=zeros(k);
for i = 1:k
for j = 1:k
sigma(i,j)=(m(:,i)-X)'*(m(:,j)-X);
end
end
return;
end
验证：
>> X=X'
X =
1 2 3
3 3 6
4 6 8
4 7 7
8 5 9
>> mycov(X)
ans =
10.4444 -2.5556 -7.8889
-2.5556 7.4444 -4.8889
-7.8889 -4.8889 12.7778
参考/引用 文章
[1] clam_clam-CSDN博主：https://blog.csdn.net/clam_clam/article/details/7335588
2020/6/16
【小感慨】
这么一个小玩意花了我一下午来捋清，好气！！差点爆粗口，真是用matlab 的基础还很不好，不管怎么说，这次终于搞清楚了（只差哪两个协方差矩阵没搞明白），希望下次不会在迷了😞😞😞