eofrot.m
function [coef,Arot,varrotated]=eofrot(z,index)
% Algorithm for eof rotation with the varimax method
% Difiere de la eofrot de Antonio, pues le agregue en el "sort" los RPC
% (coef)
%
%Inputs:
% z Data Matrix
% index EOF to be rotated
%
%Outputs:
% coef Cofficient for the rotated EOFs
% Arot Rotated EOF in ascending order
% varrotated Variance explained by rotated EOFs
% Time and space points
[npoints,ntime]=size(z);
% Compute unrotated EOF for variance calculation
[uneof,ss,vneof]=svd(z,0); diag(ss);
totvar = sum(diag(ss.^2));
lding = uneof(:,index);
sl = diag(ss(index,index).^2);
% variance explained by the unrotated modes
varexpl = sum(sl)/totvar;
[n,nf]=size(lding);
b=lding;
hl=lding*(diag(sl));
hjsq=diag(hl*hl');
hj=sqrt(hjsq); % Normalize by the communalities
bh=lding./(hj*ones(1,nf));
Vtemp=n*sum(sum(bh.^4))-sum(sum(bh.^2).^2); % VARIMAX functional to be minimized
V0=Vtemp;
for it=1:10; % Number of iterations
for i=1:nf-1; % Program cycles through 2 factors
jl=i+1; % at a time
for j=jl:nf;
xj=lding(:,i)./hj; % notation here closely
yj=lding(:,j)./hj; % follows harman
uj=xj.*xj-yj.*yj;
vj=2*xj.*yj;
A=sum(uj);
B=sum(vj);
C=uj'*uj-vj'*vj;
D=2*uj'*vj;
num=D-2*A*B/n;
den=C-(A^2-B^2)/n;
tan4p=num/den;
phi=atan2(num,den)/4;
angle=phi*180/pi;
[i j it angle];
if abs(phi)> eps;
Xj=cos(phi)*xj+sin(phi)*yj;
Yj=-sin(phi)*xj+cos(phi)*yj;
bj1=Xj.*hj;
bj2=Yj.*hj;
b(:,i)=bj1;
b(:,j)=bj2;
lding(:,i)=b(:,i);
lding(:,j)=b(:,j);
end
end
end;
lding=b;
bh=lding./(hj*ones(1,nf));
Vtemp=n*sum(sum(bh.^4))-sum(sum(bh.^2).^2);
V=Vtemp;
fprintf(' Converging V V0 \n');
fprintf(' %f\n',[V V0])
if abs(V-V0)<.0001;break;else V0=V;end;
end;
% Reflect vectors with negative sums
for i = 1:nf
if sum(lding(:,i)) < 0
lding(:,i) = -lding(:,i);
end
end
Arot=lding ; % rotated eof
coef=z'*Arot(:,1:nf); % time series for rotated eof
whos Arot coef
for i=1:nf
varex(i) = sum(var(coef(:,i)*Arot(:,i)')*(ntime-1));
end
varexpl
varexplrot = sum(varex)/totvar
zvar=sum(var(z')*(ntime-1))
totvar
% Sort in decreasing order of variance
whos varex
varex
[varex,I]=sort(varex);
Arot=Arot(:,I);
Arot = fliplr(Arot);
coef=coef(:,I);
coef=fliplr(coef);
varex = flipud(varex');
varunrotated = sl/totvar
varrotated = varex/totvar
return
PS:我从别的地方copy的{B.Aires, 20-24/02/06 - Centro de Investigaciones del Mar y la Atmosfera, Department of Atmospheric and Oceanic Sciences (UBA)
E.Scoccimarro, A.F.Carril }