% Program Kendall.m
% This program finds the best estimates of mean and environmental 
% variance for beta-binomial vital rates (see Kendall 1998),  
% using a brute force search for the best adjusted
% estimates from a very large number of combinations of different
% possible mean and variance values.
% Note that it will deliver messages of 'divide by zero,'
% indicating use of very low values for variance but not  
% a malfunction.  
clear all; % clears all preexisting variables; this speeds up
% the program and can prevent errors

%*****************User-defined Parameters***********************
% First, copy in the basic data to use. Here we show data for 
% only the last three desert tortoise growth rates (g4	
% g5, g6), coded as rates (4 5 6). 
% Enter the data in five columns: Rate identifier, 
% Year, Total number of starting individuals, 	Number
% growing (or surviving):
grow=...
[4	1	17	8	;
4	2	15	1	;
4	3	7	0	;
5	1	22	5	;
5	2	19	5	;
5	3	4	0	;
6	1	32	2	;
6	2	31	1	;
6	3	10	0	];

rates = grow; 	% rates is the data set used below (you
% could input several different data sets above here and 
% this allows a choice of which to use).			
times = 3; % How many time periods for each rate (must be 
% the same for all rates for the program as written
classes = 3; % how many different rates or classes for the 
% rates are there?
grades = 1000; % the number of different levels of means and
% variances to try; more levels is better: 1000 is quite good
% and doesn't take too long
maxvar = 0.2; % the maximum variance to search over. The maximum 
% ever possible  is 0.25, but a reasonable maximum will often be 
% much smaller. Searching a narrower range will improve the 
% accuracy of the answer.
minvar = 0.00001; % The minimum variance to search. You can't set
% this minimum to zero, as this is un-computable. 
maxmean = 1; % maximum limit on the mean values to search
minmean = 0.010; % minimum limit on the mean values to search
%*************************************************************

results = []; results2 = [];
means = repmat(linspace(minmean,maxmean,grades),grades,1);			
% makes all the sets of mean values to search over
vars =  repmat(linspace(minvar,maxvar,grades)',1,grades);  
% makes all the sets of variance values to search over

aa = means.*( ((means.*(1-means))./vars) - 1);			
bb = (1 - means).*( ((means.*(1-means))./vars) - 1);
% using the means and variances to compute the a and b
% (aa and bb) parameters of a beta distribution

aa(find(aa<=0)) = NaN;  
bb(find(bb<=0)) = NaN;
% eliminate impossible combinations of mean and variance values

for class = 1:classes % loop through each rate or class
      % find the min and max rows of the data matrices to use
	minrow = (class-1)*times +1; 
	maxrow =  (class)*times; 
	data = rates(minrow:maxrow,:); % fetch the data to use
	estmn = mean(data(:,4)./data(:,3));  % raw estimate of mean
	estvar = var(data(:,4)./data(:,3));  % raw est. of variance
	loglikli = zeros(grades,grades); %initialize -loglikelihoods
% As is often done, we use the negative of log liklihoods,
% or -log likelihoods, in the following computations.
	for tt = 1:times	 % calculate -loglikelihood for each year
		newlogL = -log(beta( aa + data(tt,4),data(tt,3) ...
 			- data(tt,4)+ bb)./beta(aa,bb)); 
		 % this is uses the -log-likelihood formula from Kendall 
		loglikli = newlogL + loglikli; 
		 % add up the -log-likelihoods for each year
	end; % tt loop

disp('Class is');, disp(class); %track where the program is

% the lines below summarize the results 
minLL = min(min(loglikli)); % what is the best log-liklihood?
[minvars,ii] = min(loglikli);  % finding best var
[minii, jj] = min(minvars);
MLEvar = vars(ii(jj),1);
[minvars,ii] = min(loglikli');  % finding best mean
[minii, jj] = min(minvars);
MLEmean = means(1,ii(jj));

clLL = loglikli - minLL;  % differences from best -LL 
hivar = max(max(vars(find(clLL<3.0))));%find conf. limits for var
lowvar = min(min(vars(find(clLL<3.0))));
himean = max(max(means(find(clLL<3.0)))); % conf limits for mean
lowmean = min(min(means(find(clLL<3.0))));

% perform correction for max. likelihood estimate of variance
% (see Kendall 1998 for explanation)
corrMLEvar = MLEvar*times/(times-1); 
corrlowvar = lowvar*times/(times-1);
corrhivar = hivar*times/(times-1);

% storing the results
results = ...
 [results; data(1,1),estmn,MLEmean,estvar,MLEvar,corrMLEvar];
results2 = ...
 [results2; data(1,1),lowmean,himean,corrlowvar,corrhivar];

end; % class

% now, display the results:
disp('Basic results are...')
disp('class, estimated mean, MLE mean, estimated var,');
disp('MLE variance, corrected MLE variance');
disp(results)
disp('the likelihood 95% conf. limits for the mean');
disp('and corrected variance are:')
disp('class, low and hi mean limits, low and hi variance limits')
disp(results2)

% finally, make a 3d mesh diagram for relative –LL for
% the last vital rate class, using values only within the 
% limits of the 95% c.l.'s for variance and mean

clLL = -(loglikli - minLL) + 3.0;		
%standardize and reverse the orientation of the -2logL's
clLL(find(clLL<0)) = 0;
vars = vars.*(times/(times-1));
mesh(means,vars,clLL);
xlabel('Mean');
ylabel('Variance (corrected)');
zlabel('2Log-Likelihood (only values within 95% C.L.s');
axis([0,1,0,maxvar]);	%AXIS([XMIN XMAX YMIN YMAX])
view(50,20);   %VIEW([AZ,EL]) AZ is the azimuth and EL is 
 % the vertical elevation (both in degrees)