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.
%
% Note that 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 (20020 chapter 6. 

if nargin < 3
   nsamples = 100;
end

% nsamples realizations for Monte Carlo
for i=1:nsamples
   disp([int2str(i)])
   % create random matrix of Normal(0,1) distributed values
   Dmc = randn(M,N);
   [U,S,V] = svd(Dmc,0);
   % fraction of variance
   fraction = diag(S*S')/trace(S*S');
   % preallocate for speed
   if i==1
      F = ones([length(fraction) nsamples]);
   end
   % save this sample
   F(:,i) = fraction;
end
F = sort(F,2);
f95 = F(:,round(100*95/nsamples),:);