function [L,Y,Uk] = eof(x,miss,p)
%EOF Empirical Othorgonal Function
% [L,Y,Uk] = eof(x,miss,p)
% 输入参数：x：m*n的矩阵，m为空间点，n为时间，按标准格式排放。
% x可以是原资料，距平资料和标准化资料。但注意相应得到的模态和时间系数不同！
% miss：输入缺测值，用于程序识别资料的缺测值，以便对有效资料做eof分析
% p：累积方差贡献率，通常取为0.8~0.9，最好不要太大，可取得最主要的模态。（如：p取0.8）
% 输出参数：L：通过显著性检验的空间模态。按列存放，第一列为方差贡献最大的模态，依此类推。
% L中还原了缺测值miss。
% Y：时间系数。按行存放，第一行为对应第一空间模态的时间系数，依此类推。
% Uk：通过显著性检验的方差贡献，与L的顺序相对应。
% Version 1.0
% By Aires, Jun 6, 2013
% 找到矩阵x中有意义的数据和非数（NaN）的位置
[m1 n]=size(x);
na=find(x(:,1)==miss);
a=find(x(:,1)~=miss);
missing_val=x(na,:);
% 用有意义的x做eof分解
x=x(a,:);
[m n]=size(x);
% 准备数据
L1=zeros(m,n);
% 做时空转换，求x的特征值和空间模态
c=x'*x; % 注意：没有乘1/n，得到的特征值不同，但得到的空间模态和时间系数是相同的。
[LQ,lambda]=eig(c);
lambda=fliplr(flipud(lambda));
lambda=diag(lambda); % 得到的是列向量
LQ=fliplr(LQ);
LR=x*LQ;
for i=1:length(lambda)
L1(:,i)=LR(:,i)/sqrt(lambda(i));
end
% 根据给出的累积方差贡献率p，取前几个主要特征值，方差贡献和空间模态
total=sum(lambda);
Uk=lambda/total; % 方差贡献率
Gp=cumsum(Uk); % 累积方差贡献率
b=find(Gp>p);
lambda=lambda(1:b(1));
Uk=Uk(1:b(1));
L1=L1(:,1:b(1));
% 对选出的前几个主要特征值做显著性检验
% 注：该显著性检验（95%置信度水平）由North（1982年）提出
% 式中n应为有效自由度，由于有效自由度公式各异，这里粗略取为样本容量n
ej=lambda(1:length(lambda)-1)*sqrt(2/n);
for i=1:length(lambda)-1
if (lambda(i)-lambda(i+1))<ej(i)
lambda=lambda(1:i);
break
end
end
Uk=Uk(1:length(lambda));
L1=L1(:,1:length(lambda));
% 求出时间系数，并恢复L（含缺测值miss）
Y=L1'*x;
L=zeros(m1,length(lambda));
missing_val=missing_val(:,1:length(lambda));
L(a,:)=L1;
L(na,:)=missing_val;
end
- 1
- 2
- 3
前往页