Nothing
## simulation of breaking time points bbs
## for a 2 state telegraph process
sim1mm.times.bbz <- function(s, lam1, lam2) {
tsum <- 0
tt <- NULL
state <- NULL
while (TRUE) {
tnew <- rexp(1, lam1) ## starting from state 1
tsum <- tsum + tnew
tt <- c(tt, tsum)
state <- c(state, 0)
if (tsum > s) break
tnew <- rexp(1, lam2)
tsum <- tsum + tnew
tt <- c(tt, tsum)
state <- c(state, 1)
if (tsum > s) break
}
cbind(tt, state)
}
## simulation of a realization given breaking times
sim1mm.bbz <- function(s, sigma1, sigma2, time, brtimes, t0moving=TRUE) {
#### time: time points in [0, s]
tt <- sort(unique(c(time, brtimes)))
nt <- length(tt)
nb <- length(brtimes)
x <- rep(NA, nt)
x[1] <- 0
status <- as.integer(t0moving)
tend <- brtimes[1]
j <- 1
for (i in 2:nt) {
if (tt[i] <= tend) { ## status unchanged
if (status == 1) ## moving1
x[i] <- x[i - 1] + rnorm(1, sd = sigma1 * sqrt(tt[i] - tt[i - 1]))
else ## moving2
x[i] <- x[i - 1] + rnorm(1, sd = sigma2 * sqrt(tt[i] - tt[i - 1]))
}
if (tt[i] == tend) {
status <- 1 - status ## switch status
j <- j + 1
tend <- brtimes[j]
}
}
x[tt %in% time]
}
## simulation a moving-moving path given a grid time
#' Sampling from a Moving-Moving Process with 2 Embedded Brownian Motion
#'
#' A moving-moving process consists of two states: moving (large) and moving (small).
#' The transition between the two states is modeled by an alternating
#' renewal process, with exponentially distributed duration. An animal
#' moves according to two Brownian motions with different volatility parameters.
#'
#' @param time time points at which observations are to be simulated
#' @param lamM1 rate parameter of the exponential duration while moving1
#' @param lamM2 rate parameter of the exponential duration while moving2
#' @param sigma1 volatility parameter of the Brownian motion while moving1
#' @param sigma2 volatility parameter of the Brownian motion while moving2
#' @param s0 the state at time 0, must be one of "m1" or "m2", for moving1 and
#' moving2, respectively
#' @param dim (integer) dimension of the Brownian motion
#'
#' @return
#' A \code{data.frame} whose first column is the time points and whose
#' other columns are coordinates of the locations.
#' @references
#' Yan, J., Chen, Y., Lawrence-Apfel, K., Ortega, I. M., Pozdnyakov, V.,
#' Williams, S., and Meyer, T. (2014) A moving-resting process with an
#' embedded Brownian motion for animal movements.
#' Population Ecology. 56(2): 401--415.
#'
#' Pozdnyakov, V., Elbroch, L., Labarga, A., Meyer, T., and Yan, J.
#' (2017) Discretely observed Brownian motion governed by telegraph
#' process: estimation. Methodology and Computing in Applied Probability.
#' doi:10.1007/s11009-017-9547-6.
#'
#' @examples
#' tgrid <- seq(0, 100, length=100)
#'
#' dat <- rMM(tgrid, 1, 0.1, 1, 0.1, "m1")
#' plot(dat[,1], dat[,2], xlab="t", ylab="X(t)", type='l')
#'
#' @export
rMM <- function(time, lamM1, lamM2, sigma1, sigma2, s0, dim = 2) {
time <- time - time[1]
stopifnot(s0 %in% c("m1", "m2"))
t0moving <- (s0 == "m1")
lam1 <- if (t0moving) lamM1 else lamM2
lam2 <- if (t0moving) lamM2 else lamM1
timeIND <- 0
if (length(time) == 1) {
timeIND <- 1
time <- c(0, time)
}
tmax <- time[length(time)]
brtimes <- sim1mm.times.bbz(tmax, lam1, lam2)
coord <- replicate(dim, sim1mm.bbz(tmax, sigma1, sigma2, time, brtimes[, 1], t0moving))
stateresult <- simmr.state(time, brtimes, t0moving = t0moving)
if (timeIND == 1) {
return(data.frame(time = time, coord)[-1, ])
}
data.frame(time = time, coord)
}
#' Fit a Moving-Moving Model with 2 Embedded Brownian Motion
#'
#' Fit a Moving-Moving Model with 2 Embedded Brownian Motion with animal
#' movement data at discretely observation times by maximizing a full
#' likelihood constructed from the marginal density of increment.
#' 'estVarMM' uses parametric bootstrap to obtain variance matrix
#' of estimators from 'fitMM'.
#'
#' @param data a data.frame whose first column is the observation time, and other
#' columns are location coordinates.
#' @param start starting value of the model, a vector of four components
#' in the order of rate for moving1, rate for moving2,
#' and volatility1(larger), volatility2(smaller).
#' @param logtr logical, if TRUE parameters are estimated on the log scale.
#' @param method the method argument to feed \code{optim}.
#' @param optim.control a list of control to be passed to \code{optim}.
#' @param integrControl a list of control parameters for the \code{integrate}
#' function: rel.tol, abs.tol, subdivision.
#'
#' @return
#' a list of the following components:
#' \item{estimate}{the esimated parameter vector}
#' \item{loglik}{maximized loglikelihood or composite loglikelihood
#' evaluated at the estimate}
#' \item{convergence}{convergence code from \code{optim}}
#' @references
#' Yan, J., Chen, Y., Lawrence-Apfel, K., Ortega, I. M., Pozdnyakov, V.,
#' Williams, S., and Meyer, T. (2014) A moving-resting process with an
#' embedded Brownian motion for animal movements.
#' Population Ecology. 56(2): 401--415.
#'
#' Pozdnyakov, V., Elbroch, L., Labarga, A., Meyer, T., and Yan, J.
#' (2017) Discretely observed Brownian motion governed by telegraph
#' process: estimation. Methodology and Computing in Applied Probability.
#' doi:10.1007/s11009-017-9547-6.
#'
#' @examples
#' \dontrun{
#' ## time consuming example
#' tgrid <- seq(0, 100, length=100)
#' set.seed(123)
#' dat <- rMM(tgrid, 1, 0.1, 1, 0.1, "m1")
#'
#' ## fit whole dataset to the MR model
#' fit <- fitMM(dat, start=c(1, 0.1, 1, 0.1))
#' fit
#'
#' var <- estVarMM(fit$estimate, dat, nBS = 10, numThreads = 6)
#' var
#' }
#' @export
fitMM <- function(data, start,
logtr = FALSE,
method = "Nelder-Mead",
optim.control = list(),
integrControl = integr.control()) {
## normal process
if (!is.matrix(data)) data <- as.matrix(data)
dinc <- apply(data, 2, diff)
integrControl <- unlist(integrControl)
fit <- optim(start, nllk_inc_mm, data = dinc, method = method,
control = optim.control,
integrControl = integrControl,
logtr = logtr)
## ## get variance estimate
## varest <- matrix(NA_real_, 3, 3)
estimate <- if (logtr) exp(fit$par) else fit$par
## if (likelihood == "full") {
## ## always do this without log transformation
## hess <- tryCatch(numDeriv::hessian(
## objfun, estimate, data = dinc,
## integrControl = integrControl, logtr = FALSE),
## error = function(e) e)
## if (!is(hess, "error")) varest <- solve(hess)
## ## not -hess because objfun is negative llk already
## }
result <- list(estimate = estimate,
#varest = varest,
loglik = -fit$value,
convergence = fit$convergence,
data = data)
attr(result, "class") <- "smam_mm"
return(result)
}
#' @param est_theta estimators of MRME model
#' @param data data used to process estimation
#' @param nBS number of bootstrap.
#' @param detailBS whether or not output estimation results during bootstrap,
#' which can be used to generate bootstrap CI.
#' @param numThreads the number of threads for parallel computation. If its value
#' is greater than 1, then parallel computation will be processed. Otherwise,
#' serial computation will be processed.
#' @param integrControl a list of control parameters for the \code{integrate}
#' function: rel.tol, abs.tol, subdivision.
#' @rdname fitMM
#' @export
estVarMM <- function(est_theta, data, nBS, detailBS = FALSE,
numThreads = 1,
integrControl = integr.control()) {
tgrid <- data[, 1]
dim <- ncol(data) - 1
lam1 <- est_theta[1]
lam2 <- est_theta[2]
sig1 <- est_theta[3]
sig2 <- est_theta[4]
p_1 <- 1/lam1/(1/lam1 + 1/lam2)
p_2 <- 1 - p_1
if (numThreads <= 1) {
result <- matrix(NA, ncol = 4, nrow = nBS)
for (i in seq_len(nBS)) {
start_state <- sample(c("m1", "m2"), size = 1, prob = c(p_1, p_2))
datBS <- rMM(tgrid, lam1, lam2, sig1, sig2, s0 = start_state, dim = dim)
result[i, ] <- fitMM(datBS, start = est_theta,
integrControl = integrControl)$estimate
}
} else {
## create parallel backend
cl = parallel::makeCluster(numThreads)
on.exit(close(pb), add = TRUE)
on.exit(parallel::stopCluster(cl), add = TRUE)
doSNOW::registerDoSNOW(cl)
pb <- utils::txtProgressBar(max = nBS, style = 3)
progress <- function(n) utils::setTxtProgressBar(pb, n)
opts <- list(progress = progress)
i = 1
result <- foreach(i = 1:nBS, .combine = rbind, .options.snow = opts) %dopar% {
start_state <- sample(c("m1", "m2"), size = 1, prob = c(p_1, p_2))
datBS <- rMM(tgrid, lam1, lam2, sig1, sig2, s0 = start_state, dim = dim)
fit <- fitMM(datBS, start = est_theta, integrControl = integrControl)$estimate
fit
}
}
if (detailBS) {
return(list(cov = cov(na.omit(result)),
BS_detail = result))
} else {
return(cov(na.omit(result)))
}
}
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.