ar1eof2.m

蒙托卡罗模拟奇异谱分析

 function [Er,Lr,Cr,g,a,mu2]=ar1eof2(x,C,E,R,k,V)
% AR1EOF2 - Calculate the red-noise basis for the noise + signal model
% Syntax: [Er,Lr,Cr,g,a,mu2]=ar1eof2(x,C,E,R,k,V);
%
% Computes the eigenfunctions of the expected lag-covariance matrix of
% a composite null-hypothesis (signal + AR(1) noise), using the method
% of Allen and Smith 1996.
%
% Note that Cr is the data + signal LCM and not the noise LCM -
% if the noise basis is desired, use [En,Vn,Cn]=ar1eof(g,a,mu2,M);. 
%
% Input:  x - the data vector.
%         C - the data lag-covariance matrix, from SSAEIG.
%         E - EOFs of x, from SSAEIG.
%         R - the matrix of RCs, from RC.
%         k - vector containing the indices of EOFs considered to be
%             signal: e. g. k=[8 9].
%         V - eigenvalues corresponding to E, from SSAEIG.
%
% Output: Er  - The red-noise + signal basis.
%         Lr  - the eigenvalues of Er.
%         Cr  - the process covariance matrix for noise + signal.
%          g  - estimate of the lag-one autocorrelation for the noise.
%          a  - estimate of the noise variance.
%         mu2 - estimated square on the mean for the noise.
%
% Written by Eric Breitenberger.      Version 3/14/97
% Please send comments and suggestions to [e-mail]       
%

M=length(E);
[g,a,mu2,c0,Q,W,K]=ar1sig(x,C,E,R,k);


% Normalization to avoid degeneracy or ill-conditioning:
% a bit different than that of A & S.

z=.9;                % arbitrary multiplier
Vmin=min(V(k));      % smallest signal eigenvalue
F=c0*Q*W*Q;
Vmax=max(eig(F));    % largest noise eigenvalue
nrm=z*Vmin/Vmax;

Cr=nrm*F+E*(eye(M)-K)*E'*C;


[Er,Lr]=eig(Cr);
[Lr,i]=sort(-diag(Lr));
Lr=-Lr';
Er=Er(:,i);

% Check to make sure that the known eigenfunctions are right:
tol=1e-12;
diff=polarize(Er(:,1:length(k)))-polarize(E(:,k));
if max(max(diff))>tol
  disp('AR1EOF2: Warning: Check eigenvectors of red-noise basis.')
end

% renormalize and resort:
Lr(length(k)+1:M)=Lr(length(k)+1:M)/nrm;

[Lr,i]=sort(-Lr);
Lr=-Lr';
Er=Er(:,i);


Cr=F+E*(eye(M)-K)*E'*C;