产品
博主最近在研究奇异值分析，参考学姐的程序自己写了一个函数。
function [lambda,L,M,u,v,SCF,CSCF,cor,cor_xu,cor_yv,cor_xv,cor_yu] = svdan(X,Y) % 奇异值分析程序 % 参考了赵梦楠学姐的程序 %X,Y为待分析的2个要素场，各输出变量的含义下面会有介绍 [p,n] = size(X);% 第一组要素X:p个地点,n次观测 [q,nn] = size(Y);% 第二组要素Y:q个地点,nn次观测 if n~=nn %SVD分析的前提是X和Y观测次数必须相等 disp(['error']); end for i = 1:p X0(i,:) = (X(i,:)-mean(X(i,:)))/std(X(i,:));%求距平后标准化 end for i = 1:q Y0(i,:) = (Y(i,:)-mean(Y(i,:)))/std(Y(i,:)); end % 求协方差阵C C= X0*Y0'/n; %SVD分解--------------------------------------------- [U,S,V]=svd(C);%这里使用了matlab自带的奇异值分解函数,U为p阶方阵，V为q阶方阵 lambda = diag(S);%求得奇异值 index=find(abs(lambda)==0 | abs(lambda)<10e-8);%找到为0或者是由于浮点计算产生的非常小但不是0的奇异值 lambda(index)=[];%去掉为0的奇异值 r1=size(lambda);r=r1(1,1);%r为非0奇异值的个数 L=U(:,1:r);%L的每列都代表了一个空间模态,又称L为左奇异向量 M=V(:,1:r);%M的每列也代表了一个空间模态，又称M为右奇异向量 %这里特别说明一下：在run程序时发现本程序算出来的L和M和赵学姐的程序算出来的结果出现了不一致，有时某些列的正负号正好相反， %但这应该影响不大，因为L和M的每列都是标准正交基 u = L'*X0;%u对应Y的主成分（时间系数） v = M'*Y0;%v对应Y的主成分（时间系数） %--------------------------------------------- % 求第k个SVD模态协方差平方贡献率SCF SCF(1:r,1) = lambda.^2./sum(lambda.^2); % 求前s个SVD模态累积协方差平方贡献率CSCF CSCF(1:r,1) = SCF(1); for s = 2:r CSCF(s) = CSCF(s-1)+SCF(s); end %------------------------------------------ % 求第k个SVD模态时间系数之间的相关系数 cor = lambda./sqrt(var(u').*var(v'))'; for i = 1:r % 求X0与u之间的同性相关系数cor_xu（p行r列，p行代表p个地点，r列代表r个模态） for j = 1:p tmp = corrcoef(X0(j,:),u(i,:));%matlab自带的求相关系数的命令 cor_xu(j,i) = tmp(1,2); end end for i = 1:r % 求Y0与v之间的同性相关系数cor_yv for j = 1:q tmp = corrcoef(Y0(j,:),v(i,:)); cor_yv(j,i) = tmp(1,2); end end for i = 1:r % 求X0与v之间的异性相关系数cor_xv for j = 1:p tmp = corrcoef(X0(j,:),v(i,:)); cor_xv(j,i) = tmp(1,2); end end for i = 1:r % 求Y0与u之间的异性相关系数cor_yu for j = 1:q tmp = corrcoef(Y0(j,:),u(i,:)); cor_yu(j,i) = tmp(1,2); end end %PS.如果你画cor_yu，cor_yv,cor_xu,cor_xv这四个场的话，你会发现cor_xu,cor_xv差不多，cor_yu，cor_yv差不多，如果你要研究要素Y随X的变化，个人认为你应该分析cor_yu，反之，如果你研究要素X随Y的变化，应该分析cor_xv.
155-2731-8020
