% Program White.m 
% This program finds the best estimate of environmental variance
% for annually collected vital rates using White's (2000) method. 
% To use the method, you must have an estimate of the mean and 
% variance in each rate for each year (the program will do the 
% analysis for each of many different rates or classes).
% This program does a brute force search for the best variance
% estimates from a very large number of  values.

%*****************Simulation Parameters***********************
% Copy in the basic data to use: mean and sampling variance 
% estimate for each year for several different classes.
% The data below are for mule deer fawn survival, used in 
% White (2000). The data columns are: 
% Class identifier, Year identifier, 
% Mean survival rate, and Within-year-sampling-variance:
rates=...
[1	81	0.326	0.0047773
1	82	0.333	0.0019493
1	83	0.0424	0.0003439
1	84	0.179	0.0013879
1	85	0.381	0.001521
1	86	0.379	0.0014617
1	87	0.129	0.0009706];

times = 7; 
% How many time periods (must be identical for all rates)
classes = 1; % how many classes or separate rates are there?
minvar = 0.0; % minimum estimate of true environmental variation
maxvar = 0.25; % Maximum possible variance estimate (0.25 for 
% survival rates, but it will often be much smaller)
trys = 10000; % how many variance estimates to try at one time
%**************************************************************

for class = 1:classes %looping through the different classes
  minrow = (class-1)*times +1; 
     % find the min and max rows of the data matrices to use
	 maxrow =  class*times; 
  data = rates(minrow:maxrow,:); % fetch the data to use
  meanvar = linspace(minvar,maxvar,trys); 
% generate a set of guesses for true env. variance
  meanrate = mean(data(:,3)); 
% calculate the grand mean rate over years 
    
  for ii = 1:length(meanvar) % a loop to calculate the sum
% that should equal 1 (see White 2000 for details)
     sumup(ii) = ...
sum( ((data(:,3)-meanrate).^2)./( (data(:,4) ... 
+ meanvar(ii)).*(times-1)) );
  end;  % ii loop
	
 [bestdiff,jj] = min(abs(sumup - 1)); %find the best sum:
% the one with the minimum difference from one
  estvar = meanvar(jj); 
% then find the corresponding best env. variance estimate

  % below, find the chi^2 values for 95% confidence limits;
  % If you don't have the Statistics Toolbox to use chi2inv,
  % look these up in a table of chi-square values.
  lowChi=chi2inv(0.975,(times-1)); 
  hiChi = chi2inv(0.025,(times-1));
  % below, find the differences from the correct Chi^2 value
  CLdiffslow=abs(sumup*(times-1)-lowChi);
  CLdiffshi=abs(sumup*(times-1)-hiChi); 
  %find differences of sums from the correct chi-square values
  [bestdifflow,jjlow] = min(CLdiffslow); 
  bestCLlow = meanvar(jjlow); % best lower CL estimate
  [bestdiffhi,jjhi] = min(CLdiffshi); 
  bestCLhi = meanvar(jjhi);% best upper CL estimate

% lines below display the best variance estimate and CLs for it
disp('for class or rate =');
disp(class);
disp('best variance estimate and lower and upper CLs =');
disp([estvar, bestCLlow,bestCLhi]);
disp('differences of rt. side of equ. 8.2 and one, ');
disp('for best, lower, and upper conf. limit estimates =');
disp([bestdiff, bestdifflow,bestdiffhi]);

% Using the estimates above, you should narrow the bounds 
% on possible variances (minvar, maxvarand rerun 
% to see if the estimates or the differences from one
% change substantially. If not, then the estimated
% variance values are good. 
end; % class loop