主成分分析及一些函数.doc

主成分分析
主成分概念首先由 Karl Parson在1901年引进，当时只对非随机变量来讨论的。1933年Hotelling将这个概念推广到随机变量。
在多数实际问题中，不同指标之间是有一定相关性。由于指标较多及指标间有一定的相关性，势必增加分析问题的复杂性。
主成分分析就是设法将原来指标重新组合成一组新的互相无关的几个综合指标来代替原来指标。同时根据实际需要从中可取几个较少的综合指标尽可能多地反映原来的指标的信息。

主成分分析（Principal Component Analysis，PCA）， 将多个变量通过线性变换以选出较少个数重要变量的一种多元统计分析方法。又称主分量分析。在实际课题中，为了全面分析问题，往往提出很多与此有关的变量（或因素），因为每个变量都在不同程度上反映这个课题的某些信息。主成分分析首先是由K.皮尔森对非随机变量引入的，尔后H.霍特林将此方法推广到随机向量的情形。信息的大小通常用离差平方和或方差来衡量。

主成分分析是设法将原来众多具有一定相关性（比如P个指标），重新组合成一组新的互相无关的综合指标来代替原来的指标。通常数学上的处理就是将原来P个指标作线性组合，作为新的综合指标。最经典的做法就是用F1（选取的第一个线性组合，即第一个综合指标）的方差来表达，即Var(F1)越大，表示F1包含的信息越多。因此在所有的线性组合中选取的F1应该是方差最大的，故称F1为第一主成分。如果第一主成分不足以代表原来P个指标的信息，再考虑选取F2即选第二个线性组合，为了有效地反映原来信息，F1已有的信息就不需要再出现再F2中，用数学语言表达就是要求Cov(F1, F2)=0，则称F2为第二主成分，依此类推可以构造出第三、第四，……，第P个主成分。

程序：
function [lambda,T,fai]=MSA2(A)
%求标准化后的协差矩阵,再求特征根和特征向量
%标准化处理
[p,n]=size(A);
for j=1:n
mju(j)=mean(A(:,j));
sigma(j)=sqrt(cov(A(:,j)));
end
for i=1:p
for j=1:n
Y(i,j)=(A(i,j)-mju(j))/sigma(j);
end
end
sigmaY=cov(Y);
%求X标准化的协差矩阵的特征根和特征向量
[T,lambda]=eig(sigmaY);
disp('特征根(由小到大):');
disp(lambda);
disp('特征向量:');
disp(T);
%方差贡献率;
Xsum=sum(sum(lambda,2),1);
for i=1:n
fai(i)=lambda(i,i)/Xsum;
end
disp('方差贡献率:');
disp(fai);
u=T(:,n);
B=[];
h=length(A(:,1));
for k=1:n
m1=mean(A(:,k));
t=(A(:,k)-m1).^2;
m2=sqrt(sum(t))/(h-1);
B=[B,(A(:,k)-m1)./m2];
End

V = var(X)
如果X是一个向量，返回向量X的方差。
如果X是一个矩阵，var(X)返回一个包含矩阵X每一列方差的行向量。
如果X是一个N维数组，var沿着第一个X的非单一维进行操作。
只要X是独立同分布的，结果V是X分布的总体方差的无偏估计。
当N>1时，var由N-1来标准化，其中N是样本大小。
只要样本是独立同分布的，它就是X分布的总体方差的无偏估计。对N=1来说，v由N来标准化。

V = var(X,1)
由N来标准化，并且生成了样本关于其均值的二阶矩，var(X,0)等价于var(X)。

V = var(X,w)
计算向量X的方差利用权重向量w，向量w中元素的数目必须和X中的列的数目相同，向量w中的元素必须全是正数。var归一化w是的总和为1。

V = var(X,w,dim)
沿着指定维数dim求X的方差，默认用N-1标准化这时w为0，w为1时用N标准化。

方差是其标准差(STD)的平方.

应用举例 编辑本段回目录x=[4,6,4,3,5,7]

var(x)
ans = 2.1667

R = corrcoef(X)
返回相关系数矩阵R，对行是观测值、列是变量的矩阵X计算相关系数得到相关系数矩阵R。
阵R=corrcoef(X)与协方差矩阵C=cov(X)有以下关系:

R = corrcoef(x,y)
如果x和y都是列向量，那么其功能和corrcoef([x y])是相同的。
如果x和y不是列向量，corrcoef函数转化他们为列向量，例如,在这种情况下R=corrcoef(x,y)是等价于R=corrcoef([x(:) y(:)])。 book.iLoveM
[R,P]=corrcoef(...)
返回R,P,含有p值的矩阵是为了测试不相关性，每个p值是相关性的概率值，真正的相关性，p值应该为0。
如果p(i,j)很小的情况下，一般设置小于0.05，我们也认为相关系数R(i,j)也是有意义的。
[R,P,RLO,RUP]=corrcoef(...)
返回矩阵RLO 和RUP，它们是和R有相同的大小，对每个系数包含一个在95%置信区间的上下界。
[...]=corrcoef(...,'param1',val1,'param2',val2,...) www.iLoveM
指定其他参数和它们的值，有效的参数如下：
'alpha' 是一个0到1之间的数，指定了100*(1 - alpha)%的置信水平，默认是0.05 相对于95%的置信区间。
'rows' 'all' (默认)作用于所有行，'complete'作用于没有非数(NaN)的行，'pairwise' 是计算R(i,j)用列i或者j中不是非数的值。

