# Data are in a table BCIdata.rdata. They consist of survival rates and wood density for individual species.
# Each row is a species, Nstart the number of individuals at the outset, S the number alive after time years.

# Wood density is wsg (wood specific gravity).

# The problem is estimating whether mortality rates are affected by wood density, that is, species
# with denser ought to have higher survival.

# Consider that each species j has mortality rate m[j], defined as -log(theta), where theta is the survival
# probability. A sampling estimate of m is -log(S/N). The mortality rates m[j] are distributed log-normally.
# with parameters logmean=a+b*w and logsd=sigma. (Those are the log-normal parameters used in R; the mean
# and SD of the log of observations.) The formula a+b*w defines the effect of wood-density (w) on mortality, that is,
# the mortality parameter varies linearly with density. 

# The number of survivors, S of Nstart, is observed with binomial error on the survival rate theta=exp(-m).

# The Gibbs sampler allows the 3 key parameters -- a, b, sigma -- to be fitted, but also works with a mortality
# parameter for every species. Ie, it fits the mortality rate for every species. The number of parameters fitted
# is thus 3 plus the number of species. The species means are called spmean in the main function.

# The data table also has growth data as well, with growth the mean diameter increment per year (of trees 10-49 mm
# in diameter only), and Ngrow the number of individuals used in calculating this mean. A growth model is given
# in WoodDensGrowthBayesWorkshop.r.

# A more difficult exercise would be to allow the parameter a to vary across taxonomic families (which
# are included in the table). This requires parameters to describe a log-normal distribution of a across families,
# and a second set parameter of parameters to describe the variation of mortality rates m within families.



# The main function sets up the data and loops through a Gibbs sampler of the specified number of steps. The data
# table must have columns Nstart, S, time, and wsg.
mortality.density=function(data,steps,showstep=25)
{
 data=subset(data,!is.na(lg)&Nstart>0&!is.na(wsg))
 nospp=dim(data)[1]
 mortparam=matrix(nrow=steps,ncol=3)
 spmean=matrix(nrow=steps,ncol=nospp)

 mortparam[1,1]=(-3)
 mortparam[1,2]=0
 mortparam[1,3]=1
 stepsize=rep(1,3)

 spmean[1,]=rep(.05,nospp)
 spstep=rep(0.01,nospp)

 llike=rep(NA,steps)

 for(i in 2:steps)
  {
   for(j in 1:3)
    {
     param=mortparam[i-1,]
     testparam=param[j]
     if(j>1) param[1:(j-1)]=mortparam[i,1:(j-1)]

     metropResult=metrop1step(func=llike.mortparam,start.param=testparam,param=param,whichtest=j,scale.param=stepsize[j],
                              spmean=spmean[i-1,],w=data$wsg,adjust=1.01,target=0.25)
     mortparam[i,j]=metropResult[1]
     stepsize[j]=metropResult[2]
    }

   for(j in 1:nospp)
    {
     metropResult=metrop1step(func=llike.mortspmean,start.param=spmean[i-1,j],param=param,scale.param=spstep[j],
                              w=data$wsg[j],N=data$Nstart[j],S=data$S[j],time=data$time[j],adjust=1.01,target=0.25)
     spmean[i,j]=metropResult[1]
     spstep[j]=metropResult[2]
    }

   llike[i]=metropResult[4]

   if(i%%showstep==0)
     cat(i, ": a ", round(mortparam[i,1],3), " b ", round(mortparam[i,2],3), " sigma ", round(mortparam[i,3],3),
         " llike ", round(llike[i],3),"\n")
  }

 return(list(mort=mortparam,spmean=spmean,llike=llike))
}


# The likelihood function for the main parameters, a, b, and sigma. In a Gibbs sampler, only one of the three is
# updated at a time, and the first argument of the function always must be the one. The argument whichtest
# indicates which of the the three is being updated by Metropolis, and the value to be tested is testparam. The other
# two parameters do not change. This method allows the same function to be used to update each of the 3 parameters (otherwise,
# 3 different functions would be needed).
llike.mortparam=function(testparam,param,whichtest,spmean,w)
{
 param[whichtest]=testparam
 a=param[1]
 b=param[2]
 sigma=param[3]

 if(sigma<=0) return(-Inf)
# browser()

 logmeanmort=a+b*w
 llike=dlnorm(spmean,meanlog=logmeanmort,sdlog=sigma,log=TRUE)

 return(sum(llike))
}


# The likelihood function for the mean mortality rates of every species, spmean. It works with a single
# species at a time, updating its estimated mean mortality rate. The likelihood requires two sections: a binomial,
# which is the probability of seeing S survivors of N individuals over time years; and a community-wide assessment
# of how likely the mortality rate is given the hyperdistribution (the distribution of mortality rates across
# the community.
llike.mortspmean=function(spmean,param,w,N,S,time)
{
 if(spmean<=0) return(-Inf)

 a=param[1]
 b=param[2]
 sigma=param[3]

 logmeanmort=a+b*w
 llike1=dlnorm(spmean,meanlog=logmeanmort,sdlog=sigma,log=TRUE)

 theta=exp(-spmean*time)
 llike2=dbinom(S,size=N,prob=theta,log=TRUE)

 return(llike1+llike2)
}