请选择 进入手机版 | 继续访问电脑版
爱气象,爱气象家园! 

气象家园

 找回密码
 立即注册

QQ登录

只需一步,快速开始

新浪微博登陆

只需一步, 快速开始

搜索
查看: 133805|回复: 126

[程序设计] REOF 的matlab代码

  [复制链接]

新浪微博达人勋

发表于 2013-7-29 11:31:45 | 显示全部楼层 |阅读模式

登录后查看更多精彩内容~

您需要 登录 才可以下载或查看,没有帐号?立即注册 新浪微博登陆

x
最近在做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



评分

参与人数 1金钱 +1 收起 理由
li_panda520 + 1 赞一个!

查看全部评分

密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2014-1-4 20:17:25 | 显示全部楼层
本帖最后由 zosuper 于 2014-1-4 20:18 编辑

在[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');即可
密码修改失败请联系微信:mofangbao
回复 支持 2 反对 0

使用道具 举报

新浪微博达人勋

发表于 2014-4-17 19:48:17 | 显示全部楼层
function [x, r] = varimax( x, normalize, tol, it_max )这里面的normallize怎么填呢~~是填矩阵还是是否呢
密码修改失败请联系微信:mofangbao
回复 支持 1 反对 0

使用道具 举报

新浪微博达人勋

发表于 2013-8-31 12:51:23 | 显示全部楼层
哈哈,虽然你这代码是直接从别人帖子里粘贴过来的,但是我还是想说一句,你这种直接盗用别人成果的行为真是太好了,少扣了我5分
密码修改失败请联系微信:mofangbao
回复 支持 1 反对 0

使用道具 举报

新浪微博达人勋

 楼主| 发表于 2013-7-29 11:33:42 | 显示全部楼层
破小财消大灾!哈哈,谢谢!
密码修改失败请联系微信:mofangbao

新浪微博达人勋

 楼主| 发表于 2013-7-29 11:34:34 | 显示全部楼层
真的是哎,刚说了,就好运来啦!
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2013-7-29 11:42:08 | 显示全部楼层
这个程序的用法,楼主介绍一下吧
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2013-7-29 11:43:37 | 显示全部楼层
谢谢分享!!
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2013-7-29 11:47:10 | 显示全部楼层
感谢分享。楼主有个例子介绍一下出来的图的结果吗?
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2013-8-16 08:10:30 | 显示全部楼层
先保存下来,谢谢楼主。
我刚看了一下网站,也得把函数下下来的,不过网站上就有很方便
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2013-8-16 09:25:14 | 显示全部楼层
还是楼主好
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2013-9-2 14:52:13 | 显示全部楼层
楼主能不能举个例子,详细介绍一下用法?
密码修改失败请联系微信:mofangbao
您需要登录后才可以回帖 登录 | 立即注册 新浪微博登陆

本版积分规则

Copyright ©2011-2014 bbs.06climate.com All Rights Reserved.  Powered by Discuz! (京ICP-10201084)

本站信息均由会员发表,不代表气象家园立场,禁止在本站发表与国家法律相抵触言论

快速回复 返回顶部 返回列表