function f95 = eofsig_montecarlo(M,N,nsamples)
% function f95 = eofsig_montecarlo(M,N)
%
% Compute the confidence limits for significant eigenvalues in an EOF
% decomposition by generating a set of random data sets of dimension M x N
% and finding the percent variance (f95) in each mode for which there would
% be only a 5% probability of this occuring by chance.
%
% For EOFs of data sets that have significant autocorrelation (say in
% the time dimension of the data matirx), it is recommended that N be
% reduced to Nstar = N/autocorrelation_scale to reflect the number of
% independent observations of the data in each row of the data matrix. See
% Joliffe (2002) chapter 6.
% Joliffe, I., 2002, Principal Component Analysis, 2nd ed., 502pp,
% Springer.

if nargin < 3
  nsamples = 100;
end

% nsamples realizations for Monte Carlo
for i=1:nsamples
  if mod(i,round(nsamples/10)) == 0
    disp(['Doing sample ' int2str(i)])
  end
  % create random matrix of Normal(0,1) distributed values
  Dmc = randn(M,N);
  [Ur,Sr,Vr] = svd(Dmc,0);
  % fraction of variance
  fraction = diag(Sr*Sr')/trace(Sr*Sr');
  % preallocate for speed
  if i==1
    F = ones([length(fraction) nsamples]);
  end
  % save this sample
  F(:,i) = fraction;
end

% Now each column of F is the fraction of variance explained by the EOF
% decomposition (via SVD) for a dataset that is completely random. Any 
% skill at describing the variability with coherent patterns (U) is by pure
% chance. Accordingly, the percent variance explained decreases slowly 
% with EOF mode number (down each column of F) but inevitably some modes
% explain more variance that others. 
%
%   plot(F)  % to see this gradual ordered decline in variance explained
%   load matlab.mat nlats nlons      % to get lon/lat for plotting
%   pcolor(reshape(Ur(:,1),[nlats nlons]))
%
% The next step is to find the 95-percentile value from the Monte Carlo 
% set. This will indicate the level of fractional variance explained by 
% EOFs that could occur by chance for a set of data that actually had no
% deterministic patterns in it.
F = sort(F,2);
f95 = F(:,round(95/100*nsamples));