R/bayesfun.R

In pttstability: Particle-Takens Stability

#' Parse parameters
#'
#' Takes in a vector of 3 or 6 parameters, and puts them into a list of the format expected by the particleFilterLL function.
#' @param param List of paramters, of length 2, 3, 4, 6, 7, or 8.
#' If 2, then in the order (obs1, proc1).
#' If 3, then in the order (obs1, proc1, proc2).
#' If 4, then in the order (obs1, obs2, proc1, proc2).
#' If 6, then in the order (obs1, proc1, pcol1, pcol2, det1, det2)
#' If 7, then in the order (obs1, proc1, proc2, pcol1, pcol2, det1, det2)
#' If 8, then in the order (obs1, obs2, proc1, proc2, pcol1, pcol2, det1, det2)
#' Note that if param is of length 2 or 3, then detparam  and colparam must be supplied. See obsfun0, procfun0, and detfun0 for more details.
#' @param colparam Optional vector of length two, including parameters for the colonization function.
#' @param detparam Optional vector of length two, including paramters for the deterministic function.
#' @keywords stability time-series particle filter
#' @return a formatted list of parameters
#' @export

parseparam0<-function(param, colparam=c(logit(0.2), log(0.1)), detparam=c(log(1.2),log(1))) {
  if(length(param)==2) {
    pars<-list(obs=c(param[1]),
               proc=c(param[2]),
               pcol=colparam[c(1:2)],
               det=detparam[c(1:2)])
  } else if(length(param)==3) {
    pars<-list(obs=c(param[1]),
               proc=c(param[2:3]),
               pcol=colparam[c(1:2)],
               det=detparam[c(1:2)])
  } else if(length(param)==4) {
    pars<-list(obs=c(param[1:2]),
               proc=c(param[3:4]),
               pcol=colparam[c(1:2)],
               det=detparam[c(1:2)])
  } else if(length(param)==6) {
    pars<-list(obs=c(param[1]),
               proc=c(param[2]),
               pcol=c(param[3:4]),
               det=c(param[5:6]))
  }  else if(length(param)==7) {
    pars<-list(obs=c(param[1]),
               proc=c(param[2:3]),
               pcol=c(param[4:5]),
               det=c(param[6:7]))
  } else if(length(param)==8) {
    pars<-list(obs=c(param[1:2]),
               proc=c(param[3:4]),
               pcol=c(param[5:6]),
               det=c(param[7:8]))
  } else {
    return("error: param must be of length 2, 3, 4, 6, 7, or 8")
  }

  return(pars)
}


#' Default likelihood function
#'
#' Calculates likelihood of vector y given parameter values in param, based on the particleFilterLL function.
#' @param param An unformatted vector of parameters, to be passed to parseparam function.
#' @param y A numeric vector of observed values, from which the likelihood of parameters and functions will be determined.
#' @param parseparam A function for transforming the vector param into a form that can be read by particleFilterLL. See particleFilterLL for details.
#' @param N Number of particles to simulate. Defaults to 1e3.
#' @param detfun A function that simulates deterministic dynamics, which takes in arguments sdet (parameters for deterministic model, taken from pars$proc), and xt, observed abundances at time t.
#' Returns estimated abundances at time t+1 based on deterministic function (either a parametric function or an EDM function). Defaults to detfun0.
#' @param edmdat A list including arguments to be passed to S_map_Sugihara1994 - see S_map_Sugihara1994 help file for details. Alternatively, the user can provide a matrix of pre-computed S-map coefficients, in element "smp_cf".
#' Default for edmdat is NULL, which implies that EDM will not be applied - instead, a detfun and pars$det must be included.
#' @param obsfun The observation error function to be used: defaults to obsfun0
#' @param procfun The process noise function to be used: defaults to procfun0
#' @param neff Should effective sample size be used to scale likelihood? Defaults to FALSE. TRUE uses automatic sample size, based on correlations in y. Otherwise, can be any positive number.
#' @param lowerbound Lower bound for log likelihood. Filter will be re-run if the value falls below this threshold. NOTE - this option may induce a bias
#' in the resulting likelihood (and subsequent parameter) estimates. Should only be set if the lower limit is indicative of filter failure (e.g. if all particles)
#' are degenerate. Defaults to (-Inf) - i.e. no lower limit.
#' @keywords stability time-series MCMC optimization
#' @return Log likelihood generated by particleFilterLL function
#' @export

