%IMPORTANT!
%******before import data are sorted as follows....
%first by unique spp number
%then by biogeozone
%then by reserve
%then by year
%The model is set up as follows:
%FIRST LEVEL-------------------------------------------------------------------------------------------------------
%MODEL:***** y~ intercept + abundance_index * (S_g + R_t +R_nt) + error ***
%where:
%S_g is the slope for a given species in a given geographic region (drawn
%from a normally distributed species specific hyperparameter)
%R_t is a species specific additive adjustment to the slope that is drawn
%from a normally distributed hyper parameter representing the reserve effect on harvested species. Activated by a dummy variable.
%R_nt is the species specific additive adjustment to the slope that is drawn from a normally distributed hyper parameter
%representing the reserve effect on non-harvested species. Activated by a dummy variable.
%both intercept and error terms are nuisance parameters and do not need hyperparameters. I suppose you could put in an
%error hyper parameter with an informative prior just to keep the model in check.
%SECOND LEVEL (hyper parameters all normally distributed)-----------------------------------------------------------
%S_g_hat -- These are slopes for a given species from which S_g's are drawn
%R_t_hat -- This is the reserve effect for targeted species from which R_t's across targeted species are drawn
%R_nt_hat -- This is the reserve effect for non-targed species from which
%R_nt's across non-targeted species are drawn
%--------------------------------------------------------------------------------------------------------------------------
%written by Brice X. Semmens (brice.semmens @ noaa.gov)
%last updated on 12/18/07
clear all;
%DECLARE MODEL PARAMETERS
nreps = 100000; % total number of evaluations in the MCMC
nburns =20000;% number of reps discarded before results start to be stored
nevery = 100; % period of replicates to be stored
gold = 0;%like stor
g1 = 0; %another like stor
prtcount=0;%print ticker (don't mess with!)
% read in data
data = load('output_selected_species.txt'); %bring in data file
%format of data is matrix, with columns as follows:
% spp | year | biogeographic_zone | reserve_status | abund_pdfest | targeted? | mean_reefscore
% break out data
spp = data(:,1);
year = data(:,2);
year = year -min(year)+1; %get rid of the big years reference.
zone = data(:,3);
reserve = data(:,4);
abdest = data(:,5);
target = data(:,6);
meanest = data(:,7);
numspps = size(unique(spp),1);%total number of spps
numzones = size(unique(zone),1);%total number of biogeographic zones
tep =unique(cat(2,spp,target),'rows'); %temp storage
spplist = tep(:,1);%list of species
tlist = tep(:,2);%target status of each species
clear tep;
zonelist = unique(zone);%list of zones
NNL=0; %storage param for neg log likelihood
dspplist= cat(1,spplist,spplist);%double up the spplist for use in regression loop
dtlist= cat(1,tlist,tlist);%double up the tlist for use in regression loop
%--------------------------------------------------------------------------
% set up 1st level parameter vectors (columns)
p.a = zeros(numzones*numspps, 1); %intercepts
p.S_g = ones(numzones*numspps, 1)*.1; %slopes
p.err = zeros(numzones*numspps,1);%error terms
p.R_t =ones(sum(tlist),1)*.01;
p.R_nt = ones(size(tlist,1)-sum(tlist),1)*.01;
%Set some inits and load data into structure fields
for i = 0:numzones-1
for j= 1:numspps
rlist(i*numspps+j,1)=i*numspps+j; %list of regressions (in reserves)
nrlist(i*numspps+j,1)=i*numspps+j; %list of regressions (outside reserves)
rlist(i*numspps+j,2)= min(intersect(find(reserve==1),intersect(find(spp==spplist(j)),find(zone==zonelist(i+1)))));%cell refs here and below
rlist(i*numspps+j,3)= max(intersect(find(reserve==1),intersect(find(spp==spplist(j)),find(zone==zonelist(i+1)))));
nrlist(i*numspps+j,2)= min(intersect(find(reserve==0),intersect(find(spp==spplist(j)),find(zone==zonelist(i+1)))));
nrlist(i*numspps+j,3)= max(intersect(find(reserve==0),intersect(find(spp==spplist(j)),find(zone==zonelist(i+1)))));
%now do regressions for each to get init values for p.a, p.err, p.S_g
%first reserve regressions...
[B,BINT,R,RINT,STATS] =...
regress(meanest(rlist(i*numspps+j,2):rlist(i*numspps+j,3),1),...
[ones(1+rlist(i*numspps+j,3)-rlist(i*numspps+j,2),1) year(rlist(i*numspps+j,2):rlist(i*numspps+j,3),1)]);
p.a(i*numspps+j) = B(1);
p.S_g(i*numspps+j) = B(2);
p.err(i*numspps+j) = STATS(4);
%now outside reserve regressions...
[B,BINT,R,RINT,STATS] =...
regress(meanest(nrlist(i*numspps+j,2):nrlist(i*numspps+j,3),1),...
[ones(1+nrlist(i*numspps+j,3)-nrlist(i*numspps+j,2),1) year(nrlist(i*numspps+j,2):nrlist(i*numspps+j,3),1)]);
%now get the average of these (min for error--MCMC should bring it up if need be)
p.a(i*numspps+j) = (p.a(i*numspps+j)+B(1))/2;
p.S_g(i*numspps+j) = (p.S_g(i*numspps+j)+B(2))/2;
p.err(i*numspps+j) = min(STATS(4), p.err(i*numspps+j));
end
end
% set up 2nd level parameters
p.S_g_hat = ones(numspps,1); %establishes and initializes here
p.S_g_hat_err = ones(numspps,1)*10; %establishes and initializes here
p.R_t_hat = .05;
p.R_t_hat_err = 1;
p.R_nt_hat = .05;
p.R_nt_hat_err = 1;
%-------------------------------------------------------------
%MCMC VECS FOR WORK
deltainit = .01; % 10% param jump size
ndelta = 50; % how long betwen updating deltas
theta = [p.a;p.S_g;p.err;p.R_t;p.R_nt;p.R_t_hat;p.R_t_hat_err; p.R_nt_hat;p.R_nt_hat_err]; %concats param vectors
npars = size(theta,1);
prtcount = 0; %ticker used to write param vectors at nevery intervals
%get MLE ests without hypers ... doesn't work. but keep for later reference.
OptArg=optimset('MaxFunEvals',100000,'MaxIter',100000,'TolX',1e-6,'TolFun',1e-6);
theta = fminsearch('hierarch_like3',theta,OptArg,rlist,nrlist,year,meanest,numzones,numspps,dtlist,tlist);
%get the initial jump size for each param in sigma
delta = theta * deltainit; %produces vector of jump sizes for each param
jchk(npars,1) = 0;
ichk(npars,1) = 0;
%now do a little parameter jumping!
repcounter=0;
gold = hierarch_like3(theta,rlist,nrlist,year,meanest,numzones,numspps,dtlist,tlist); %calculate the joint-like oval all params as a start
adjust = gold-350; %this used to make gold small enough to take the exp
gold=exp(-gold+adjust);
for irep = 1:nreps %how many reps in MCMC
tic
repcounter=repcounter+1;
%PICK A NEW VECTOR OF PARAMETERS
%-----------------------------------------------------------------------
for ipar = 1:npars %loop over params
%save old vals for spp-pair parameters
asave = theta(ipar);
%take the jump
theta(ipar) = theta(ipar) + ((rand - 0.5) * 2 * delta(ipar));
%re-evaluate the nlike after the jumps
g1 = hierarch_like3(theta,rlist,nrlist,year,meanest,numzones,numspps,dtlist,tlist);
g1=exp(-g1+adjust);
% if the new like is greater than the old
%or the ratio is > than a random number between 0 and 1 then accept
if (g1/gold) > rand %if g1 is bigger than gold, or that ratio bigger than a 0-1 unif. random number...
gold=-log(g1) + adjust; %the true NNL without adjustment
%now recalc the adjustment
adjust = gold-350; %this used to make gold small enough to take the exp
gold=exp(-gold+adjust);%keeps the new (adjusted) like
ichk(ipar) = ichk(ipar) + 1; %keeps track of numbers of accepted
else % don't accept
theta(ipar) = asave;
jchk(ipar) = jchk(ipar) + 1; %keeps track of number of rejections
end
end %ends ipar loop
%NOW ADJUST JUMP SIZES DURING THE BURN IN
%--------------------------------------------------------------------
if mod(irep,ndelta) == 0 & irep < nburns
delta = delta*1.05 .* ceil((ichk - jchk) /irep)... %adjust up if more accepts than rejects
+ delta*0.95 .* abs(floor((ichk - jchk -.01)/irep));%adjust down (the .01 accounts for when =)
delta
%now reset the vectors for keeping track of accepts and rejects
ichk(:,:) = 0; jchk(:,:) = 0;
save temp_hie_stor; %saves the workspace
end %ends burn adjust if statement
%WRITE RESULTS
%--------------------------------------------------------------------
if mod(irep,nevery) == 0 & irep > nburns
save temp_hie_stor; %saves the workspace
prtcount = prtcount + 1; %advance ticker
a_out(prtcount,:) = theta(1:numzones*numspps,1)';
S_g_out(prtcount,:) = theta(numzones*numspps+1:numzones*numspps*2,1)';
err_out(prtcount,:) = theta(numzones*numspps*2+1:numzones*numspps*3,1)';
R_t_out(prtcount,:) = theta(numzones*numspps*3+1:numzones*numspps*3+sum(tlist),1)';
R_nt_out(prtcount,:) = theta(numzones*numspps*3+sum(tlist):size(theta,1),1)';
p.R_t_hat_out(prtcount,1) = theta(size(theta,1)-3);
p.R_t_hat_err_out(prtcount,1) = theta(size(theta,1)-2);
p.R_nt_hat_out(prtcount,1) = theta(size(theta,1)-1);
p.R_nt_hat_err_out(prtcount,1) = theta(size(theta,1));
end
if mod(irep,50) == 0 %print reps to screen every 10 reps
irep
end
toc
end %end irep