R/simulate.R In SimoneHermann/BaPreStoPro: Bayesian Prediction of Stochastic Processes

Defines functions simNdrawSDE

########
#' Simulation of diffusion process
#'
#' @description Simulation of a stochastic process
#'   \eqn{dY_t = b(\phi,t,Y_t)dt + \gamma \widetilde{s}(t,Y_t)dW_t}.
#' @param object class object of parameters: "Diffusion"
#' @param nsim number of trajectories to simulate. Default is 1.
#' @param seed optional: seed number for random number generator
#' @param t vector of time points
#' @param y0 starting point of the process
#' @param mw mesh width for finer Euler approximation to simulate time-continuity
#' @param plot.series logical(1), if TRUE, simulated series are depicted grafically
#' @examples
#' model <- set.to.class("Diffusion", parameter = list(phi = 0.5, gamma2 = 0.01))
#' t <- seq(0, 1, by = 0.01)
#' data <- simulate(model, t = t, y0 = 0.5, plot.series = TRUE)
#' @export
setMethod(f = "simulate", signature = "Diffusion",
definition = function(object, nsim = 1, seed = NULL, t, y0, mw = 1, plot.series = TRUE) {
set.seed(seed)
if(nsim > 1){
result <- matrix(0, nsim, length(t))
for(i in 1:nsim){
result[i,] <- drawSDE(object@phi, object@gamma2, t, object@b.fun, y0.fun = function(phi, t) y0,
sT = object@sT.fun, mw = mw, strictly.positive = FALSE)
}
if(plot.series){
plot(t, result[1,], type = "l", ylab = "Y", ylim = range(result))
for(i in 2:nsim) lines(t, result[i,], col = i)
}
}else{
result <- drawSDE(object@phi, object@gamma2, t, object@b.fun, y0.fun = function(phi, t) y0,
sT = object@sT.fun, mw = mw, strictly.positive = FALSE)
if(plot.series){
plot(t, result, type = "l", ylab = "Y")
}
}

return(result)
})

########
#' Simulation of hierarchical (mixed) diffusion model
#'
#' @description Simulation of the stochastic process model
#'   \eqn{dY_t = b(\phi_j,t,Y_t)dt + \gamma \widetilde{s}(t,Y_t)dW_t, \phi_j~N(\mu, \Omega)}.
#' @param object class object of parameters: "mixedDiffusion"
#' @param nsim number of data sets to simulate. Default is 1.
#' @param seed optional: seed number for random number generator
#' @param t vector of time points
#' @param mw mesh width for finer Euler approximation to simulate time-continuity
#' @param plot.series logical(1), if TRUE, simulated series are depicted grafically
#' @examples
#' mu <- 2; Omega <- 0.4; phi <- matrix(rnorm(21, mu, sqrt(Omega)))
#' model <- set.to.class("mixedDiffusion", y0.fun = function(phi, t) 0.5,
#'   parameter = list(phi = phi, mu = mu, Omega = Omega, gamma2 = 0.1),
#'   b.fun = function(phi, t, x) phi*x, sT.fun = function(t, x) x)
#' t <- seq(0, 1, by = 0.01)
#' data <- simulate(model, t = t, plot.series = TRUE)
#' @export
setMethod(f = "simulate", signature = "mixedDiffusion",
definition = function(object, nsim = 1, seed = NULL, t, mw = 1, plot.series = TRUE) {
set.seed(seed)
if(nsim > 1){
result <- list()
for(j in 1:nsim){
Y <- matrix(0, nrow(object@phi), length(t))
for(i in 1:nrow(object@phi)){
Y[i,] <- drawSDE(object@phi[i,], object@gamma2, t, object@b.fun, y0.fun = object@y0.fun,
sT = object@sT.fun, mw = mw, strictly.positive = FALSE)
}
result[[j]] <- Y
}
if(plot.series){
plot(t, result[[1]][1,], type = "l", ylab = "Y", ylim = range(result[[1]]))
for(i in 2:nrow(object@phi)) lines(t, result[[1]][i,])
}
}else{
result <- matrix(0, nrow(object@phi), length(t))
for(i in 1:nrow(object@phi)){
result[i,] <- drawSDE(object@phi[i,], object@gamma2, t, object@b.fun, y0.fun = object@y0.fun,
sT = object@sT.fun, mw = mw, strictly.positive = FALSE)
}
if(plot.series){
plot(t, result[1,], ylim = range(result), type = "l", ylab="Y")
for(i in 2:nrow(object@phi)) lines(t, result[i,])
}
}

return(result)
})