likelihood0 <- function(param, y=y, parseparam=parseparam0, N=1e3, detfun=detfun0, edmdat=NULL, obsfun = obsfun0, procfun=procfun0, neff=FALSE, lowerbound=(-999)) {
  pars<-parseparam(param)

  tmp<-particleFilterLL(y = y, pars = pars, N=N, detfun = detfun, edmdat = edmdat, procfun = procfun)
  LL<-tmp$LL
  if(LL<=lowerbound) {
    LL<-lowerbound
  }

  if(isTRUE(neff)) {
    neff<-length(y)/(1+(length(y)-1)*(abs(cor(y[-length(y)], y[-1]))))
    LL<-LL*neff
  } else if(neff>0) {
    LL<-LL*neff
  }

  return(LL)
}


#' Default density function for prior
#'
#' Default density function, following the syntax for priors in the BayesianTools package. Uses
#' flat priors for all paramters, within the given interval. Density function integrates to one.
#' @param param a vector model parameters
#' @param minv a vector of minimum values for the interval
#' @param maxv a vector of maximum values for the interval
#' @keywords stability time-series MCMC optimization
#' @return returns log likelihood of parameters given priors.
#' @export

density_fun0 = function(param, minv, maxv){
  dsum<-0
  dwidth<-maxv-minv

  dsum<-0
  for(i in 1:length(param)) {
    if(param[i]<=maxv[i] & param[i]>=minv[i]) {
      dsum<-dsum+log(1/dwidth[i])
    } else {
      dsum<-(-Inf)
    }
  }

  return(dsum)
}


#' Default sampler function for prior
#'
#' Draws samples from a flat prior
#' @param n number of random draws to take from the priors
#' @param minv Vector of minimum values to return for each parameter
#' @param maxv Vector of maximum values to return for each parameter
#' @keywords stability time-series MCMC optimization
#' @return returns random draws from the priors
#' @import stats
#' @export

sampler_fun0 = function(n=1, minv, maxv){
  sampout<-matrix(nrow=n, ncol=length(minv))
  for(i in 1:length(minv)) {
    sampout[,i]<-runif(n, min = minv[i], max = maxv[i])
  }
  sampout
}


#' Default inverse transormation function
#'
#' Takes in a matrix, where each column represents a parameter. Returns parameters in untransformed space.
#' If length = 2, then in the order (obs1, proc1).
#' If 3, then in the order (obs1, proc1, proc2).
#' If 4, then in the order (obs1, obs2, proc1, proc2).
#' If 6, then in the order (obs1, proc1, pcol1, pcol2, det1, det2)
#' If 7, then in the order (obs1, proc1, proc2, pcol1, pcol2, det1, det2)
#' If 8, then in the order (obs1, obs2, proc1, proc2, pcol1, pcol2, det1, det2)

#' @param x an nxm matrix with
#' @keywords stability time-series MCMC optimization
#' @return returns back-transformed values of parameters
#' @export

inv_fun0<-function(x) {
  if(ncol(x)==2) {
    cbind(exp(x[,1]), exp(x[,2]))
  } else if(ncol(x)==3) {
    cbind(exp(x[,1]), exp(x[,2]), exp(x[,3]))
  } else if(ncol(x)==4) {
    cbind(exp(x[,1]), exp(x[,2]), exp(x[,3]), exp(x[,4]))
  } else if(ncol(x)==6) {
    cbind(exp(x[,1]), exp(x[,2]), ilogit(x[,3]), exp(x[,4]), exp(x[,5]), exp(x[,6]))
  } else if(ncol(x)==7) {
    cbind(exp(x[,1]), exp(x[,2]), exp(x[,3]), ilogit(x[,4]), exp(x[,5]), exp(x[,6]), exp(x[,7]))
  } else if(ncol(x)==8) {
      cbind(exp(x[,1]), exp(x[,2]), exp(x[,3]), exp(x[,4]), ilogit(x[,5]), exp(x[,6]), exp(x[,7]), exp(x[,8]))
  } else {
    return("error: x must be a matrix with 2, 3, 4, 6, 7, or 8 columns.")
  }
}