hierarch_like3.m

%negloglike calc for statespace model
function NNL = hierarch_like2(theta,rlist,nrlist,year,meanest,numzones,numspps,dtlist,tlist)

t = size(theta,1);
%break the param vector back up into components for like calc
p.a = theta(1:numzones*numspps,1);
p.S_g = theta(numzones*numspps+1:numzones*numspps*2,1);
p.err = theta(numzones*numspps*2+1:numzones*numspps*3,1);
p.R_t = theta(numzones*numspps*3+1:numzones*numspps*3+sum(tlist),1);
p.R_nt = theta(numzones*numspps*3+sum(tlist):t-4,1);
p.R_t_hat = theta(t-3); 
p.R_t_hat_err = theta(t-2);
p.R_nt_hat = theta(t-1);
p.R_nt_hat_err = theta(t);
datasize = size(meanest,1);

NNL=0;
y=zeros(datasize,1);
E = zeros(datasize,1);%vector of errors for use in normlike below


for i = 0:numzones-1
    for j=1:numspps

        if dtlist(i*numspps+j)==1

            %now do regression calcs in reserves for TARGETED SPPS
            y(rlist(i*numspps+j,2):rlist(i*numspps+j,3),1) = ...
                p.a(i*numspps+j) + year(rlist(i*numspps+j,2):rlist(i*numspps+j,3)) .* (p.S_g(i*numspps+j)+ p.R_t(sum(tlist(1:j)))); 
            %now spread out error terms across data bins
            E(rlist(i*numspps+j,2):rlist(i*numspps+j,3),1) = p.err(i*numspps+j);
            
            
            %now do regression calcs outside reserves for TARGETED SPPS
            y(nrlist(i*numspps+j,2):nrlist(i*numspps+j,3),1) = ...
                p.a(i*numspps+j) + year(nrlist(i*numspps+j,2):nrlist(i*numspps+j,3)).* (p.S_g(i*numspps+j)); 
            E(nrlist(i*numspps+j,2):nrlist(i*numspps+j,3),1) = p.err(i*numspps+j);
            
        else

             %now do regression calcs in reserves FOR NON-TARGETED SPPS
            y(rlist(i*numspps+j,2):rlist(i*numspps+j,3),1) = ...
                p.a(i*numspps+j) + year(rlist(i*numspps+j,2):rlist(i*numspps+j,3)).* (p.S_g(i*numspps+j)+ p.R_nt(sum(tlist(1:j)))); 
            E(rlist(i*numspps+j,2):rlist(i*numspps+j,3),1) = p.err(i*numspps+j);
            
            %now do regression calcs outside reserves FOR NON-TARGETED SPPS
            y(nrlist(i*numspps+j,2):nrlist(i*numspps+j,3),1) = ...
                p.a(i*numspps+j) + year(nrlist(i*numspps+j,2):nrlist(i*numspps+j,3)).* (p.S_g(i*numspps+j)); 
            E(nrlist(i*numspps+j,2):nrlist(i*numspps+j,3),1) = p.err(i*numspps+j);
            
        end
    end  
end

q=normpdf(meanest,y,E) + 1e-300;
rt=normpdf(p.R_t,ones(size(p.R_t,1),1)*p.R_t_hat,ones(size(p.R_t,1),1)*p.R_t_hat_err);
rnt=normpdf(p.R_nt,ones(size(p.R_nt,1),1)*p.R_nt_hat,ones(size(p.R_nt,1),1)*p.R_nt_hat_err);
NNL=sum(-log(q))+sum(-log(rt))+sum(-log(rnt));

%NNL = normlike([y,E],meanest); %calc regression specific NNL


       
%note that err may need to be fixed!!!!!!!!!! Currently shared by the same
%species in both reserve and non reserve sites. this is probably not right