########
#' Simulation of regression model
#'
#' @description Simulation of the regression model
#'   \eqn{y_i = f(\phi, t_i) + \epsilon_i, \epsilon_i\sim N(0,\gamma^2\widetilde{s}(t_i))}.
#' @param object class object of parameters: "Diffusion"
#' @param nsim number of trajectories to simulate. Default is 1.
#' @param seed optional: seed number for random number generator
#' @param t vector of time points
#' @param plot.series logical(1), if TRUE, simulated series are depicted grafically
#' @examples
#' model <- set.to.class("Regression", parameter = list(phi = 5, gamma2 = 0.1),
#'    fun = function(phi, t) phi*t)
#' t <- seq(0, 1, by = 0.01)
#' data <- simulate(model, t = t, plot.series = TRUE)
#' @export
setMethod(f = "simulate", signature = "Regression",
definition = function(object, nsim = 1, seed = NULL, t, plot.series = TRUE) {
set.seed(seed)

if(nsim > 1){
result <- matrix(0, nsim, length(t))
for(i in 1:nsim){
result[i,] <- object@fun(object@phi, t) + rnorm(length(t), 0, sqrt(object@gamma2*object@sT.fun(t)))
}
if(plot.series){
plot(t, result[1,], ylab = "Y", ylim = range(result))
for(i in 2:nsim) points(t, result[i,], col = i)
}
}else{
result <- object@fun(object@phi, t) + rnorm(length(t), 0, sqrt(object@gamma2*object@sT.fun(t)))
if(plot.series){
plot(t, result, ylab = "Y")
}

}
return(result)
})

########
#' Simulation of hierarchical (mixed) regression model
#'
#' @description Simulation of regression model
#'   \eqn{y_{ij} = f(\phi_j, t_{ij}) + \epsilon_{ij}, \phi_j\sim N(\mu, \Omega),
#'   \epsilon_{ij}\sim N(0,\gamma^2\widetilde{s}(t_{ij}))}.
#' @param object class object of parameters: "mixedRegression"
#' @param nsim number of data sets to simulate. Default is 1.
#' @param seed optional: seed number for random number generator
#' @param t vector of time points
#' @param plot.series logical(1), if TRUE, simulated series are depicted grafically
#' @examples
#' mu <- 2; Omega <- 0.4; phi <- matrix(rnorm(21, mu, sqrt(Omega)))
#' model <- set.to.class("mixedRegression",
#'    parameter = list(phi = phi, mu = mu, Omega = Omega, gamma2 = 0.1),
#'    fun = function(phi, t) phi*t, sT.fun = function(t) t)
#' t <- seq(0, 1, by = 0.01)
#' data <- simulate(model, t = t, plot.series = TRUE)
#' @export
setMethod(f = "simulate", signature = "mixedRegression",
definition = function(object, nsim = 1, seed = NULL, t, plot.series = TRUE) {
set.seed(seed)

if(nsim > 1){
result <- list()
for(j in 1:nsim){
Y <- matrix(0, nrow(object@phi), length(t))
for(i in 1:nrow(object@phi)){
Y[i,] <- object@fun(object@phi[i,], t) + rnorm(length(t), 0, sqrt(object@gamma2*object@sT.fun(t)))
}
result[[j]] <- Y
}
if(plot.series){
plot(t, result[[1]][1,], type = "l", ylab = "Y", ylim = range(result[[1]]))
for(i in 2:nrow(object@phi)) lines(t, result[[1]][i,])
}
}else{
result <- matrix(0, nrow(object@phi), length(t))
for(i in 1:nrow(object@phi)){
result[i,] <- object@fun(object@phi[i,], t) + rnorm(length(t), 0, sqrt(object@gamma2*object@sT.fun(t)))
}
if(plot.series){
plot(t, result[1,], ylim = range(result), type = "l", ylab="Y")
for(i in 2:nrow(object@phi)) lines(t, result[i,])
}
}

return(result)
})

