varimax.m

本文件夹中包含了克立格插值

function [x, r] = varimax( x, normalise, tol, it_max )

% HEADER_START
%
% NAME:
%
%   varimax
%
% DESCRIPTION:
%
%   VARIMAX - Rotate EOF's according to varimax algorithm
%
% REQUIRED_ARGS:
%
%   loadings:   the matrix of EOFs (or loadings) or PCs to be rotated.
%
% OPTIONAL_ARGS:
%
%   normalise:  Logical: T or F for normalising (or not) tthe rows or columns
%               of the loadings before rotation.
%   tolerance: default 1e-10
%   it_max:    maximum number of iterations, default 1000.
%
% OUTPUT_ARGS:
%
%   x : the new ratated loadings
%   r : the rotation matrix used.
%
% USAGE:
%
%   [new_loads, rotmax] = varimax(loadings, normalise, tolerance, it_max)
%
% HEADER_END
% 
%%%%%%%%%%%%%%%%%%
%
% 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, normalise, tolerance, it_max )
%
% where all but the loadings are optional.  rotmax is the rotation matrix used.
%
% normalise determines whether or not to normalise the rows or columns of
% the loadings before performing the rotation.  If normalise is true, then
% the rows are normalised by there individual lengths.  Otherwise, no
% normalisation is performed (default).  After rotation, the matrix is
% renormalised. Normalising over the rows corresponds to the Kaiser
% normalisation 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
  normalise = 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 normalise 
  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;

% Renormalise.

if normalise
  x = x .* rl;
end

if nargout > 1
  r = TT;
end