% PROGRAM correct_sigma2
%   Corrects the estimate of sigma^2 for sampling variation when 
%   census counts represent means of replicate samples


%************** USER-SUPPLIED INFORMATION *********************

% The matrix "samples" stores the counts from each sample; 
% it has censuses (years) as its columns and individual samples 
% as its rows; -1 is used as a place holder at the bottom of 
% columns corresponding to censuses with fewer samples
samples=[ 
 10 11  7  0 14  0  6 26 22 12  6 10 22 20 39 31 23 21 14 14 19;
  8  4  6 13  1 14 23 10 11 11 10  8 11 27 19 18 16 15 32 38 21;
  1 10 11 27 36 16 18 13 29  7  5  2 16  7 25 15 21 12 15 18  9;
  7  8  1  9 17 14  9 11 24 13  4  8 12 28  8 15 21 11 30 28 15;
  0 16 10  5 16 17  4 19 13 14 21 10 14  6 21 12  8  8  6 14 32;
  7  6  5  4  2 15  3 10 13  9 15 22  6 33 14 23  8 14 27 28 17;
  5 11 10 21 10 23 12  7 14 12 11 11 22 10 44  5 12 17 21 42 40;
  2 16 15  8  6  6 12 14 15 18  9  5  4  7 20 22 14  5 16 13 30;
  6  5 12 11 10  6 18  8 17  7 22 13 -1  5 15  7  8 29 25 16 23;
  6  6 -1 25 10 14 19  4 16  7 19  6 -1 10 16 18 14 15 11 26 18];
 
% To calculate confidence limits using the method of White (see
% below), users must supply the critical values of the chi-
% square distribution with degrees of freedom equal to the 
% number of censuses minus one; see Chapter 3 of Morris and Doak, 
% Quantitative Conservation Biology, for how to 
% compute these values.  The following are critical chi-square
% values for p=0.025 and p=0.975 with 20 degrees of freedom:
chi2crit=[34.16958143, 9.590772474];

% Set tolerance for computing corrected estimate and its 
% confidence limits using White's method; smaller tolerance means 
% more accuracy but a slower program
tolerance=1E-8;

%***************************************************************

% Set tolerance for the function fzero below
options=optimset('TolX',tolerance);

% The vector n stores the number of samples for each census
n=sum(samples>=0);
    
% q=number of censuses minus 1
q=size(samples,2) - 1;          

% Compute means and variances of the samples from each census
for t=1:q+1
    samplest=samples(1:n(t),t);
    Nbar(t)=mean(samplest);  
    Vs(t)=var(samplest);
end;
 
% Compute raw estimates of mu and sigma^2 using the conventional method
for t=1:q
   loglam(t)=log(Nbar(t+1)/Nbar(t));
end;
muest=mean(loglam);
disp('Raw estimate of mu:');
disp(muest);
s2raw=var(loglam);
disp('Raw estimate of sigma^2:');
disp(s2raw);

% Compute the component of total variance due to sampling variation
for t=1:q
    SampleVar(t)=Vs(t)/(n(t)*Nbar(t)^2) + ...
                    Vs(t+1)/(n(t+1)*Nbar(t+1)^2); 
end;
MeanVar=mean(SampleVar);  % mean sampling variation across years

% Compute the simple corrected estimate of sigma^2
s2corr=s2raw-MeanVar;
disp('Estimate of sigma^2 with simple correction for sampling variation');
disp(s2corr);

% Compute the corrected estimate of sigma^2 using the method of White,
% which weights by the sum of the environmental and sampling
% variance in each year;
% See G.C. White, Population viability analysis: data requirements 
% and essential analyses. Pages 288-331 in L. Boitani and T.K. Fuller (eds.), 
% Research Techniques in Animal Ecology: Controversies and Consequences. 
% Columbia University Press, New York.(2000)
Dev2=(loglam-muest).^2;
s2corrwhite=fzero(inline('sum(Dev2./(s2+SampleVar))-q+1', 's2', ...
    'Dev2','SampleVar', 'q'),s2raw,options,Dev2,SampleVar,q);
disp('Estimate of sigma^2 with White''s correction for sampling variation');
disp(s2corrwhite);

% Compute confidence interval for White's version of the corrected
% sigma^2;
% WARNING: Users may need to try different values of 'start' below to 
% get reasonable values for the confidence limits, as there may be more than
% one root to the equation being solved by the function fzero

start=0.001;
s2lower=fzero(inline('sum(Dev2./(s2+SampleVar))-chi2', 's2', ...
    'Dev2','SampleVar', 'chi2'),...
    start,options,Dev2,SampleVar,chi2crit(1));

start=0.1;
s2upper=fzero(inline('sum(Dev2./(s2+SampleVar))-chi2', 's2', ...
    'Dev2','SampleVar', 'chi2'),...
    start,options,Dev2,SampleVar,chi2crit(2));

CI=[s2lower s2upper];
disp('Confidence interval for White''s correction:');
disp(CI);