########
#' Simulation of hidden diffusion process
#'
#' @description Simulation of a hidden stochastic process model
#'   \eqn{Z_i = Y_{t_i} + \epsilon_i, dY_t = b(\phi,t,Y_t)dt + \gamma \widetilde{s}(t,Y_t)dW_t,
#'   \epsilon_i\sim N(0,\sigma^2), Y_{t_0}=y_0(\phi, t_0)}.
#' @param object class object of parameters: "hiddenDiffusion"
#' @param nsim number of trajectories to simulate. Default is 1.
#' @param seed optional: seed number for random number generator
#' @param t vector of time points
#' @param mw mesh width for finer Euler approximation to simulate time-continuity
#' @param plot.series logical(1), if TRUE, simulated series are depicted grafically
#' @examples
#' model <- set.to.class("hiddenDiffusion", parameter = list(phi = 0.5, gamma2 = 0.01, sigma2 = 0.1))
#' t <- seq(0, 1, by = 0.01)
#' data <- simulate(model, t = t, plot.series = TRUE)
#' @export
setMethod(f = "simulate", signature = "hiddenDiffusion",
definition = function(object, nsim = 1, seed = NULL, t, mw = 10, plot.series = TRUE) {
set.seed(seed)
if(nsim > 1){
result <- matrix(0, nsim, length(t))
for(i in 1:nsim){
result[i,] <- drawSDE(object@phi, object@gamma2, t, object@b.fun, y0.fun = object@y0.fun,
sT = object@sT.fun, mw = mw, strictly.positive = FALSE)
}
obs <- apply(result, 2, function(vec) rnorm(length(vec), vec, sqrt(object@sigma2)))
if(plot.series){
plot(t, result[1,], type = "l", ylab = "Y", ylim = range(obs))
for(i in 2:nsim) lines(t, result[i,], col = i)
for(i in 1:nsim) points(t, obs[i,], col = i)
}
}else{
result <- drawSDE(object@phi, object@gamma2, t, object@b.fun, y0.fun = object@y0.fun,
sT = object@sT.fun, mw = mw, strictly.positive = FALSE)

obs <- rnorm(length(t), result, sqrt(object@sigma2))
if(plot.series){
plot(t, obs, ylab="simulations"); lines(t, result)
legend("topleft", "hidden process", col = 1, lty = 1, box.lty = 0, inset = 0.01)
}
result <- list(Z = obs, Y = result)

}

return(result)
})

