chi2conf.m

蒙托卡罗模拟奇异谱分析

 function [c,L]=chi2conf(E,Cn,N,p)
% CHI2CONF - calculate approximate confidence limits for an eigenspectrum
% Syntax: [c,L]=chi2conf(E,Cn,N,p);
%
% CHI2CONF implements the method described in the Appendix of Allen and 
% Smith 1996 for estimating confidence limits without resorting to
% Monte Carlo testing. The algorithm assumes known or Gaussian noise,
% (known noise covariance matrix), and approximately sinusoidal EOFs.
% It is thus best suited for use in testing against the red-noise basis
% rather than the signal basis.
%
% Input: E  - eigenbasis (normally from SSAEIG).
%
%        Cn - expected noise covariance matrix, possibly from AR1EOF.  
%
%        N  - number of points in the data series.
%
%        p  - confidence limit, in a two-tailed form, e. g. for p=.95 the
%             confidence limits will be at 2.5% and 97.5%. The default is p=.95.
%             p may be a vector, e. g. p=[.9 .95 .99].
%
% Output: c - An M by 2*length(p) matrix containing upper and 
%             lower confidence limits.
%         L - The eigenspectrum (projection of the noise LCM on the data basis).
%
% Written by Eric Breitenberger.      Version 5/22/96
% Please send comments and suggestions to [e-mail]       
%

if nargin==3, p=.95; end
p=p(:)';
M=size(E,1);

L=diag(E'*Cn*E);

nu=3*N/M; % Approximate degrees of freedom

% WILHIL takes single-tailed p, rather than 
% double-tailed, so convert: (e. g. p=.95 => p=[.025 .975])

p=[ (1-p)/2 1-(1-p)/2 ];

% Rearrange p so that the confidence limits are paired,
% i. e. the first two elements are for p(1), etc.

index=1:length(p)/2:length(p);
for i=2:length(p)/2
  index=[index i:length(p)/2:length(p)];
end

p=p(index);

[chi2]=wilhil(nu,p);

c=(L*ones(1,length(p))).*(ones(M,1)*chi2)/nu;