Nothing
#' 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.")
}
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.