########
#' Simulation of hierarchical (mixed) hidden diffusion model
#'
#' @description Simulation of a stochastic process
#'   \eqn{Z_{ij} = Y_{t_{ij}} + \epsilon_{ij}, dY_t = b(\phi_j,t,Y_t)dt + \gamma \widetilde{s}(t,Y_t)dW_t, \phi_j\sim N(\mu, \Omega),
#'   Y_{t_0}=y_0(\phi, t_0), \epsilon_{ij}\sim N(0,\sigma^2)}.
#' @param object class object of parameters: "hiddenmixedDiffusion"
#' @param nsim number of data sets to simulate. Default is 1.
#' @param seed optional: seed number for random number generator
#' @param t vector of time points
#' @param mw mesh width for finer Euler approximation to simulate time-continuity
#' @param plot.series logical(1), if TRUE, simulated series are depicted grafically
#' @examples
#' mu <- c(5, 1); Omega <- c(0.9, 0.04)
#' phi <- cbind(rnorm(21, mu[1], sqrt(Omega[1])), rnorm(21, mu[2], sqrt(Omega[2])))
#' y0.fun <- function(phi, t) phi[2]
#' model <- set.to.class("hiddenmixedDiffusion", y0.fun = y0.fun,
#'    b.fun = function(phi, t, y) phi[1],
#'    parameter = list(phi = phi, mu = mu, Omega = Omega, gamma2 = 1, sigma2 = 0.01))
#' t <- seq(0, 1, by = 0.01)
#' data <- simulate(model, t = t)
#' @export
setMethod(f = "simulate", signature = "hiddenmixedDiffusion",
definition = function(object, nsim = 1, seed = NULL, t, mw = 10, plot.series = TRUE) {
set.seed(seed)
if(nsim > 1){
result1 <- list()
for(j in 1:nsim){
Y <- matrix(0, nrow(object@phi), length(t))
for(i in 1:nrow(object@phi)){
Y[i,] <- drawSDE(object@phi[i,], object@gamma2, t, object@b.fun, y0.fun = object@y0.fun,
sT = object@sT.fun, mw = mw, strictly.positive = FALSE)
}
result1[[j]] <- Y
}
result2 <- lapply(result1, function(Y) apply(Y, 2, function(vec) rnorm(length(vec), vec, sqrt(object@sigma2))))
if(plot.series){

par(mfrow = c(2,1))
plot(t, result1[[1]][1,], type = "l", ylab = "Y", ylim = range(result2[[1]]))
for(i in 2:nrow(object@phi)) lines(t, result1[[1]][i,])
plot(t, result2[[1]][1,], type = "l", ylab = "Z", ylim = range(result2[[1]]))
for(i in 2:nrow(object@phi)) lines(t, result2[[1]][i,])
par(old.settings)

}
}else{
result1 <- matrix(0, nrow(object@phi), length(t))
for(i in 1:nrow(object@phi)){
result1[i,] <- drawSDE(object@phi[i,], object@gamma2, t, object@b.fun, y0.fun = object@y0.fun,
sT = object@sT.fun, mw = mw, strictly.positive = FALSE)
}

result2 <- apply(result1, 2, function(vec) rnorm(length(vec), vec, sqrt(object@sigma2)))
if(plot.series){

par(mfrow = c(2,1))
plot(t, result1[1,], ylim = range(result2), type = "l", ylab = "Y")
for(i in 2:nrow(object@phi)) lines(t, result1[i,])
plot(t, result2[1,], ylim = range(result2), type = "l", ylab = "Z")
for(i in 2:nrow(object@phi)) lines(t, result2[i,])
par(old.settings)

}
result <- list(Z = result2, Y = result1)

}

return(result)
})

########
#' Simulation of Poisson process
#'
#' @description Simulation of non-homogeneous Poisson process with cumulative intensity function \eqn{\Lambda(t, \xi)}.
#' @param object class object of parameters: "NHPP"
#' @param nsim number of trajectories to simulate. Default is 1.
#' @param seed optional: seed number for random number generator
#' @param t vector of time points
#' @param plot.series logical(1), if TRUE, simulated series are depicted grafically
#' @examples
#' model <- set.to.class("NHPP", parameter = list(xi = c(5, 1/2)),
#'                       Lambda = function(t, xi) (t/xi[2])^xi[1])
#' t <- seq(0, 1, by = 0.01)
#' data <- simulate(model, t = t)
#' @export
setMethod(f = "simulate", signature = "NHPP",
definition = function(object, nsim = 1, seed = NULL, t, plot.series = TRUE) {
set.seed(seed)

result <- simN(t, object@xi, len = nsim, Lambda = object@Lambda)
if(plot.series){
if(nsim > 1){
plot(t, result$N[1,], ylab="N", type="l", ylim = range(result$N))
for(i in 2:nsim) lines(t, result$N[i,]) }else{ plot(t, result$N, ylab="N", type="l")
}
}

return(result)
})

