REOF 的matlab代码

简单幸福

最近在做REOF分析，虽然网上有很多Fortran代码，但对于已经习惯了matlab的童鞋们，更倾向于用matlab来搞定。下面程序为REOF的matlab代码，供大家参考使用（新建一个名为varimax的m文件，复制下面代码，使用时然后调用varimax函数即可）。送给那些金钱跟我一样少的同学！更多内容请参考网址：http://pmc.ucsc.edu/~dmk/notes/EOFs/


function [x, r] = varimax( x, normalize, tol, it_max )
% VARIMAX - Rotate EOF's according to varimax algorithm
%
% This is actually a generic varimax routine and knows nothing special about
% EOFs.  It expects a matrix of "loadings".  Typically (in state space
% rotation), these loadings are the expansion coefficients (aka Principal
% Component Time Series) for the truncated basis of eigenvectors (EOFs), but
% they could also be the EOFs*diag(L)^(1/2) (in the case of rotation in
% sample space).
%
% Usage: [new_loads, rotmax] = varimax( loadings, normalize, tolerance, it_max )
%
% where all but the loadings are optional.  rotmax is the rotation matrix used.
%
% normalize determines whether or not to normalize the rows or columns of
% the loadings before performing the rotation.  If normalize is true, then
% the rows are normalized by there individual lengths.  Otherwise, no
% normalization is performed (default).  After rotation, the matrix is
% renormalized. Normalizing over the rows corresponds to the Kaiser
% normalization often used in factor analysis.
%
% tolerance defaults to 1e-10 if not given.  it_max specifies the maximum
% number of iterations to do - defaults to 1000.
%
% After the varimax rotation is performed, the new EOFs (in the case that
% the EC's were rotated - state space) can be found by new_eofs =
% eofs*rotmax.
%
% This function is derived from the R function varimax in the mva
% library.   
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  $Id: varimax.m,v 1.4 2002/10/10 00:28:30 dmk Exp $
%
% Copyright (C) 2002 David M. Kaplan
% Licence: GPL
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if nargin < 2
  normalize = 0;
end
if nargin < 3
  tol = 1e-10;
end
if nargin < 4
  it_max = 1000;
end
[p, nc] = size(x);
if nc < 2, return; end
if normalize
  rl = repmat( sqrt(diag( x*x' )), [1,nc] ); % By rows.
  % rl = repmat( sqrt(diag( x'*x ))', [p,1] ); % By columns.
  x = x ./ rl;
end
TT = eye( nc );
d = 0;
for i = 1 : it_max
  z = x * TT;
  B = x' * ( z.^3 - z * diag(squeeze( ones(1,p) * (z.^2) )) / p );
  
  [U,S,V] = svd(B);
  
  TT = U * V';
  d2 = d;
  d = sum(diag(S));
  
  % End if exceeded tolerance.
  if d < d2 * (1 + tol), break; end
  
end
% Final matrix.
x = x * TT;
% Renormalize.
if normalize
  x = x .* rl;
end
if nargout > 1
  r = TT;
end

zosuper

在[p, nc] = size(x);这一句之前给x一个值，例如eof分解出来空间场的前4个场，然后计算旋转后的结果。
与直接用[xr,T] = rotatefactors(x,'method','varimax','Normalize','off');计算出来的一样呀！
也就是说如果有了x,直接调用这个函数就可以实现帖子给出的功能了。xr就是旋转后的结果。如果y是eof分解得到的空间场，用前4个特征向量的话，[yr,T] = rotatefactors(y(:,1:4),'method','varimax','Normalize','off');即可

xiaowan

function [x, r] = varimax( x, normalize, tol, it_max )这里面的normallize怎么填呢~~是填矩阵还是是否呢

龙行天下

哈哈，虽然你这代码是直接从别人帖子里粘贴过来的，但是我还是想说一句，你这种直接盗用别人成果的行为真是太好了，少扣了我5分

简单幸福

破小财消大灾！哈哈，谢谢！

简单幸福

真的是哎，刚说了，就好运来啦！

维欻?

这个程序的用法，楼主介绍一下吧

又是那隻貓

谢谢分享！！

kongfeng0824

感谢分享。楼主有个例子介绍一下出来的图的结果吗？

阿吉木柯

先保存下来，谢谢楼主。
我刚看了一下网站，也得把函数下下来的，不过网站上就有很方便

1079652110

还是楼主好

龙行天下

楼主能不能举个例子，详细介绍一下用法？