function NLL=logregA(params)
% This function gives the likelihood of a logistic regression 
% model for a metapopulation, given a particular set of 
% parameter values. 
% This function is called by a program (logregB) that fits the 
% maximum likelihood values for the parameters.
% This function is tailored to the particular example
% used in Morris and Doak Chapter 11 (Delphinium
% uliginosum data from the study of Harrison et al. 2000)

%*****************User-defined Parameters*********************
global FF1 FF2 AA sitenum Dists AAmx FF1mx NNL randtime ... 
mintime bestparams bestNLL;
ee = params(1:5)'; 			 % separate the extinction 
cc = [params(6:8)'; params(4:5)']; % and colonization parameters
totProbff1 = 0;
simcount = 0;
%************************************************************
% setting these parameters equal to 0 eliminates size effects:
ee(2)=0; cc(2)=0; ee(5)=0; cc(5)=0; % remove to include size

for reps = 1:randtime  % multiple simulations of ehe metapop'n
% On the first time through, it predicts the probability of
% seeing the second year of data -- other simulations are for 
% predictions of the first year occurrence pattern.
 FFnow = FF1; 	% set initial occupancies 
 for tt=1:(mintime +1) % time loop

%----------Start of basic calculations---------------------
 % matrix of occupancies for all neighbors:
 FFnowmx = FFnow(:,[ones(1,sitenum)])'; 
 % Calculation of neighbor effects for extinction probability:
 bit1=-ee(4)*Dists; % first part of extinction equation
 bit1(find(bit1>100))=100;	% check for values too big or small
 bit1(find(bit1<-100))=-100;
 connmxe = exp(bit1).*(AAmx.^(ee(5))).*FFnowmx;  
 connmxe = connmxe-diag(diag(connmxe)); % remove site’s self    
 connmxe(find(connmxe>20))=20; % check for too big or small
 connmxe(find(connmxe<(-20)))=-20;
 % summed, dist.-discounted areas of neighbors:
 conne = sum(connmxe')'; 

 % Calculate neighbor effects for colonization probability:
 bit1=-cc(4)*Dists;
 bit1(find(bit1>100))=100;
 bit1(find(bit1<-100))=-100;
 connmxc = exp(bit1).*(AAmx.^(cc(5))).*FFnowmx;
 connmxc = connmxc-diag(diag(connmxc));  
 connmxc(find(connmxc>20))=20;
 connmxc(find(connmxc<(-20)))=-20;
 % summed, dist.-discounted areas of neighbors:
 connc = sum(connmxc')'; 

 % make the linear fns ue and uc to go into logit. fns:
 uue = ee(1) + ee(2)*AA + ee(3)*conne;
 uuc = cc(1) + cc(2)*AA + cc(3)*connc;
 % initialize logistic equations:
 loge=zeros(sitenum,1); logc=zeros(sitenum,1);

 % find the predicted probabilities of ext and col:
 uue(find(uue >10))=10; uue(find(uue<-10))=-10;
 loge = exp(uue)./(exp(uue)+1);
 uuc(find(uuc >10))=10; uuc(find(uuc<-10))=-10;
 logc = exp(uuc)./(exp(uuc)+1);

 % calculate Negative Log-Likelihood for the observed
 % extinctions and colonizations:
 if ((tt == 1) & (reps == 1))
    LLec=FFnow.*( (1-FF2).*log(loge) +FF2.*log(1-loge) ) ...
    + (1-FFnow).*( (FF2).*log(logc) +(1-FF2).*log(1-logc) );
    NLLec = -sum(LLec); % neg LL of extinctions and col.'s
 end; % if 

 % randomly make next year's simulated occupancy:
 transprobs=rand(sitenum,1); 
 FFnew = FFnow.*(ceil(transprobs-loge)) + ...
(1-FFnow).*(ceil(logc-transprobs));
%----------End of basic calculations---------------------

 % calculate the probability of going from the simulated
 if tt == mintime +1 % occupancy to the initial real one:
    logProbff1= sum(FFnew.*( (1-FF1).*log(loge) + ...
FF1.*log(1-loge) ) + (1-FFnew).*( (FF1).*log(logc) ...
+(1-FF1).*log(1-logc)));
    totProbff1 = totProbff1 +exp(logProbff1);
    simcount = simcount+1;
 end; % tt == mintime + 1
 FFnow = FFnew; % advance the simulated occupancies
 end; % tt loop over the multi-year simulation
end; % reps loop

% negative log-likelihood of initial occupancy data
NLLint=-log(totProbff1/simcount); 
% to use both initial occupancy and ext/col data for fitting: 
% NLL= NLLec+NLLint; % final, total negative log-likelihood
% to use only ext/col data for fitting: 
NLL= NLLec; % final negative log-likelihood
disp(params') % displays each set of parameter values used
disp([NLLec,NLLint, NLL]); % displays the 3 NLL values
if NLL <bestNLL % keeps a running tab of best estimates found
    bestNLL = NLL % remove ';' on this line to show estimates
    bestparams=params'
end;