########
#' Simulation of jump diffusion process
#'
#' @description Simulation of jump diffusion process
#'   \eqn{dY_t = b(\phi,t,Y_t)dt + s(\gamma,t,Y_t)dW_t + h(\eta,t,Y_t)dN_t}.
#' @param object class object of parameters: "jumpDiffusion"
#' @param nsim number of trajectories to simulate. Default is 1.
#' @param seed optional: seed number for random number generator
#' @param t vector of time points
#' @param y0 starting point of process
#' @param start vector: start[1] starting point time, start[2] starting point for Poisson process
#' @param mw mesh width for finer Euler approximation to simulate time-continuity
#' @param plot.series logical(1), if TRUE, simulated series are depicted grafically
#' @examples
#' model <- set.to.class("jumpDiffusion",
#'    parameter = list(theta = 0.1, phi = 0.05, gamma2 = 0.1, xi = c(3, 1/4)),
#'    Lambda = function(t, xi) (t/xi[2])^xi[1])
#' t <- seq(0, 1, by = 0.01)
#' data <- simulate(model, t = t, y0 = 0.5)
#' @export
setMethod(f = "simulate", signature = "jumpDiffusion",
definition = function(object, nsim = 1, seed = NULL, t, y0, start = c(0,0), mw = 1, plot.series = TRUE) {
set.seed(seed)
N <- simN(t, object@xi, len = nsim, start = start, Lambda = object@Lambda)$N phi <- object@phi theta <- object@theta gamma2 <- object@gamma2 b.fun <- object@b.fun s.fun <- object@s.fun h.fun <- object@h.fun start <- y0 lt <- length(t) lt2 <- (lt-1)*mw+1 t2 <- seq(min(t), max(t), length = lt2) dt2 <- t2[2] - t2[1] if(nsim > 1){ dN <- t(apply(N, 1, diff)) number <- nrow(N) dN2 <- matrix(0, number, lt2-1) dN2[, seq(1, lt2-1, by = mw)] <- dN Y <- matrix(y0, number, lt2) for(i in 2:lt2){ W <- rnorm(number, 0, sqrt(dt2)) Y[,i] <- Y[,i-1] + b.fun(phi, t2[i-1], Y[,i-1])*dt2 + s.fun(gamma2, t2[i-1], Y[,i-1])*W + h.fun(theta, t2[i-1], Y[,i-1])*dN2[,i-1] } Y <- Y[, seq(1, lt2, by = mw)] } if(nsim == 1){ Y <- rep(y0, lt2) dN2 <- numeric(lt2-1) dN2[seq(1, lt2-1, by = mw)] <- diff(N) for(i in 2:lt2){ W <- rnorm(1, 0, sqrt(dt2)) Y[i] <- Y[i-1] + b.fun(phi, t2[i-1], Y[i-1])*dt2 + s.fun(gamma2, t2[i-1], Y[i-1])*W + h.fun(theta, t2[i-1], Y[i-1])*dN2[i-1] } Y <- Y[seq(1, lt2, by = mw)] } result <- list(N = N, Y = Y) if(plot.series){ old.settings <- par(no.readonly = TRUE) if(nsim > 1){ par(mfrow = c(2,1)) plot(t, N[1,], type = "l", ylab = "N", ylim = range(N)) for(i in 2:nsim) lines(t, N[i,]) plot(t, result$Y[1,], type = "l", ylab = "Y", ylim = range(result$Y)) for(i in 1:nsim) lines(t, result$Y[i,])

}else{

par(mfrow = c(2,1))
plot(t, N, type = "l")
plot(t, result$Y, type = "l", ylab = "Y") } par(old.settings) } return(result) }) ######## #' Simulation of jump diffusion process #' #' @description Simulation of jump diffusion process #' \eqn{Y_t = y_0 \exp( \phi t - \gamma^2/2 t+\gamma W_t + \log(1+\theta) N_t)}. #' @param object class object of parameters: "Merton" #' @param nsim number of trajectories to simulate. Default is 1. #' @param seed optional: seed number for random number generator #' @param t vector of time points #' @param y0 starting point of process #' @param start vector: start[1] starting point time, start[2] starting point for Poisson process #' @param plot.series logical(1), if TRUE, simulated series are depicted grafically #' @examples #' model <- set.to.class("Merton", parameter = list(thetaT = 0.1, phi = 0.05, gamma2 = 0.1, xi = 10)) #' t <- seq(0, 1, by = 0.01) #' data <- simulate(model, t = t, y0 = 0.5) #' @export setMethod(f = "simulate", signature = "Merton", definition = function(object, nsim = 1, seed = NULL, t, y0, start = c(0,0), plot.series = TRUE) { set.seed(seed) N <- simN(t, object@xi, len = nsim, start = start, Lambda = object@Lambda)$N