应用举例 编辑本段回目录生成随机数据其中第4列和其他列相关。
x = randn(30,4); % Uncorrelated data
x(:,4) = sum(x,2); % Introduce correlation.
[r,p] = corrcoef(x) % Compute sample correlation and p-values.
[i,j] = find(p<0.05); % Find significant correlations.
[i,j] % Display their (row,col) indices.
r =
1.0000 -0.3566 0.1929 0.3457
-0.3566 1.0000 -0.1429 0.4461
0.1929 -0.1429 1.0000 0.5183
0.3457 0.4461 0.5183 1.0000
p =
1.0000 0.0531 0.3072 0.0613
0.0531 1.0000 0.4511 0.0135
0.3072 0.4511 1.0000 0.0033
0.0613 0.0135 0.0033 1.0000
ans =
4 2
4 3
2 4
3 4

cov(x)
如果x是一个向量，cov(x)返回x中元素的方差。
如果x是一个矩阵，并且是各行为观测值、各列为变量的矩阵，cov(x)计算x的协方差。diag(cov(X))是包含每列的方差的向量，sqrt(diag(cov(X)))是包含标准差的向量。当X和Y是大小相同的矩阵时，cov(X,Y)等价于cov([X(:) Y(:)])。 Matlab中文论坛

cov(x) 或 cov(x,y)
用N-1规范化如果N>1，其中N是观测值的个数。如果观测值属于正态分布，那么cov(x)可以达到协方差最好的无偏估计。如果N=1用N去规范化。

cov(x,1)或 cov(x,y,1) www.iLoveM
用N来规范化，cov(X,Y,0)和cov(X,Y)是功能相同的，cov(X,0)和cov(X)是功能相同的。

备注 编辑本段回目录cov在计算结果前对每一列删除平均值。
协方差函数定义为：

其中E是数学期望并且有
应用举例 编辑本段回目录考虑A = [-1 1 2 ; -2 3 1 ; 4 0 3]，对A每一列求方差：
v = diag(cov(A))'
v =
10.3333 2.3333 1.0000
比较v和协方差C：
C=cov(A)
C =
10.3333 -4.1667 3.0000
-4.1667 2.3333 -1.5000
3.0000 -1.5000 1.0000

n是元素的个数，上面两个等式的不同仅仅在于被除数，一个是n，一个是n-1。

s = std(X)
如果X是一个向量，用上面公式(1)返回标准差。
结果s是数据X分布的无偏估计方差的平方根，只要数据X是有独立同分布样本组成。
如果X是一个矩阵，std(X)返回一个包含每一列所有元素的标准差的行向量。
如果A是一个多元数组，std(X)将沿着数组X中第一个非单一维的值来计算标准差。

s = std(X,flag)
若flag = 0，这个命令和s = std(X)功能相同。
若flag = 1，std(X,1)用上面公式2来计算标准差。

s = std(X,flag,dim)
用指定的标量dim沿着数据X的维数来计算标准差，设置flag为0用n-1去标准化X，设置flag为1用n去标准化。

应用举例 编辑本段回目录对于矩阵X
X =
1 5 9
7 15 22

s = std(X,0,1)
s =
4.2426 7.0711 9.1924

s = std(X,0,2)
s =
4.000
7.5056

d=eig(A)

返回矩阵A特征值的一个向量d。
d=eig(A,B)

如果A和B是方阵的，返回包含广义特征的向量。
注意：如果S是稀疏对称的，用d = eig(S) 可以返回S的特征值。如果S是稀疏但不对称，若要想得到S的特征向量，用函数eigs代替eig。
[V,D]=eig(A)

计算矩阵A的特征值D和特征矩阵V，满足A*V = V*D。矩阵D是矩阵A的规范形式：主对角线上的元素是矩阵A的特征值的对角矩阵。
[V,D] = eig(A,'nobalance')
当矩阵A中有与截断误差数量级相差不远的值时，改命令可能更加精确，'nobalance'起误差调节作用。

[V,D] = eig(A,B)
计算广义特征向量矩阵V和广义特征值矩阵D，满足A*V = B*V*D 。
[V,D] = eig(A,B,flag)
由flag指定算法来计算特征值D和特征向量V，flag的可能值为：
'chol'：表示对B进行cholesky分解算法，这里A是对称的Hermitian矩阵，B为正定矩阵。

'qz'：表示使用qz分解算法，这里A、B是非对称或非Hermitian矩阵。