phi <- object@phi
thetaT <- object@thetaT
gamma2 <- object@gamma2

l <- length(t)
if(nsim > 1){
Y <- t(sapply(1:nsim, function(i){
y0*exp( (phi-gamma2/2)*(t-start[1])+sqrt(gamma2)*cumsum(c(0, rnorm(l-1, 0, sqrt(diff(t[-l])))))+thetaT*(N[i,]-start[2]) )
}) )
}
if(nsim == 1){
Y <- y0*exp( (phi-gamma2/2)*(t-start[1])+sqrt(gamma2)*cumsum(c(0, rnorm(l-1, 0, sqrt(diff(t[-l])))))+thetaT*(N-start[2]) )
}

result <- list(N = N, Y = Y)

if(plot.series){

if(nsim > 1){
par(mfrow = c(2,1))
plot(t, N[1,], type = "l", ylab = "N", ylim = range(N))
for(i in 2:nsim) lines(t, N[i,])
plot(t, result$Y[1,], type = "l", ylab = "Y", ylim = range(result$Y))
for(i in 1:nsim) lines(t, result$Y[i,]) }else{ par(mfrow = c(2,1)) plot(t, N, type = "l") plot(t, result$Y, type = "l", ylab = "Y")
}
par(old.settings)

}
return(result)
})

########
#' Simulation of regression model dependent on Poisson process
#'
#' @description Simulation of of the regression model
#'   \eqn{y_i = f(t_i, N_{t_i}, \theta) + \epsilon_i} with
#'   \eqn{N_t\sim Pois(\Lambda(t, \xi)), \epsilon_i\sim N(0,\gamma^2\widetilde{s}(t))}.
#' @param object class object of parameters: "jumpRegression"
#' @param nsim number of trajectories to simulate. Default is 1.
#' @param seed optional: seed number for random number generator
#' @param t vector of time points
#' @param plot.series logical(1), if TRUE, simulated series are depicted grafically
#' @examples
#' model <- set.to.class("jumpRegression", fun = function(t, N, theta) theta[1]*t + theta[2]*N,
#'    parameter = list(theta = c(1,2), gamma2 = 0.1, xi = 10))
#' t <- seq(0, 1, by = 0.01)
#' data <- simulate(model, t = t)
#' @export
setMethod(f = "simulate", signature = "jumpRegression",
definition = function(object, nsim = 1, seed = NULL, t, plot.series = TRUE) {
set.seed(seed)

N <- simN(t, object@xi, len = nsim, Lambda = object@Lambda)$N if(nsim > 1){ result <- apply(N, 2, function(Nt) object@fun(t, Nt, object@theta) + rnorm(length(t), 0, sqrt(object@gamma2*object@sT.fun(t)))) }else{ result <- object@fun(t, N, object@theta) + rnorm(length(t), 0, sqrt(object@gamma2*object@sT.fun(t))) } result <- list(N = N, Y = result) if(plot.series){ old.settings <- par(no.readonly = TRUE) if(nsim > 1){ par(mfrow = c(2,1)) plot(t, N[1,], type = "l", ylab = "N", ylim = range(N)) for(i in 2:nsim) lines(t, N[i,]) plot(t, result$Y[1,], type = "l", ylab = "Y", ylim = range(result$Y)) for(i in 1:nsim) lines(t, result$Y[i,])
}else{
par(mfrow = c(2,1))
plot(t, N, type = "l")
plot(t, result\$Y, ylab = "Y")
}
par(old.settings)

}
return(result)
})

drawSDE <- function(phi, gamma2, t, b, y0.fun, sT, mw = 10, strictly.positive = TRUE){
lt <- length(t)
lt2 <- (lt-1)*mw+1
t2 <- seq(min(t), max(t), length = lt2)
dt2 <- t2[2] - t2[1]
X <- numeric(lt2)
X[1] <- -1

if(strictly.positive){
while(any(X < 0)){
err <- rnorm(lt2-1, 0, sqrt(dt2))
X[1] <- y0.fun(phi, t2[1])
for(j in 2:lt2){
X[j] <- X[j-1] + b(phi, t2[j-1], X[j-1])*dt2 + sqrt(gamma2)*sT(t2[j-1], X[j-1])*err[j-1]
}
}
} else {
err <- rnorm(lt2-1, 0, sqrt(dt2))
X[1] <- y0.fun(phi, t2[1])
for(j in 2:lt2){
X[j] <- X[j-1] + b(phi, t2[j-1], X[j-1])*dt2 + sqrt(gamma2)*sT(t2[j-1], X[j-1])*err[j-1]
}
}

X[seq(1, lt2, by = mw)]
}

simN <- function(t, xi, len, start = c(0,0), Lambda, int = c("Weibull","Exp")){  # start=(T_n,n)
if(missing(Lambda)){
int <- match.arg(int)
if(int == "Weibull"){
Lambda <- function(t){
(t/xi[2])^xi[1]
}
F_1 <- function(u, Tn_1){
xi[2]*(Lambda(Tn_1)-log(1-u))^(1/xi[1])
}
}else{
F_1 <- function(u, Tn_1){
(log(exp(xi[1]*Tn_1+xi[2])-log(1-u))-xi[2])/xi[1]
}
}
drawTn <- function(Tn_1){
help <- F_1(runif(1), Tn_1)
if(help < t[1]){
return(t[1])
}else{
return(t[which(min(abs(help-t)) == abs(help-t))[1]])
}
}
} else {

drawTn <- function(Tn_1){
prob <- 1-exp(-(Lambda(t, xi)-Lambda(Tn_1, xi)))
u <- runif(1)
t[which(abs(u-prob) == min(abs(u-prob)))]
}
}
lt <- length(t)
drawTn_vec <- function(Tn_1){
sapply(Tn_1, drawTn)
}

Tn <- drawTn_vec(rep(start[1], len))
times <- Tn

while(min(Tn, na.rm = T) < t[lt]){
ind <- which(Tn < t[lt])
Tn[ind] <- drawTn_vec(Tn[ind])
Tn[-ind] <- NA
if(all(is.na(Tn))) break
times <- cbind(times, Tn)
}
if(len > 1){
N_out <- t(apply(times, 1, TimestoN, t)) + start[2]
Times <- times
}else{
N_out <- TimestoN(times, t) + start[2]
Times <- as.numeric(times)
}
list(Times = Times, N = N_out)
}

SimoneHermann/BaPreStoPro documentation built on May 10, 2017, 1:42 p.m.