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

#' Bayesian estimation
#'
#' @description Estimation method for the S4 classes.
#' @param model.class class object with model informations, see \code{\link{set.to.class}}
#' @param t vector or list of time points
#' @param data vector or list or matrix of observation variables
#' @param nMCMC length of Markov chain
#' @param propSd vector of proposal variances
#' @param proposal proposal density: "normal" (default) or "lognormal" (for positive parameters)
#' @param ... parameters dependent on the model class
#' @return class object \code{est.}\code{model.class} containing Markov chains, data input and model informations
#' @references
#' Hermann, S. (2016). BaPreStoPro: an R Package for Bayesian Prediction of Stochastic Processes.
#' SFB 823 discussion paper 28/16.
#'
#' Robert, C. P. and G. Casella (2004). Monte Carlo Statistical Methods. Springer, New York.
#'
#' Rosenthal, J. S. (2011). Optimal Proposal Distributions and Adaptive MCMC. In: Handbook of Markov Chain Monte Carlo, pp. 93-112.
#'
setGeneric("estimate", function(model.class, t, data, nMCMC, propSd, adapt = TRUE, proposal = c("normal", "lognormal"), ...) {
standardGeneric("estimate")
})

########
#' Estimation for diffusion process
#'
#' @description Bayesian estimation of the parameters \eqn{\phi} and \eqn{\gamma^2} of the stochastic process
#'   \eqn{dY_t = b(\phi,t,Y_t)dt + \gamma \widetilde{s}(t,Y_t)dW_t}.
#' @param model.class class of the diffusion process model including all required information, see \code{\link{Diffusion-class}}
#' @param t vector of time points
#' @param data vector of observation variables
#' @param nMCMC length of Markov chain
#' @param propSd vector of proposal variances for \eqn{\phi}
#' @param proposal proposal density: "normal" (default) or "lognormal" (for positive parameters)
#'
#' @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)
#' est_diff <- estimate(model, t, data, 1000)
#' plot(est_diff)
#'
#' @references
#' Hermann, S., K. Ickstadt and C. H. Mueller (2016).
#' Bayesian Prediction of Crack Growth Based on a Hierarchical Diffusion Model.
#' Applied Stochastic Models in Business and Industry, DOI: 10.1002/asmb.2175.
#'
#' @export
setMethod(f = "estimate", signature = "Diffusion",
definition = function(model.class, t, data, nMCMC, propSd, adapt = TRUE, proposal = c("normal", "lognormal")){
proposal <- match.arg(proposal)

if (!is.vector(t, mode = "numeric"))
stop(
"t has to be a vector"
)
if (!is.vector(data, mode = "numeric"))
stop(
"data has to be a vector"
)
if (length(t) != length(data))
stop(
"t and data must have the same length"
)
if(!missing(propSd) && length(model.class@start$phi) != length(propSd)) stop( "propSd must have length of phi" ) if(!is.numeric(nMCMC) || length(nMCMC) > 1 || nMCMC < 1) stop( "nMCMC has to be a natural number" ) if(any(model.class@start$phi < 0) && proposal == "lognormal")
stop(
"lognormal proposal density has positive support and starting value is negative"
)

X <- data
prior <- model.class@prior
start <- model.class@start
bSDE <- model.class@b.fun
sVar <- model.class@sT.fun
len <- nMCMC

priorDensity <- function(phi) dnorm(phi, prior$m.phi, sqrt(prior$v.phi))

if(missing(propSd)) propSd <- abs(prior$m.phi)/5 lt <- length(t) dt <- t[-1] - t[-lt] lphi <- length(propSd) if(proposal == "lognormal"){ proposals <- list() proposals$draw <- function(old, propSd){
proposal(old, propSd)
}
proposals$ratio <- function(drawn, old, propSd){ proposalRatio(old, drawn, propSd) } }else{ proposals <- list() proposals$draw <- function(old, propSd){
rnorm(length(old), old, propSd)
}
proposals$ratio <- function(drawn, old, propSd) 1 } postPhi <- function(lastPhi, gamma2, propSd){ phi_old <- lastPhi phi_drawn <- proposals$draw(phi_old, propSd)
ratio <- prod(priorDensity(phi_drawn) / priorDensity(phi_old))
ratio <- ratio* prod( dnorm(X[-1], X[-lt] + bSDE(phi_drawn, t[-lt], X[-lt])*dt, sqrt(gamma2*sVar(t[-lt], X[-lt])^2*dt))/dnorm(X[-1], X[-lt] + bSDE(phi_old, t[-lt], X[-lt])*dt, sqrt(gamma2*sVar(t[-lt], X[-lt])^2*dt)))
ratio <- ratio* proposals$ratio(phi_drawn, phi_old, propSd) if(is.na(ratio)){ratio <- 0} if(runif(1) < ratio){ phi_old <- phi_drawn } phi_old } postGamma2 <- function(phi){ alphaPost <- prior$alpha.gamma + (lt-1)/2
betaPost <-  prior$beta.gamma + sum( (X[-1] - X[-lt] - bSDE(phi, t[-lt], X[-lt])*dt)^2/(sVar(t[-lt], X[-lt])^2*dt) )/2 1/rgamma(1, alphaPost, betaPost) } phi_out <- matrix(0, len, lphi) gamma2_out <- numeric(len) phi <- start$phi
gamma2 <- start$gamma2 for(count in 1:len){ phi <- postPhi(phi, gamma2, propSd) gamma2 <- postGamma2(phi) phi_out[count, ] <- phi gamma2_out[count] <- gamma2 if (adapt && count%%50 == 0){ propSd <- sapply(1:length(phi), function(i){ ad.propSd(phi_out[(count-50+1):count, i], propSd[i], count/50) }) } } result <- list(phi = phi_out, gamma2 = gamma2_out) if(nMCMC > 100){ he <- matrix(0, ncol(result$phi) + 1, 2)
he[1, ] <- diagnostic(result$gamma2) for(i in 2:(ncol(result$phi)+1)) he[i, ] <- diagnostic(result$phi[,i-1]) burnIn <- max(he[, 1]) thinning <- min( max(he[, 2]), ceiling((nMCMC-burnIn)/100) ) } else { burnIn <- 0; thinning <- 1 } result <- new(Class = "est.Diffusion", phi = result$phi, gamma2 = result$gamma2, model = class.to.list(model.class), t = t, Y = data, burnIn = burnIn, thinning = thinning) return(result) }) ######## #' Estimation for hierarchical (mixed) diffusion model #' #' @description Bayesian estimation of a model #' \eqn{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)}. #' @param model.class class of the hierarchical diffusion model including all required information, see \code{\link{mixedDiffusion-class}} #' @param t list or vector of time points #' @param data list or matrix of observation variables #' @param nMCMC length of Markov chain #' @param propSd vector of proposal variances for \eqn{\phi} #' @param adapt if TRUE (default), proposal variance is adapted #' @param proposal proposal density: "normal" (default) or "lognormal" (for positive parameters) #' @examples #' mu <- 2; Omega <- 0.4; phi <- matrix(rnorm(21, mu, sqrt(Omega))) #' model <- set.to.class("mixedDiffusion", #' 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) #' est <- estimate(model, t, data[1:20,], 100) # nMCMC should be much larger #' plot(est) #' #' # OU #' b.fun <- function(phi, t, y) phi[1]-phi[2]*y; y0.fun <- function(phi, t) phi[3] #' mu <- c(10, 5, 0.5); Omega <- c(0.9, 0.01, 0.01) #' phi <- sapply(1:3, function(i) rnorm(21, mu[i], sqrt(Omega[i]))) #' model <- set.to.class("mixedDiffusion", #' parameter = list(phi = phi, mu = mu, Omega = Omega, gamma2 = 0.1), #' y0.fun = y0.fun, b.fun = b.fun, sT.fun = function(t, x) 1) #' t <- seq(0, 1, by = 0.01) #' data <- simulate(model, t = t, plot.series = TRUE) #' est <- estimate(model, t, data[1:20,], 100) # nMCMC should be much larger #' plot(est) #' #' ## #' t.list <- list() #' for(i in 1:20) t.list[[i]] <- t #' t.list[[21]] <- t[1:50] #' data.list <- list() #' for(i in 1:20) data.list[[i]] <- data[i,] #' data.list[[21]] <- data[21, 1:50] #' est <- estimate(model, t.list, data.list, 100) #' pred <- predict(est, t = t[50:101], which.series = "current", ind.pred = 21, #' b.fun.mat = function(phi, t, y) phi[,1]-phi[,2]*y) #' @references #' Hermann, S., K. Ickstadt and C. H. Mueller (2016). #' Bayesian Prediction of Crack Growth Based on a Hierarchical Diffusion Model. #' Applied Stochastic Models in Business and Industry, DOI: 10.1002/asmb.2175. #' @export setMethod(f = "estimate", signature = "mixedDiffusion", definition = function(model.class, t, data, nMCMC, propSd, adapt = TRUE, proposal = c("normal", "lognormal")) { proposal <- match.arg(proposal) if (!(is.vector(t))) stop( "t has to be a vector or a list" ) if (!(is.matrix(data) || is.list(data))) stop( "data has to be a matrix or a list" ) if (is.vector(t, mode = "numeric") && (length(t) != nrow(data) && length(t) != ncol(data))) stop( "length of t has to be equal to the columns/rows of data" ) if (is.list(t) && ( length(t) != length(data) || !all(sapply(t, length) == sapply(data, length)))) stop( "data must match with t" ) if(!missing(propSd) && ncol(model.class@start$phi) != length(propSd))
stop(
"propSd must have length of phi_j"
)
if(!is.numeric(nMCMC) || length(nMCMC) > 1 || nMCMC < 1)
stop(
"nMCMC has to be a natural number"
)
if(any(model.class@start$phi < 0) && proposal == "lognormal") stop( "lognormal proposal density has positive support and starting value is negative" ) prior <- model.class@prior start <- model.class@start y0.fun <- model.class@y0.fun bSDE <- model.class@b.fun sVar <- model.class@sT.fun len <- nMCMC if(is.matrix(data)){ if(nrow(data) == length(t)){ y <- t(data) } times <- list() y <- list() for(i in 1:nrow(data)){ times[[i]] <- t y[[i]] <- data[i,] } }else{ y <- data times <- t } postOm <- function(phi, mu){ postOmega(prior$alpha.omega, prior$beta.omega, phi, mu) } if(proposal == "lognormal"){ proposals <- list() proposals$draw <- function(old, propSd){
proposal(old, propSd)
}
proposals$ratio <- function(drawn, old, propSd){ proposalRatio(old, drawn, propSd) } }else{ proposals <- list() proposals$draw <- function(old, propSd){
rnorm(length(old), old, propSd)
}
proposals$ratio <- function(drawn, old, propSd) 1 } postPhii_old <- function(lastPhi, mu, Omega, gamma2, X, t, propSd){ # constant y0 lt <- length(t) dt <- diff(t) phi_old <- lastPhi phi_drawn <- proposals$draw(phi_old, propSd)
ratio <- prod(dnorm(phi_drawn, mu, sqrt(Omega)) / dnorm(phi_old, mu, sqrt(Omega)) )
ratio <- ratio* prod( dnorm(X[-1], X[-lt] + bSDE(phi_drawn, t[-lt], X[-lt])*dt, sqrt(gamma2*sVar(t[-lt], X[-lt])^2*dt))/dnorm(X[-1], X[-lt] + bSDE(phi_old, t[-lt], X[-lt])*dt, sqrt(gamma2*sVar(t[-lt], X[-lt])^2*dt)))
ratio <- ratio* proposals$ratio(phi_drawn, phi_old, propSd) if(is.na(ratio)) ratio <- 0 if(runif(1) <= ratio){ phi_old <- phi_drawn } phi_old } postPhii <- function(lastPhi, mu, Omega, gamma2, X, t, propSd){ # X, t vektoren lt <- length(t) dt <- diff(t) phi_old <- lastPhi; lphi <- length(lastPhi) phi_drawn <- phi_old + rnorm(length(mu), 0, propSd) for(k in 1:lphi){ phitest <- phi_drawn; phitest[k] <- rnorm(1, phitest[k], 0.1) if(y0.fun(phi_drawn, t[1]) != y0.fun(phitest, t[1])){ fun <- function(theta){ phi <- phi_drawn; phi[k] <- theta abs(y0.fun(phi, t[1]) - X[1]) } phi_drawn[k] <- optimize(f = fun, phi_drawn[k] + c(-1,1)*start$mu[k], maximum = FALSE)$minimum } } ratio <- prod(dnorm(phi_drawn, mu, sqrt(Omega)) / dnorm(phi_old, mu, sqrt(Omega)) ) ratio <- ratio* prod( dnorm(X[-1], X[-lt] + bSDE(phi_drawn, t[-lt], X[-lt])*dt, sqrt(gamma2*sVar(t[-lt], X[-lt])^2*dt))/dnorm(X[-1], X[-lt] + bSDE(phi_old, t[-lt], X[-lt])*dt, sqrt(gamma2*sVar(t[-lt], X[-lt])^2*dt))) if(is.na(ratio)) ratio <- 0 if(runif(1) <= ratio){ phi_old <- phi_drawn } phi_old } n <- length(y) N.all <- sum(sapply(y, length)-1) postGamma2 <- function(phi){ alphaPost <- prior$alpha.gamma + N.all/2
help <- numeric(n)
for(i in 1:n){
ni <- length(times[[i]])
delta <- diff(times[[i]])
help[i] <- sum( (y[[i]][-1] - y[[i]][-ni] - bSDE(phi[i,], times[[i]][-ni], y[[i]][-ni])*delta)^2/(sVar(times[[i]][-ni], y[[i]][-ni])^2*delta) )
}
betaPost <-  prior$beta.gamma + sum(help)/2 1/rgamma(1, alphaPost, betaPost) } phi_out <- list() mu_out <- matrix(0, len, length(start$mu))
Omega_out <- matrix(0, len, length(start$mu)) gamma2_out <- numeric(len) phi <- start$phi
gamma2 <- start$gamma2 mu <- start$mu
Omega <- postOm(phi, mu)

if(missing(propSd)) propSd <- abs(start$mu)/5 for(count in 1:len){ for(i in 1:n){ phi[i,] <- postPhii(phi[i,], mu, Omega, gamma2, y[[i]], times[[i]], propSd) } mu <- postmu(phi, prior$m.mu, prior$v.mu, Omega) Omega <- postOm(phi, mu) gamma2 <- postGamma2(phi) phi_out[[count]] <- phi mu_out[count,] <- mu Omega_out[count,] <- Omega gamma2_out[count] <- gamma2 if (adapt && count%%50 == 0){ propSd <- sapply(1:length(phi[1,]), function(i){ ad.propSd(sapply(phi_out[(count-50+1):count], function(mat) mat[1,i]), propSd[i], count/50) }) } } result <- list(phi = phi_out, mu = mu_out, Omega = Omega_out, gamma2 = gamma2_out) if(nMCMC > 100){ he <- matrix(0, ncol(result$mu) + 1, 2)
he[1, ] <- diagnostic(result$gamma2) for(i in 2:(ncol(result$mu)+1)) he[i, ] <- diagnostic(result$mu[,i-1]) burnIn <- max(he[, 1]) thinning <- min( max(he[, 2]), ceiling((nMCMC-burnIn)/100) ) } else { burnIn <- 0; thinning <- 1 } if(is.matrix(data)){ result <- new(Class = "est.mixedDiffusion", phi = result$phi, mu = result$mu, Omega = result$Omega, gamma2 = result$gamma2, model = class.to.list(model.class), t = t, Y = data, burnIn = burnIn, thinning = thinning) }else{ result <- new(Class = "est.mixedDiffusion", phi = result$phi, mu = result$mu, Omega = result$Omega, gamma2 = result$gamma2, model = class.to.list(model.class), t = t[[1]], Y = matrix(data[[1]], 1), t.list = t, Y.list = data, burnIn = burnIn, thinning = thinning) } return(result) }) ######## #' Estimation for hidden diffusion process #' #' @description Bayesian estimation of the 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)} with a particle Gibbs sampler. #' @param model.class class of the hidden diffusion model including all required information, see \code{\link{hiddenDiffusion-class}} #' @param t vector of time points #' @param data vector of observation variables #' @param nMCMC length of Markov chain #' @param propSd vector of proposal variances for \eqn{\phi} #' @param adapt if TRUE (default), proposal variance is adapted #' @param proposal proposal density: "normal" (default) or "lognormal" (for positive parameters) #' @param Npart number of particles in the particle Gibbs sampler #' #' @references #' Andrieu, C., A. Doucet and R. Holenstein (2010). Particle Markov Chain Monte Carlo Methods. #' Journal of the Royal Statistical Society B 72, pp. 269-342. #' @examples #' model <- set.to.class("hiddenDiffusion", y0.fun = function(phi, t) 0.5, #' parameter = list(phi = 5, gamma2 = 1, sigma2 = 0.1)) #' t <- seq(0, 1, by = 0.01) #' data <- simulate(model, t = t, plot.series = TRUE) #' est <- estimate(model, t, data$Z, 100)  # nMCMC should be much larger!
#' plot(est)
#'
#' \dontrun{
#' # OU
#' b.fun <- function(phi, t, y) phi[1]-phi[2]*y
#' model <- set.to.class("hiddenDiffusion", y0.fun = function(phi, t) 0.5,
#'                parameter = list(phi = c(10, 1), gamma2 = 1, sigma2 = 0.1),
#'                b.fun = b.fun, sT.fun = function(t, x) 1)
#' t <- seq(0, 1, by = 0.01)
#' data <- simulate(model, t = t, plot.series = TRUE)
#' est <- estimate(model, t, data$Z, 1000) #' plot(est) #' } #' @export setMethod(f = "estimate", signature = "hiddenDiffusion", definition = function(model.class, t, data, nMCMC, propSd, adapt = TRUE, proposal = c("normal", "lognormal"), Npart = 100) { proposal <- match.arg(proposal) if (!is.vector(t, mode = "numeric")) stop( "t has to be a vector" ) if (!is.vector(data, mode = "numeric")) stop( "data has to be a vector" ) if (length(t) != length(data)) stop( "t and data must have the same length" ) if(!missing(propSd) && length(model.class@start$phi) != length(propSd))
stop(
"propSd must have length of phi"
)
if(!is.numeric(nMCMC) || length(nMCMC) > 1 || nMCMC < 1)
stop(
"nMCMC has to be a natural number"
)
if(any(model.class@start$phi < 0) && proposal == "lognormal") stop( "lognormal proposal density has positive support and starting value is negative" ) y <- data prior <- model.class@prior start <- model.class@start len <- nMCMC sigmaTilde <- model.class@sT.fun y0.fun <- model.class@y0.fun b.fun <- model.class@b.fun maxIt <- 1 # input parameter ? if(missing(propSd)) propSd <- abs(start$phi)/5

lt <- length(t)
lphi <- length(start$phi) phi_out <- matrix(0,len,length(start$phi))
sigma2_out <- numeric(len)
gamma2_out <- numeric(len)
X_out <- matrix(0,len,length(t))

phi <- start$phi sigma2 <- start$sigma2
gamma2 <- start$gamma2 postSigma2 <- function(X){ alphaPost <- prior$alpha.sigma + lt/2
betaPost <-  prior$beta.sigma + sum((y-X)^2)/2 1/rgamma(1, alphaPost, betaPost) } postGamma2 <- function(phi, X){ alphaPost <- prior$alpha.gamma + (lt-1)/2
help <- sum( (X[-1] - X[-lt] - b.fun(phi,t[-lt],X[-lt])*(t[-1]-t[-lt]))^2/(sigmaTilde(t[-lt],X[-lt])^2*(t[-1]-t[-lt])) )
betaPost <-  prior$beta.gamma + help/2 1/rgamma(1, alphaPost, betaPost) } priorDensity <- function(phi) dnorm(phi, prior$m.phi, sqrt(prior$v.phi)) if(proposal == "lognormal"){ if(any(start$phi < 0)) warning("Attention: proposal density has positive support")
proposals <- list()
proposals$draw <- function(old, propSd){ proposal(old, propSd) } proposals$ratio <- function(drawn, old, propSd){
proposalRatio(old, drawn, propSd)
}
}else{
proposals <- list()
proposals$draw <- function(old, propSd){ rnorm(length(old), old, propSd) } proposals$ratio <- function(drawn, old, propSd) 1
}

estPhi_and_X <- function(X, lastPhi, gamma2, sigma2, B.fixed, propSd){
likeli <- function(phi,X){
prod(dnorm(X[-1],X[-lt]+b.fun(phi, t[-lt], X[-lt])*diff(t), sqrt(gamma2*diff(t))*sigmaTilde(t[-lt],X[-lt])))
}
# output variables
phi_out <- lastPhi
phi <- lastPhi
for(k in 1:lphi){
for(count in 1:maxIt){
phi[k] <- proposals$draw(phi_out[k], propSd[k]) ratio <- priorDensity(phi)[k]/priorDensity(phi_out)[k] ratio <- ratio * proposals$ratio(phi[k], phi_out[k], propSd[k])
ratio <- ratio*likeli(phi, X)/likeli(phi_out, X)
ratio <- ratio*dnorm(y[1], y0.fun(phi, t[1]), sqrt(sigma2))/dnorm(y[1], y0.fun(phi_out,t[1]), sqrt(sigma2))
if(is.na(ratio)) ratio <- 0

if(runif(1) <= ratio){
phi_out[k] <- phi[k]
break
}
}
}
res_SMC <- SMC(phi_out, gamma2, sigma2, Npart, t, y, b.fun, y0.fun, sigmaTilde, Y.fixed = X, B.fixed = B.fixed)
lign.B <- A.to.B(res_SMC$parents) indice <- sample(1:Npart, 1, prob = res_SMC$W[,lt])
X <- res_SMC$y[indice,] list(phi = phi_out, X = X, B.fixed = lign.B[indice,]) } res_SMC <- SMC(phi, gamma2, sigma2, Npart, t, y, b.fun, y0.fun, sigmaTilde, conditional = FALSE) lign.B <- A.to.B(res_SMC$parents)
indice <- sample(1:Npart, 1, prob =  res_SMC$W[,lt]) X <- res_SMC$y[indice,]
B.fixed <- lign.B[indice,]

for(count in 1:len){

# Particle Gibbs
result <- estPhi_and_X(X, phi, gamma2, sigma2, B.fixed, propSd)
phi <- result$phi X <- result$X
B.fixed <- result$B.fixed gamma2 <- postGamma2(phi, X) sigma2 <- postSigma2(X) phi_out[count,] <- phi sigma2_out[count] <- sigma2 gamma2_out[count] <- gamma2 X_out[count,] <- X if (adapt && count%%50 == 0){ propSd <- sapply(1:lphi, function(i){ ad.propSd(phi_out[(count-50+1):count, i], propSd[i], count/50) }) } if (count%%1000 == 0) message(paste(count, "iterations done")) } result <- list(phi = phi_out, sigma2 = sigma2_out, gamma2 = gamma2_out, X = X_out) if(nMCMC > 100){ he <- matrix(0, ncol(result$phi) + 2, 2)
he[1, ] <- diagnostic(result$gamma2); he[2,] <- diagnostic(result$sigma2)
for(i in 3:(ncol(result$phi)+2)) he[i, ] <- diagnostic(result$phi[,i-2])
burnIn <- max(he[, 1])
thinning <- min( max(he[, 2]), ceiling((nMCMC-burnIn)/100) )
} else {
burnIn <- 0; thinning <- 1
}

result <- new(Class = "est.hiddenDiffusion", phi = result$phi, gamma2 = result$gamma2, sigma2 = result$sigma2, Y.est = result$X,
model = class.to.list(model.class), t = t, Z = data, burnIn = burnIn, thinning = thinning)
return(result)

})

########
#' Estimation for hierarchical (mixed) hidden diffusion process
#'
#' @description Bayesian estimation of the parameters in the hierarchical model:
#'   \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)} with the particle Gibbs sampler.
#' @param model.class class of the hierarchical hidden diffusion model including all required information, see \code{\link{hiddenmixedDiffusion-class}}
#' @param t list or vector of time points
#' @param data list or matrix of observation variables
#' @param nMCMC length of Markov chain
#' @param propSd vector of proposal variances for \eqn{\phi}
#' @param proposal proposal density: "normal" (default) or "lognormal" (for positive parameters)
#' @param Npart number of particles in the particle Gibbs sampler
#' @references
#' Andrieu, C., A. Doucet and R. Holenstein (2010). Particle Markov Chain Monte Carlo Methods.
#' Journal of the Royal Statistical Society B 72, pp. 269-342.
#' @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, plot.series = TRUE)
#'
#' \dontrun{
#' est <- estimate(model, t, data$Z[1:20,], 2000) #' plot(est) #' } #' @export setMethod(f = "estimate", signature = "hiddenmixedDiffusion", definition = function(model.class, t, data, nMCMC, propSd, adapt = TRUE, proposal = c("normal", "lognormal"), Npart = 100) { proposal <- match.arg(proposal) if (!(is.vector(t))) stop( "t has to be a vector or a list" ) if (!(is.matrix(data) || is.list(data))) stop( "data has to be a matrix or a list" ) if (is.vector(t, mode = "numeric") && (length(t) != nrow(data) && length(t) != ncol(data))) stop( "length of t has to be equal to the columns/rows of data" ) if (is.list(t) && ( length(t) != length(data) || !all(sapply(t, length) == sapply(data, length)))) stop( "data must match with t" ) if(!missing(propSd) && ncol(model.class@start$phi) != length(propSd))
stop(
"propSd must have length of phi_j"
)
if(!is.numeric(nMCMC) || length(nMCMC) > 1 || nMCMC < 1)
stop(
"nMCMC has to be a natural number"
)
if(any(model.class@start$phi < 0) && proposal == "lognormal") stop( "lognormal proposal density has positive support and starting value is negative" ) if(is.matrix(data)){ if(nrow(data) == length(t)){ y <- t(data) }else{ if(ncol(data) != length(t)){ stop("length of t has to be equal to the columns of y") } } times <- list() y <- list() for(i in 1:nrow(data)){ times[[i]] <- t y[[i]] <- data[i,] } }else{ y <- data times <- t } prior <- model.class@prior start <- model.class@start len <- nMCMC sigmaTilde <- model.class@sT.fun y0.fun <- model.class@y0.fun b.fun <- model.class@b.fun maxIt <- 1 # input parameter ? if(missing(propSd)) propSd <- abs(start$mu)/5
if(proposal == "lognormal"){
proposals <- list()
proposals$draw <- function(old, propSd){ proposal(old, propSd) } proposals$ratio <- function(drawn, old, propSd){
proposalRatio(old, drawn, propSd)
}
}else{
proposals <- list()
proposals$draw <- function(old, propSd){ rnorm(length(old), old, propSd) } proposals$ratio <- function(drawn, old, propSd) 1
}

lphi <- length(start$mu) n <- length(y) n_all1 <- sum(sapply(y, length) - 1) n_all2 <- sum(sapply(y, length)) postSigma2 <- function(X){ alphaPost <- prior$alpha.sigma + n_all2/2
betaPost <-  prior$beta.sigma + sum((unlist(y)-unlist(X))^2)/2 1/rgamma(1, alphaPost, betaPost) } postGamma2 <- function(X, phi){ # phi matrix, X, t lists, alpha, beta prior parameters alphaPost <- prior$alpha.gamma + n_all1/2

help <- numeric(n)
for(i in 1:n){
ni <- length(times[[i]])
dt <- diff(times[[i]])
help[i] <- sum( (X[[i]][-1] - X[[i]][-ni] - b.fun(phi[i,],times[[i]][-ni],X[[i]][-ni])*dt)^2/(sigmaTilde(times[[i]][-ni],X[[i]][-ni])^2*dt) )
}
betaPost <-  prior$beta.gamma + sum(help)/2 1/rgamma(1, alphaPost, betaPost) } postPhii_Xi <- function(y, t, X, lastPhi, mu, Omega, gamma2, sigma2, B.fixed, propSd){ lt <- length(t) likeli <- function(phi,X){ prod(dnorm(X[-1], X[-lt]+b.fun(phi, t[-lt], X[-lt])*diff(t), sqrt(gamma2*diff(t))*sigmaTilde(t[-lt],X[-lt]))) # true transition density??? } # output variables phi_out <- lastPhi phi <- lastPhi for(k in 1:lphi){ for(count in 1:maxIt){ phi[k] <- proposals$draw(phi_out[k], propSd[k])
ratio <- dnorm(phi[k], mu[k], sqrt(Omega[k]))/dnorm(phi_out[k], mu[k], sqrt(Omega[k]))
ratio <- ratio * proposals$ratio(phi[k], phi_out[k], propSd[k]) ratio <- ratio*likeli(phi, X)/likeli(phi_out, X) ratio <- ratio*dnorm(y[1], y0.fun(phi, t[1]), sqrt(sigma2))/dnorm(y[1], y0.fun(phi_out, t[1]), sqrt(sigma2)) if(is.na(ratio)) ratio <- 0 if(runif(1) <= ratio){ phi_out[k] <- phi[k] break } } } res_SMC <- SMC(phi_out, gamma2, sigma2, Npart, t, y, b.fun, y0.fun, sigmaTilde, Y.fixed = X, B.fixed = B.fixed) lign.B <- A.to.B(res_SMC$parents)
indice <- sample(1:Npart, 1, prob = res_SMC$W[,lt]) X <- res_SMC$y[indice,]

list(phi = phi_out, X = X, B.fixed = lign.B[indice,])

}

n <- length(y)

phi_out <- list()
mu_out <- matrix(0, len, length(start$mu)) Omega_out <- matrix(0, len, length(start$mu))
sigma2_out <- numeric(len)
gamma2_out <- numeric(len)
X_out <- list()

phi <- start$phi sigma2 <- start$sigma2
gamma2 <- start$gamma2 mu <- start$mu
Omega <- postOmega(prior$alpha.omega, prior$beta.omega, phi, mu)
X <- list()
B.fixed <- list()

for(i in 1:n){

if(dnorm(y[[i]][1], y0.fun(phi[i,], times[[i]][1]), sqrt(sigma2)) == 0){
}

result <- SMC(phi[i,], gamma2, sigma2, Npart, times[[i]], y[[i]], b.fun, y0.fun, sigmaTilde, conditional = FALSE)
lign.B <- A.to.B(result$parents) indice <- sample(1:Npart, 1, prob = result$W[,length(times[[i]])])
X[[i]] <- result$y[indice,] B.fixed[[i]] <- lign.B[indice,] } for(count in 1:len){ for(i in 1:n){ help <- postPhii_Xi(y[[i]], times[[i]], X[[i]], phi[i,], mu, Omega, gamma2, sigma2, B.fixed[[i]], propSd) X[[i]] <- help$X
phi[i,] <- help$phi B.fixed[[i]] <- help$B.fixed
}

mu <- postmu(phi, prior$m.mu, prior$v.mu, Omega)
Omega <- postOmega(prior$alpha.omega, prior$beta.omega, phi, mu)

sigma2 <- postSigma2(X)
gamma2 <- postGamma2(X, phi)

phi_out[[count]] <- phi
mu_out[count,] <- mu
Omega_out[count,] <- Omega
sigma2_out[count] <- sigma2
gamma2_out[count] <- gamma2
X_out[[count]] <- X

if (adapt && count%%50 == 0){
propSd <- sapply(1:lphi, function(i){
ad.propSd(sapply(phi_out[(count-50+1):count], function(mat) mat[1, i]), propSd[i], count/50) })
}

if (count%%100 == 0) message(paste(count, "iterations done"))
}
result <- list(phi = phi_out, mu = mu_out, Omega = Omega_out, sigma2 = sigma2_out, gamma2 = gamma2_out, X = X_out)

if(nMCMC > 100){
he <- matrix(0, ncol(result$mu) + 2, 2) he[1, ] <- diagnostic(result$gamma2); he[2,] <- diagnostic(result$sigma2) for(i in 3:(ncol(result$mu)+2)) he[i, ] <- diagnostic(result$mu[,i-2]) burnIn <- max(he[, 1]) thinning <- min( max(he[, 2]), ceiling((nMCMC-burnIn)/100) ) } else { burnIn <- 0; thinning <- 1 } if(is.matrix(data)){ result <- new(Class = "est.hiddenmixedDiffusion", phi = result$phi, mu = result$mu, Omega = result$Omega, gamma2 = result$gamma2, sigma2 = result$sigma2, Y.est = result$X, model = class.to.list(model.class), t = t, Z = data, burnIn = burnIn, thinning = thinning) }else{ result <- new(Class = "est.hiddenmixedDiffusion", phi = result$phi, mu = result$mu, Omega = result$Omega, gamma2 = result$gamma2, sigma2 = result$sigma2, Y.est = result$X, model = class.to.list(model.class), t = t[[1]], Z = matrix(data[[1]], 1), t.list = t, Z.list = data, burnIn = burnIn, thinning = thinning) } return(result) }) ######## #' Estimation for a non-homogeneous Poisson process #' #' @description Bayesian estimation of a non-homogeneous Poisson process (NHPP) with cumulative intensity function \eqn{\Lambda(t, \xi)}. #' @param model.class class of the NHPP model including all required information, see \code{\link{NHPP-class}} #' @param t vector of time points #' @param data vector of observation variables #' @param nMCMC length of Markov chain #' @param propSd vector of proposal variances for \eqn{\xi} #' @param adapt if TRUE (default), proposal variance is adapted #' @param proposal proposal density: "normal" (default) or "lognormal" (for positive parameters) #' @references #' Hermann, S., K. Ickstadt and C. H. Mueller (2015). #' Bayesian Prediction for a Jump Diffusion Process with Application to Crack Growth in Fatigue Experiments. #' SFB 823 discussion paper 30/15. #' @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, plot.series = TRUE) #' est <- estimate(model, t, data$Times, 10000, proposal = "lognormal")
#' plot(est)
#'
#' ##
#' model <- set.to.class("NHPP", parameter = list(xi = 5),
#'                    Lambda = function(t, xi) t*xi)
#' t <- seq(0, 1, by = 0.01)
#' data <- simulate(model, t = t, plot.series = TRUE)
#' est <- estimate(model, t, data$N, 10000) #' plot(est, par.options = list(mfrow = c(1,1))) #' @export setMethod(f = "estimate", signature = "NHPP", definition = function(model.class, t, data, nMCMC, propSd, adapt = TRUE, proposal = c("normal", "lognormal")) { proposal <- match.arg(proposal) if (!is.vector(t, mode = "numeric")) stop( "t has to be a vector" ) if (!is.vector(data, mode = "numeric")) stop( "data has to be a vector" ) if(!missing(propSd) && ncol(model.class@start) != length(propSd)) stop( "propSd must have length of xi" ) if(!is.numeric(nMCMC) || length(nMCMC) > 1 || nMCMC < 1) stop( "nMCMC has to be a natural number" ) if(any(model.class@start < 0) && proposal == "lognormal") stop( "lognormal proposal density has positive support and starting value is negative" ) if(length(t) != length(data)){ jumpTimes <- data }else{ jumpTimes <- dNtoTimes(diff(data), t[-1]) } Tend <- max(t) start <- model.class@start n <- nMCMC Lambda <- model.class@Lambda priorDensity <- model.class@priorDensity lambda <- function(t, xi){ h <- 1e-05 (Lambda(t+h,xi)-Lambda(t,xi))/h } Lik <- function(xi){ lambda_vec <- lambda(jumpTimes, xi) prod(lambda_vec)*exp(-Lambda(Tend, xi)) } if(proposal == "lognormal"){ proposals <- list() proposals$draw <- function(xi_old, propSd){
proposal(xi_old, propSd)
}
proposals$ratio <- function(xi_drawn, xi_old, propSd){ proposalRatio(xi_old, xi_drawn, propSd) } }else{ proposals <- list() proposals$draw <- function(xi_old, propSd){
rnorm(length(xi_old), xi_old, propSd)
}
proposals$ratio <- function(xi_drawn, xi_old, propSd) 1 } if(missing(propSd)) propSd <- (abs(start)+0.1)/2 xi_old <- start xi_out <- matrix(0, length(start), n) LikOld <- Lik(xi_old) for(count in 1:n){ xi_drawn <- proposals$draw(xi_old, propSd)

LikNew <- Lik(xi_drawn)
ratio <- proposals$ratio(xi_drawn, xi_old, propSd) ratio <- ratio*prod(priorDensity(xi_drawn)/priorDensity(xi_old)) ratio <- ratio*LikNew/LikOld if(is.na(ratio)) ratio <- 0 if(runif(1) <= ratio){ xi_old <- xi_drawn LikOld <- LikNew } xi_out[,count] <- xi_old if (adapt && count%%50 == 0){ propSd <- sapply(1:length(start), function(i){ ad.propSd(xi_out[i, (count-50+1):count], propSd[i], count/50) }) } } res <- xi_out if(nMCMC > 100){ if(is.vector(res)){ he <- diagnostic(res); burnIn <- he[1]; thinning <- min( he[2], ceiling((nMCMC-burnIn)/100) ) }else{ he <- matrix(0, ncol(res), 2) for(i in 1:nrow(res)) he[i,] <- diagnostic(res[i,]) burnIn <- max(he[, 1]) thinning <- min( max(he[, 2]), ceiling((nMCMC-burnIn)/100) ) } } else { burnIn <- 0; thinning <- 1 } if(length(t) != length(data)){ result <- new(Class = "est.NHPP", xi = t(res), jumpTimes = as.numeric(data), t = t, N = TimestoN(data, t), model = class.to.list(model.class), burnIn = burnIn, thinning = thinning) } else { result <- new(Class = "est.NHPP", xi = t(res), t = t, N = data, jumpTimes = dNtoTimes(diff(data), t), model = class.to.list(model.class), burnIn = burnIn, thinning = thinning) } return(result) }) ######## #' Estimation for jump diffusion process #' #' @description Bayesian estimation of a stochastic process #' \eqn{dY_t = b(\phi,t,Y_t)dt + s(\gamma^2,t,Y_t)dW_t + h(\theta,t,Y_t)dN_t}. #' @param model.class class of the jump diffusion model including all required information, see \code{\link{jumpDiffusion-class}} #' @param t vector of time points #' @param data vector of observation variables #' @param nMCMC length of Markov chain #' @param propSd vector of proposal variances for \eqn{(\phi, \theta, \gamma^2, \xi)} #' @param adapt if TRUE (default), proposal variance is adapted #' @param proposal proposal density for phi, theta: "normal" (default) or "lognormal" (for positive parameters), see description below #' @param it.xi number of iterations for MH step for \eqn{\xi} inside the Gibbs sampler #' @section Proposal densities: #' For \eqn{\gamma^2}, always the lognormal density is taken, since the parameter is always positive. #' For \eqn{\theta} and \eqn{\phi}, there is the possibility to choose "normal" or "lognormal" (for both together). #' The proposal density for \eqn{\xi} depends on the starting value of \eqn{\xi}. If all components are positive, the proposal density is lognormal, and normal otherwise. #' @examples #' # non-informative #' model <- set.to.class("jumpDiffusion", Lambda = function(t, xi) (t/xi[2])^xi[1], #' parameter = list(theta = 0.1, phi = 0.05, gamma2 = 0.1, xi = c(3, 1/4))) #' t <- seq(0, 1, by = 0.01) #' data <- simulate(model, t = t, y0 = 0.5, plot.series = TRUE) #' est <- estimate(model, t, data, 1000) #' plot(est) #' #' # informative #' model <- set.to.class("jumpDiffusion", Lambda = function(t, xi) (t/xi[2])^xi[1], #' parameter = list(theta = 0.1, phi = 0.05, gamma2 = 0.1, xi = c(3, 1/4)), #' priorDensity = list(phi = function(phi) dnorm(phi, 0.05, 0.01), #' theta = function(theta) dgamma(1/theta, 10, 0.1*9), #' gamma2 = function(gamma2) dgamma(1/gamma2, 10, 0.1*9), #' xi = function(xi) dnorm(xi, c(3, 1/4), c(1,1)))) #' t <- seq(0, 1, by = 0.01) #' data <- simulate(model, t = t, y0 = 0.5, plot.series = TRUE) #' est <- estimate(model, t, data, 1000) #' plot(est) #' #' \dontrun{ #' est_hidden <- estimate(model, t, data$Y, 1000)
#' plot(est_hidden)
#' }
#' @export
setMethod(f = "estimate", signature = "jumpDiffusion",
definition = function(model.class, t, data, nMCMC, propSd, adapt = TRUE, proposal = c("normal", "lognormal"), it.xi = 5) {
proposal <- match.arg(proposal)

if (!is.vector(t, mode = "numeric"))
stop(
"t has to be a vector"
)
if (!is.vector(data))
stop(
"data has to be a vector or a list with entries N and Y"
)
if (is.vector(data, mode = "numeric") && length(t) != length(data))
stop(
"t and data must have the same length"
)
if (is.list(data) &&  !(length(t) == length(data$Y) && length(t) == length(data$N)) )
stop(
"t, data$N and data$Y must have the same length"
)
if(!missing(propSd) && length(model.class@start$xi) + 3 != length(propSd)) stop( "propSd must have length of xi + 3" ) if(!is.numeric(nMCMC) || length(nMCMC) > 1 || nMCMC < 1) stop( "nMCMC has to be a natural number" ) if(any(c(model.class@start$theta, model.class@start$phi) < 0) && proposal == "lognormal") stop( "lognormal proposal density has positive support and starting value is negative" ) if(is.list(data)){ X <- data$Y
N <- data$N jumpTimes <- dNtoTimes(diff(N), t[-1]) }else{ X <- data } Tend <- max(t) n <- nMCMC start <- model.class@start Lambda <- model.class@Lambda lambda <- function(t, xi, h = 1e-05) (Lambda(t+h, xi)-Lambda(t, xi))/h Lik.N <- function(xi, jumpTimes){ lambda_vec <- lambda(jumpTimes, xi) prod(lambda_vec)*exp(-Lambda(Tend, xi)) } b <- model.class@b.fun s <- model.class@s.fun h <- model.class@h.fun priorDensity <- model.class@priorDensity dX <- diff(X) dt <- diff(t) lt <- length(t) # starting values phi <- start$phi
theta <- start$theta gamma2 <- start$gamma2
xi <- start$xi if(proposal == "lognormal"){ proposals <- list() proposals$draw <- function(old, propSd){
proposal(old, propSd)
}
proposals$ratio <- function(drawn, old, propSd){ proposalRatio(old, drawn, propSd) } }else{ proposals <- list() proposals$draw <- function(old, propSd){
rnorm(length(old), old, propSd)
}
proposals$ratio <- function(drawn, old, propSd) 1 } if(all(xi > 0)){ proposals_xi <- list() proposals_xi$draw <- function(xi_old, propSd){
proposal(xi_old, propSd)
}
proposals_xi$ratio <- function(xi_drawn, xi_old, propSd){ proposalRatio(xi_old, xi_drawn, propSd) } } else { proposals_xi <- list() proposals_xi$draw <- function(xi_old, propSd){
rnorm(length(xi_old), xi_old, propSd)
}
proposals_xi$ratio <- function(xi_drawn, xi_old, propSd) 1 } if(is.numeric(data)){ rangeN <- 2 post_dN <- function(dN, i, phi, gamma2, theta, xi){ Di <- dX[i]-b(phi,t[i],X[i])*dt[i]-h(theta,t[i],X[i])*dN dLambda <- Lambda(t[i+1],xi)-Lambda(t[i],xi) exp(-dLambda)/prod(1:max(dN,1))*exp(-Di^2/(2*s(gamma2,t[i],X[i])^2*dt[i])+dN*log(dLambda))/sqrt(2*pi*s(gamma2,t[i],X[i])^2*dt[i]) } drawN <- function(phi, gamma2, theta, xi, dN_old){ dN_new <- dN_old cands <- lapply(dN_old, function(n) 0:(n*rangeN+5)) prob <- lapply(1:(lt-1), function(i) sapply(cands[[i]], post_dN, i, phi, gamma2, theta, xi)) diFu <- lapply(prob, cumsum) ind <- sapply(diFu, function(vec) any(is.na(vec)) | any(is.infinite(vec)) ) u <- numeric(length(diFu)) if(sum(ind) < length(diFu)){ u[!ind] <- runif(length(diFu[!ind]), 0, sapply(diFu[!ind], max)) for(j in (1:(lt-1))[!ind]){ dN_new[j] <- cands[[j]][which(diFu[[j]] >= u[j])[1]] } } dN_new } if(is.null(start$N)){
N <- simN(t, xi, 1, start = c(t[1], 0), Lambda)$N }else{ N <- start$N
}
N_out <- matrix(0, lt, n)
sample.N <- TRUE
}else{
sample.N <- FALSE
}
dN <- diff(N)

likeli <- function(phi, gamma2, theta, dN){
dnorm(dX, b(phi, t[-lt], X[-lt])*dt + h(theta, t[-lt], X[-lt])*dN, s(gamma2, t[-lt], X[-lt])*sqrt(dt))
}

# storage variables
phi_out <- rep(0, n)
theta_out <- rep(0, n)
gamma2_out <- rep(0, n)
xi_out <- matrix(0, n, length(xi))

if(missing(propSd)){
propSd_phi <- start$phi propSd_theta <- start$theta/5
propSd_gamma2 <- start$gamma2/5 propSd_xi <- (abs(start$xi)+0.1)/5
}else{
propSd_phi <- propSd[1]
propSd_theta <- propSd[2]
propSd_gamma2 <- propSd[3]
propSd_xi <- propSd[4:length(propSd)]
}

for(count in 1:n){
if(sample.N){
dN <- drawN(phi, gamma2, theta, xi, dN)
N <- cumsum(c(0, dN))
jumpTimes <- dNtoTimes(dN, t[-1])
}
for(count2 in 1:it.xi){
xi_drawn <- proposals_xi$draw(xi, propSd_xi) ratio <- proposals_xi$ratio(xi_drawn, xi, propSd_xi)
ratio <- ratio*Lik.N(xi_drawn, jumpTimes)/Lik.N(xi, jumpTimes)
ratio <- ratio*prod(priorDensity$xi(xi_drawn)/priorDensity$xi(xi))
if(is.na(ratio)) ratio <- 0

if(runif(1) <= ratio){
xi <- xi_drawn
}
}

phi_drawn <- proposals$draw(phi, propSd_phi) ratio <- prod(likeli(phi_drawn, gamma2, theta, dN)/likeli(phi, gamma2, theta, dN)) ratio <- ratio*priorDensity$phi(phi_drawn)/priorDensity$phi(phi) ratio <- ratio*proposals$ratio(phi_drawn, phi, propSd_phi)
phi[runif(1) <= ratio] <- phi_drawn

theta_drawn <- proposals$draw(theta, propSd_theta) ratio <- prod(likeli(phi, gamma2, theta_drawn, dN)/likeli(phi, gamma2, theta, dN)) ratio <- ratio*priorDensity$theta(theta_drawn)/priorDensity$theta(theta) ratio <- ratio*proposals$ratio(theta_drawn, theta, propSd_theta)
theta[runif(1) <= ratio] <- theta_drawn

gamma2_drawn <- proposal(gamma2, propSd_gamma2)
ratio <- prod(likeli(phi, gamma2_drawn, theta, dN)/likeli(phi, gamma2, theta, dN))
ratio <- ratio*proposalRatio(gamma2, gamma2_drawn, propSd_gamma2)
ratio <- ratio*priorDensity$gamma2(gamma2_drawn)/priorDensity$gamma2(gamma2)
gamma2[runif(1) <= ratio] <- gamma2_drawn

# storage
phi_out[count] <- phi
theta_out[count] <- theta
gamma2_out[count] <- gamma2
xi_out[count, ] <- xi

if(sample.N){
N_out[, count] <- N
if(count %% 1000 == 0) message(paste(count, "iterations are calculated"))
}

if (adapt && count%%50 == 0){
propSd_xi <- sapply(1:length(xi), function(i){
}
}

result <- list(phi = phi_out, gamma2 = gamma2_out, theta = theta_out, xi = xi_out)

if(nMCMC > 100){
he <- matrix(0, 3 + ncol(result$xi), 2) he[1, ] <- diagnostic(result$gamma2); he[2,] <- diagnostic(result$phi); he[3,] <- diagnostic(result$theta)
for(i in 4:(3+ncol(result$xi))) he[i,] <- diagnostic(result$xi[, i-3])
burnIn <- max(he[, 1])
thinning <- min( max(he[, 2]), ceiling((nMCMC-burnIn)/100) )
} else {
burnIn <- 0; thinning <- 1
}

if(is.list(data)){
result <- new(Class = "est.jumpDiffusion", theta = result$theta, phi = result$phi, gamma2 = result$gamma2, xi = result$xi,
model = class.to.list(model.class), t = t, Y = data$Y, N = data$N, burnIn = burnIn, thinning = thinning)
} else {
result <- new(Class = "est.jumpDiffusion", theta = result$theta, phi = result$phi, gamma2 = result$gamma2, xi = result$xi,
N.est = N_out, model = class.to.list(model.class), t = t, Y = data, burnIn = burnIn, thinning = thinning)
}

return(result)
})

########
#' Estimation for jump diffusion process
#'
#' @description Bayesian estimation of a stochastic process
#'   \eqn{Y_t = y_0 \exp( \phi t - \gamma^2/2 t+\gamma W_t + \log(1+\theta) N_t)}.
#' @param model.class class of the jump diffusion model including all required information, see \code{\link{Merton-class}}
#' @param t vector of time points
#' @param data vector of observation variables
#' @param nMCMC length of Markov chain
#' @param propSd vector of proposal variances for \eqn{\xi}
#' @param proposal proposal density for xi: "normal" (default) or "lognormal"
#' @param it.xi number of iterations for MH step for \eqn{\xi} inside the Gibbs sampler
#' @references
#' Hermann, S. and F. Ruggeri (2016). Modelling Wear Degradation in Cylinder Liners.
#' SFB 823 discussion paper 06/16.
#'
#' Hermann, S., K. Ickstadt and C. H. Mueller (2015).
#' Bayesian Prediction for a Jump Diffusion Process with Application to Crack Growth in Fatigue Experiments.
#' SFB 823 discussion paper 30/15.
#' @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, plot.series = TRUE)
#' est <- estimate(model, t, data, 1000)
#' plot(est)
#' \dontrun{
#' est_hidden <- estimate(model, t, data$Y, 1000) #' plot(est_hidden) #' } #' @export setMethod(f = "estimate", signature = "Merton", definition = function(model.class, t, data, nMCMC, propSd, adapt = TRUE, proposal = c("normal", "lognormal"), it.xi = 10) { proposal <- match.arg(proposal) if (!is.vector(t, mode = "numeric")) stop( "t has to be a vector" ) if (!is.vector(data)) stop( "data has to be a vector or a list with entries N and Y" ) if (is.vector(data, mode = "numeric") && length(t) != length(data)) stop( "t and data must have the same length" ) if (is.list(data) && !(length(t) == length(data$Y) && length(t) == length(data$N)) ) stop( "t, data$N and data$Y must have the same length" ) if(!missing(propSd) && length(model.class@start$xi) != length(propSd))
stop(
"propSd must have length of phi"
)
if(!is.numeric(nMCMC) || length(nMCMC) > 1 || nMCMC < 1)
stop(
"nMCMC has to be a natural number"
)

if(is.list(data)){
X <- data$Y N <- data$N
jumpTimes <- dNtoTimes(diff(N), t[-1])
}else{
X <- data
}
Tend <- max(t)

n <- nMCMC
start <- model.class@start
prior <- model.class@prior
Lambda <- model.class@Lambda
lambda <- function(t, xi, h = 1e-05) (Lambda(t+h, xi)-Lambda(t, xi))/h
Lik.N <- function(xi, jumpTimes){
lambda_vec <- lambda(jumpTimes, xi)
prod(lambda_vec)*exp(-Lambda(Tend, xi))
}
priorDensity <- model.class@priorDensity

Delta <- diff(t)
t.l <- t

dlX <- diff(log(X))
x0 <- X[1]
X <- X[-1]
t <- t[-1]
if(is.list(data)){
dN <- diff(N)
N <- N[-1]
}
logX <- log(X)
logx0 <- log(x0)
l <- length(t)
if(any(t == 0)){
help <- matrix(rep(t, l-1), nrow = l-1, ncol = l-1)
Th <- matrix(0, l-1, l-1)
Th[lower.tri(Th, diag = TRUE)] <- t(help)[lower.tri(Th, diag = TRUE)]
Th[upper.tri(Th)] <- help[upper.tri(help)]
T_1 <- solve(Th)
T_1 <- cbind(rep(0, l), rbind(rep(0, l-1), T_1))
}else{
help <- matrix(rep(t, l), nrow = l, ncol = l)
Th <- matrix(0, l, l)
Th[lower.tri(Th, diag = TRUE)] <- t(help)[lower.tri(Th, diag = TRUE)]
Th[upper.tri(Th)] <- help[upper.tri(help)]
T_1 <- solve(Th)
}
postPhi <- function(gamma2, thetaT, N){
Vpost <- 1/( t%*%T_1%*%t/gamma2 + 1/prior$v.phi ) mpost <- Vpost*( prior$m.phi/prior$v.phi + 1/gamma2*t%*%T_1%*%(logX - logx0 + gamma2*t/2 - thetaT*N) ) rnorm(1, mpost, sqrt(Vpost)) } postThetaT <- function(gamma2, phi, N){ Vpost <- 1/( N%*%T_1%*%N/gamma2 + 1/prior$v.thetaT )
mpost <- Vpost*( prior$m.thetaT/prior$v.thetaT + 1/gamma2*N%*%T_1%*%(logX - logx0 - phi*t + gamma2*t/2) )

rnorm(1, mpost, sqrt(Vpost))
}
proposalGamma <- function(gamma2, phi, thetaT, N){
lt <- length(t)
aPost <- prior$alpha.gamma+lt/2 u_X <- logx0 + (phi-gamma2/2)*t + thetaT*N bPost <- prior$beta.gamma+(logX-u_X)%*%T_1%*%(logX-u_X)/2
1/rgamma(1, aPost, bPost)
}
RatioGamma <- function(gamma2, gamma2_drawn, phi, thetaT, N){
h1 <- t%*%T_1%*%t
h2 <- (logX - logx0 - phi*t - thetaT*N)%*%T_1%*%t
exp((1/gamma2_drawn - 1/gamma2)*( (gamma2-gamma2_drawn)*h2/2 + (gamma2-gamma2_drawn)*h1/8))
}

# starting values
phi <- start$phi thetaT <- start$thetaT
gamma2 <- start$gamma2 xi <- start$xi
if(missing(propSd)) propSd <- (abs(start$xi)+0.1)/10 if(proposal == "lognormal"){ proposals <- list() proposals$draw <- function(xi_old, propSd){
proposal(xi_old, propSd)
}
proposals$ratio <- function(xi_drawn, xi_old, propSd){ proposalRatio(xi_old, xi_drawn, propSd) } } else { proposals <- list() proposals$draw <- function(xi_old, propSd){
rnorm(length(xi_old), xi_old, propSd)
}
proposals$ratio <- function(xi_drawn, xi_old, propSd) 1 } if(is.numeric(data)){ rangeN <- 2 post_dN <- function(dN, i, phi, gamma2, thetaT, xi){ Di <- dlX[i]-(phi-gamma2/2)*Delta[i]-thetaT*dN dLambda <- Lambda(t.l[i], xi)-Lambda(t.l[i-1], xi) exp(-dLambda)/prod(1:max(dN,1))*exp(-Di^2/(2*gamma2*Delta[i])+dN*log(dLambda))/sqrt(2*pi*gamma2*Delta[i]) } drawN <- function(phi, gamma2, thetaT, xi, dN_old){ dN_new <- dN_old cands <- lapply(dN_old, function(n) 0:(n*rangeN+5)) prob <- lapply(1:l, function(i) sapply(cands[[i]], post_dN, i, phi, gamma2, thetaT, xi)) diFu <- lapply(prob, cumsum) ind <- sapply(diFu, function(vec) any(is.na(vec)) | any(is.infinite(vec)) ) u <- numeric(length(diFu)) if(sum(ind) < length(diFu)){ u[!ind] <- runif(length(diFu[!ind]), 0, sapply(diFu[!ind], max)) for(j in (1:l)[!ind]){ dN_new[j] <- cands[[j]][which(diFu[[j]]>=u[j])[1]] } } dN_new } if(is.null(start$dN)){
dN <- diff(simN(t.l, xi, 1, start = c(t[1], 0), Lambda = Lambda)$N) }else{ dN <- start$dN
}
N <- cumsum(dN)
sample.N <- TRUE

N_out <- matrix(0, l, n)

}else{
sample.N <- FALSE
}

# storage variables
phi_out <- rep(0, n)
thetaT_out <- rep(0, n)
gamma2_out <- rep(0, n)
xi_out <- matrix(0, n, length(xi))

for(count in 1:n){

if(sample.N){
dN <- drawN(phi, gamma2, thetaT, xi, dN)
N <- cumsum(dN)
jumpTimes <- dNtoTimes(dN, t)
if(count %% 1000 == 0) message(paste(count, "iterations are calculated"))
}
for(count2 in 1:it.xi){
xi_drawn <- proposals$draw(xi, propSd) ratio <- proposals$ratio(xi_drawn, xi, propSd)
ratio <- ratio*Lik.N(xi_drawn, jumpTimes)/Lik.N(xi, jumpTimes)
ratio <- ratio*prod(priorDensity(xi_drawn)/priorDensity(xi))
if(is.na(ratio)) ratio <- 0

if(runif(1) <= ratio){
xi <- xi_drawn
}
}

phi <- postPhi(gamma2, thetaT, N)

thetaT <- postThetaT(gamma2, phi, N)

gamma2_drawn <- proposalGamma(gamma2, phi, thetaT, N)
gamma2[runif(1) <= RatioGamma(gamma2, gamma2_drawn, phi, thetaT, N)] <- gamma2_drawn

# storage
phi_out[count] <- phi
thetaT_out[count] <- thetaT
gamma2_out[count] <- gamma2
xi_out[count, ] <- xi
if(sample.N){
N_out[, count] <- N
}

if (adapt && count%%50 == 0){
propSd <- sapply(1:length(xi), function(i){
}
}

result <- list(phi = phi_out, gamma2 = gamma2_out, thetaT = thetaT_out, xi = xi_out)

if(nMCMC > 100){
he <- matrix(0, 3 + ncol(result$xi), 2) he[1, ] <- diagnostic(result$gamma2); he[2,] <- diagnostic(result$phi); he[3,] <- diagnostic(result$thetaT)
for(i in 4:(3+ncol(result$xi))) he[i,] <- diagnostic(result$xi[, i-3])
burnIn <- max(he[, 1])
thinning <- min( max(he[, 2]), ceiling((nMCMC-burnIn)/100) )
} else {
burnIn <- 0; thinning <- 1
}

if(is.list(data)){
result <- new(Class = "est.Merton", thetaT = result$thetaT, phi = result$phi, gamma2 = result$gamma2, xi = result$xi,
model = class.to.list(model.class), t = t.l, Y = data$Y, N = data$N, burnIn = burnIn, thinning = thinning)
} else {
result <- new(Class = "est.Merton", thetaT = result$thetaT, phi = result$phi, gamma2 = result$gamma2, xi = result$xi, N.est = N_out,
model = class.to.list(model.class), t = t.l, Y = data, burnIn = burnIn, thinning = thinning)
}
return(result)

})

########
#' Estimation for regression model dependent on Poisson process
#'
#' @description Bayesian estimation of the parameter 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 model.class class of the regression model based on the NHPP including all required information, see \code{\link{jumpRegression-class}}
#' @param t vector of time points
#' @param data vector of observation variables
#' @param nMCMC length of Markov chain
#' @param propSd vector of proposal variances for \eqn{(\theta, \xi)}
#' @param proposal proposal density for \eqn{\theta}: "normal" (default) or "lognormal" (for positive parameters)
#' @param it.xi number of iterations for MH step for \eqn{\xi} inside the Gibbs sampler
#' @section Proposal densities:
#' For \eqn{\theta}, there is the possibility to choose "normal" or "lognormal".
#' The proposal density for \eqn{\xi} depends on the starting value of \eqn{\xi}. If all components are positive, the proposal density is lognormal, and normal otherwise.
#'
#' @references
#' Heeke, G., S. Hermann, R. Maurer, K. Ickstadt, and C. H. Mueller (2015).
#' Stochastic Modeling and Statistical Analysis of Fatigue Tests on Prestressed Concrete Beams under Cyclic Loadings.
#' SFB 823 discussion paper 25/15.
#' @examples
#' t <- seq(0,1, by = 0.01)
#' model <- set.to.class("jumpRegression", fun = function(t, N, theta) exp(theta[1]*t) + theta[2]*N,
#'                    parameter = list(theta = c(2, 2), gamma2 = 0.25, xi = c(3, 0.5)),
#'                    Lambda = function(t, xi) (t/xi[2])^xi[1])
#' data <- simulate(model, t = t, plot.series = FALSE)
#' est <- estimate(model, t, data, 1000)
#' plot(est)
#' \dontrun{
#' # work in progress
#' est_hid <- estimate(model, t, data$Y, 1000) #' plot(est_hid) #' } #' @export setMethod(f = "estimate", signature = "jumpRegression", definition = function(model.class, t, data, nMCMC, propSd, adapt = TRUE, proposal = c("normal", "lognormal"), it.xi = 10) { proposal <- match.arg(proposal) if (!is.vector(t, mode = "numeric")) stop( "t has to be a vector" ) if (!is.vector(data)) stop( "data has to be a vector or a list with entries N and Y" ) if (is.vector(data, mode = "numeric") && length(t) != length(data)) stop( "t and data must have the same length" ) if (is.list(data) && !(length(t) == length(data$Y) && length(t) == length(data$N)) ) stop( "t, data$N and data$Y must have the same length" ) if(!missing(propSd) && length(model.class@start$theta) != length(propSd))
stop(
"propSd must have length of phi"
)
if(!is.numeric(nMCMC) || length(nMCMC) > 1 || nMCMC < 1)
stop(
"nMCMC has to be a natural number"
)
if(any(model.class@start$theta < 0) && proposal == "lognormal") stop( "lognormal proposal density has positive support and starting value is negative" ) if(is.list(data)){ Y <- data$Y
N <- data$N jumpTimes <- dNtoTimes(diff(N), t[-1]) }else{ Y <- data } Tend <- max(t) fun <- model.class@fun n <- nMCMC start <- model.class@start prior <- model.class@prior Lambda <- model.class@Lambda lambda <- function(t, xi, h = 1e-05) (Lambda(t+h, xi)-Lambda(t, xi))/h Lik.N <- function(xi, jumpTimes){ lambda_vec <- lambda(jumpTimes, xi) prod(lambda_vec)*exp(-Lambda(Tend, xi)) } sVar <- model.class@sT.fun lt <- length(t) # starting values theta <- start$theta
gamma2 <- start$gamma2 xi <- start$xi

if(proposal == "lognormal"){
proposals <- list()
proposals$draw <- function(old, propSd){ proposal(old, propSd) } proposals$ratio <- function(drawn, old, propSd){
proposalRatio(old, drawn, propSd)
}
}else{
proposals <- list()
proposals$draw <- function(old, propSd){ rnorm(length(old), old, propSd) } proposals$ratio <- function(drawn, old, propSd) 1
}
if(all(xi > 0)){
proposals_xi <- list()
proposals_xi$draw <- function(xi_old, propSd){ proposal(xi_old, propSd) } proposals_xi$ratio <- function(xi_drawn, xi_old, propSd){
proposalRatio(xi_old, xi_drawn, propSd)
}
} else {
proposals <- list()
proposals_xi$draw <- function(xi_old, propSd){ rnorm(length(xi_old), xi_old, propSd) } proposals_xi$ratio <- function(xi_drawn, xi_old, propSd) 1
}

if(is.numeric(data)){

Npart <- 10
rangeN <- 2
A.to.B = function(A){
B <- matrix(0, Npart, lt)
B[,lt] <- 1:Npart;
for (t in seq(lt, 2, by = -1)){
B[,t-1] <- A[B[,t], t-1]}
return(B)
}

CSMC = function(theta, gamma2, xi, N.cond, B.fixed, conditional = TRUE){# conditional SMC
# N.cond = the old samples
x <- matrix(0, Npart, lt)
w <- matrix(1, Npart, lt)
W <- matrix(1, Npart, lt)
parents <-  matrix(1, Npart, lt-1)

# initialisation
x[,1] <- 0  # poisson always starts in 0
if(conditional){
x[B.fixed[1], 1] <- N.cond[1]
}else{
N.cond <- numeric(lt)
}
w[,1] <- dnorm(Y[1], mean = fun(t[1], x[,1], theta), sd = sqrt(gamma2))
W[,1] <- w[,1]/sum(w[,1])

for (n in 2:lt){
if(conditional){
set.parents  <- (1:Npart)[-B.fixed[n]]

On_1 <- rmultinom(1, Npart-1, W[,n-1])
O <- On_1[B.fixed[n-1]] + 1
he <- sample(set.parents, O-1)

parents[B.fixed[n], n-1] <- B.fixed[n-1]
parents[he, n-1] <- B.fixed[n-1]
parents[-c(B.fixed[n],he), n-1] <- sample(set.parents, Npart-O, replace = TRUE, prob =  W[-B.fixed[n], n-1])

}else{
set.parents <- 1:Npart
parents[, n-1] <- sample(1:Npart, Npart, replace = TRUE, prob = W[,n-1])
}

x[,1:(n-1)] <- x[parents[,n-1], 1:(n-1)]
x.past <- x[,n-1]

dr <- function(Ni, dNi_old){
cands <- 0:(dNi_old*rangeN + 5)
prob <- dnorm(Y[n], fun(t[n], Ni+cands, theta), sqrt(2*gamma2))*
dpois(Ni+cands, Lambda(t[n], xi))
diFu <- cumsum(prob)
u <- runif(1, 0, max(diFu))
Ni + cands[which(diFu >= u)[1]]
}
x.new <- sapply(x.past, dr, diff(N.cond)[n-1])

x[set.parents, n] <- x.new[set.parents]
x[-set.parents, n] <- N.cond[n]
wn <- dpois(x[, n], Lambda(t[n], xi))*
dnorm(Y[n], fun(t[n], x[, n], theta), sqrt(gamma2))/
(dnorm(Y[n], fun(t[n], x[, n], theta), sqrt(2*gamma2))*
dpois(x[, n], Lambda(t[n], xi)))
if(sum(wn) == 0) wn <- rep(1, length(wn))
w[, n] <- wn
W[, n] <- w[, n]/sum(w[, n])
}
lign.B <- A.to.B(parents)
indice <- sample(1:Npart, 1, prob = W[, lt])
X <- x[indice, ]
B.fixed <- lign.B[indice,]
return(list(N = X, B.fixed = B.fixed) )
}

if(is.null(start$N)){ # N <- simN(t, xi, 1, start = c(t[1], 0), Lambda = Lambda)$N
he <- CSMC(theta, gamma2, xi, conditional = FALSE)
N <- he$N B.fixed <- he$B.fixed
}else{
N <- start$N } N_out <- matrix(0, lt, n) sample.N <- TRUE }else{ sample.N <- FALSE } ltheta <- length(start$theta)

if(missing(propSd)){
propSd_theta <- abs(prior$m.theta)/20 propSd_xi <- (abs(start$xi)+0.1)/10
}else{
propSd_theta <- propSd[1:ltheta]
propSd_xi <- propSd[(ltheta+1):length(propSd)]
}

postTheta <- function(N, gamma2, lastPhi, propSd){
phi_old <- lastPhi
phi_drawn <- proposals$draw(phi_old, propSd) ratio <- prod(dnorm(phi_drawn, prior$m.theta, sqrt(prior$v.theta)))/prod(dnorm(phi_old, prior$m.theta, sqrt(prior$v.theta))) ratio <- ratio* prod(dnorm(Y, fun(t, N, phi_drawn), sqrt(gamma2*sVar(t)))/dnorm(Y, fun(t, N, phi_old), sqrt(gamma2*sVar(t)))) ratio <- ratio * proposals$ratio(phi_drawn, phi_old, propSd)
if(is.na(ratio)) ratio <- 0
if(runif(1) < ratio){
phi_old <- phi_drawn
}
phi_old
}

postGamma2 <- function(theta, N){
alphaPost <- prior$alpha.gamma + lt/2 betaPost <- prior$beta.gamma + sum((Y-fun(t, N, theta))^2/sVar(t))/2
1/rgamma(1, alphaPost, betaPost)
}

theta_out <- matrix(0, n, ltheta)
gamma2_out <- numeric(n)
xi_out <- matrix(0, n, length(xi))

for(count in 1:n){
if(sample.N){
he <- CSMC(theta, gamma2, xi, N, B.fixed)
N <- he$N jumpTimes <- dNtoTimes(diff(N), t[-1]) B.fixed <- he$B.fixed
}
for(count2 in 1:it.xi){
xi_drawn <- proposals_xi$draw(xi, propSd_xi) ratio <- proposals_xi$ratio(xi_drawn, xi, propSd_xi)
ratio <- ratio*Lik.N(xi_drawn, jumpTimes)/Lik.N(xi, jumpTimes)
if(is.na(ratio)) ratio <- 0

if(runif(1) <= ratio){
xi <- xi_drawn
}
}

theta <- postTheta(N, gamma2, theta, propSd_theta)

gamma2 <- postGamma2(theta, N)

# storage
theta_out[count, ] <- theta
gamma2_out[count] <- gamma2
xi_out[count, ] <- xi
if(sample.N){
N_out[, count] <- N
}
if (adapt && count%%50 == 0){
propSd <- sapply(1:ltheta, function(i){
propSd_xi <- sapply(1:length(xi), function(i){
}

}

result <- list(gamma2 = gamma2_out, theta = theta_out, xi = xi_out)

if(nMCMC > 100){
he <- matrix(0, ncol(result$theta) + ncol(result$xi) + 1, 2)
he[1, ] <- diagnostic(result$gamma2) for(i in 2:(ncol(result$theta)+1)) he[i, ] <- diagnostic(result$theta[,i-1]) for(i in (ncol(result$theta)+2):(ncol(result$theta) + ncol(result$xi) + 1)) he[i, ] <- diagnostic(result$xi[,i - ncol(result$theta) - 1])
burnIn <- max(he[, 1])
thinning <- min( max(he[, 2]), ceiling((nMCMC-burnIn)/100) )
} else {
burnIn <- 0; thinning <- 1
}

if(is.list(data)){
result <- new(Class = "est.jumpRegression", theta = result$theta, gamma2 = result$gamma2, xi = result$xi, model = class.to.list(model.class), t = t, Y = data$Y, N = data$N, burnIn = burnIn, thinning = thinning) }else{ result <- new(Class = "est.jumpRegression", theta = result$theta, gamma2 = result$gamma2, xi = result$xi, N.est = N_out,
model = class.to.list(model.class), t = t, Y = data, burnIn = burnIn, thinning = thinning)
}
return(result)
})

########
#' Estimation for regression model
#'
#' @description Bayesian estimation of the parameter 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 model.class class of the regression model including all required information, see \code{\link{Regression-class}}
#' @param t vector of time points
#' @param data vector of observation variables
#' @param nMCMC length of Markov chain
#' @param propSd vector of proposal variances for \eqn{\phi}
#' @param proposal proposal density: "normal" (default) or "lognormal" (for positive parameters)
#'
#' @references
#' Hermann, S., K. Ickstadt, and C. H. Mueller (2016).
#' Bayesian Prediction of Crack Growth Based on a Hierarchical Diffusion Model.
#' Applied Stochastic Models in Business and Industry, DOI: 10.1002/asmb.2175.
#' @examples
#' t <- seq(0,1, by = 0.01)
#' model <- set.to.class("Regression", fun = function(phi, t) phi[1]*t + phi[2],
#'                    parameter = list(phi = c(1,2), gamma2 = 0.1))
#' data <- simulate(model, t = t, plot.series = TRUE)
#' est <- estimate(model, t, data, 1000)
#' plot(est)
#' @export
setMethod(f = "estimate", signature = "Regression",
definition = function(model.class, t, data, nMCMC, propSd, adapt = TRUE, proposal = c("normal", "lognormal")) {
proposal <- match.arg(proposal)

if (!is.vector(t, mode = "numeric"))
stop(
"t has to be a vector"
)
if (!is.vector(data, mode = "numeric"))
stop(
"data has to be a vector"
)
if (length(t) != length(data))
stop(
"t and data must have the same length"
)
if(!missing(propSd) && length(model.class@start$phi) != length(propSd)) stop( "propSd must have length of phi" ) if(!is.numeric(nMCMC) || length(nMCMC) > 1 || nMCMC < 1) stop( "nMCMC has to be a natural number" ) if(any(model.class@start$phi < 0) && proposal == "lognormal")
stop(
"lognormal proposal density has positive support and starting value is negative"
)

y <- data
prior <- model.class@prior
start <- model.class@start
fODE <- model.class@fun
sVar <- model.class@sT.fun
len <- nMCMC

if(missing(propSd)) propSd <- abs(prior$m.phi)/50 lt <- length(t) lphi <- length(start$phi)

if(proposal == "lognormal"){
proposals <- list()
proposals$draw <- function(old, propSd){ proposal(old, propSd) } proposals$ratio <- function(drawn, old, propSd){
proposalRatio(old, drawn, propSd)
}
}else{
proposals <- list()
proposals$draw <- function(old, propSd){ rnorm(length(old), old, propSd) } proposals$ratio <- function(drawn, old, propSd) 1
}
postPhi <- function(lastPhi, gamma2, propSd){
phi_old <- lastPhi
phi_drawn <- proposals$draw(phi_old, propSd) ratio <- prod(dnorm(phi_drawn, prior$m.phi, sqrt(prior$v.phi)) / dnorm(phi_old, prior$m.phi, sqrt(prior$v.phi))) ratio <- ratio* prod(dnorm(y, fODE(phi_drawn, t), sqrt(gamma2*sVar(t)))/dnorm(y, fODE(phi_old, t), sqrt(gamma2*sVar(t)))) ratio <- ratio * proposals$ratio(phi_drawn, phi_old, propSd)
if(is.na(ratio)) ratio <- 0
if(runif(1) < ratio){
phi_old <- phi_drawn
}
phi_old
}

postGamma2 <- function(phi){
alphaPost <- prior$alpha.gamma + lt/2 betaPost <- prior$beta.gamma + sum((y-fODE(phi, t))^2/sVar(t))/2
1/rgamma(1, alphaPost, betaPost)
}

phi_out <- matrix(0, len, length(prior$m.phi)) gamma2_out <- numeric(len) phi <- start$phi
gamma2 <- start$gamma2 for(count in 1:len){ phi <- postPhi(phi, gamma2, propSd) gamma2 <- postGamma2(phi) phi_out[count,] <- phi gamma2_out[count] <- gamma2 if (adapt && count%%50 == 0){ propSd <- sapply(1:length(phi), function(i){ ad.propSd(phi_out[(count-50+1):count, i], propSd[i], count/50) }) } } result <- list(phi = phi_out, gamma2 = gamma2_out) if(nMCMC > 100){ he <- matrix(0, ncol(result$phi) + 1, 2)
he[1, ] <- diagnostic(result$gamma2) for(i in 2:(ncol(result$phi)+1)) he[i, ] <- diagnostic(result$phi[,i-1]) burnIn <- max(he[, 1]) thinning <- min( max(he[, 2]), ceiling((nMCMC-burnIn)/100) ) } else { burnIn <- 0; thinning <- 1 } result <- new(Class = "est.Regression", phi = result$phi, gamma2 = result$gamma2, model = class.to.list(model.class), t = t, Y = data, burnIn = burnIn, thinning = thinning) return(result) }) ######## #' Estimation for the hierarchical (mixed) regression model #' #' @description Bayesian estimation of the parameter of the hierarchical 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 model.class class of the hierarchical regression model including all required information, see \code{\link{mixedRegression-class}} #' @param t list or vector of time points #' @param data list or matrix of observation variables #' @param nMCMC length of Markov chain #' @param propSd vector of proposal variances for \eqn{\phi} #' @param adapt if TRUE (default), proposal variance is adapted #' @param proposal proposal density: "normal" (default) or "lognormal" (for positive parameters) #' @references #' Hermann, S., K. Ickstadt, and C. H. Mueller (2016). #' Bayesian Prediction of Crack Growth Based on a Hierarchical Diffusion Model. #' Applied Stochastic Models in Business and Industry, DOI: 10.1002/asmb.2175. #' @examples #' mu <- c(10, 5); Omega <- c(0.9, 0.01) #' phi <- cbind(rnorm(21, mu[1], sqrt(Omega[1])), rnorm(21, mu[2], sqrt(Omega[2]))) #' model <- set.to.class("mixedRegression", #' parameter = list(phi = phi, mu = mu, Omega = Omega, gamma2 = 0.1), #' fun = function(phi, t) phi[1]*t + phi[2], sT.fun = function(t) 1) #' t <- seq(0, 1, by = 0.01) #' data <- simulate(model, t = t, plot.series = FALSE) #' est <- estimate(model, t, data[1:20,], 1000) #' plot(est) #' #' @export setMethod(f = "estimate", signature = "mixedRegression", definition = function(model.class, t, data, nMCMC, propSd, adapt = TRUE, proposal = c("normal", "lognormal")) { proposal <- match.arg(proposal) if (!(is.vector(t))) stop( "t has to be a vector or a list" ) if (!(is.matrix(data) || is.list(data))) stop( "data has to be a matrix or a list" ) if (is.vector(t, mode = "numeric") && (length(t) != nrow(data) && length(t) != ncol(data))) stop( "length of t has to be equal to the columns/rows of data" ) if (is.list(t) && ( length(t) != length(data) || !all(sapply(t, length) == sapply(data, length)))) stop( "data must match with t" ) if(!missing(propSd) && ncol(model.class@start$phi) != length(propSd))
stop(
"propSd must have length of phi_j"
)
if(!is.numeric(nMCMC) || length(nMCMC) > 1 || nMCMC < 1)
stop(
"nMCMC has to be a natural number"
)
if(any(model.class@start$phi < 0) && proposal == "lognormal") stop( "lognormal proposal density has positive support and starting value is negative" ) prior <- model.class@prior start <- model.class@start fODE <- model.class@fun sVar <- model.class@sT.fun len <- nMCMC if(missing(propSd)) propSd <- abs(start$mu)/5

if(is.matrix(data)){
if(nrow(data) == length(t)){
y <- t(data)
}else{
if(ncol(data) != length(t)){
stop("length of t has to be equal to the columns of y")
}
}
times <- list()
y <- list()
for(i in 1:nrow(data)){
times[[i]] <- t
y[[i]] <- data[i,]
}
}else{
y <- data
times <- t
}

postOm <- function(phi,mu){
postOmega(prior$alpha.omega, prior$beta.omega, phi, mu)
}
if(proposal == "lognormal"){
proposals <- list()
proposals$draw <- function(old, propSd){ proposal(old, propSd) } proposals$ratio <- function(drawn, old, propSd){
proposalRatio(old, drawn, propSd)
}
}else{
proposals <- list()
proposals$draw <- function(old, propSd){ rnorm(length(old), old, propSd) } proposals$ratio <- function(drawn, old, propSd) 1
}
postPhii <- function(lastPhi, mu, Omega, gamma2, y, t, propSd){
lt <- length(t)
phi_old <- lastPhi

phi_drawn <- proposals$draw(phi_old, propSd) ratio <- prod(dnorm(phi_drawn, mu, sqrt(Omega))/dnorm(phi_old, mu, sqrt(Omega))) ratio <- ratio* prod(dnorm(y, fODE(phi_drawn,t), sqrt(gamma2*sVar(t)))/dnorm(y, fODE(phi_old,t), sqrt(gamma2*sVar(t)))) ratio <- ratio * proposals$ratio(phi_drawn, phi_old, propSd)
if(is.na(ratio)){ratio <- 0}
if(runif(1) < ratio){
phi_old <- phi_drawn
}
phi_old
}
N_all <- sum(sapply(y,length))
n <- length(y)
postGamma2 <- function(phi){
alphaPost <- prior$alpha.gamma + N_all/2 help <- numeric(n) for(i in 1:n){ help[i] <- sum((y[[i]]-fODE(phi[i,], times[[i]]))^2/sVar(times[[i]])) } betaPost <- prior$beta.gamma + sum(help)/2
1/rgamma(1, alphaPost, betaPost)
}

phi_out <- list()
mu_out <- matrix(0, len, length(start$mu)) Omega_out <- matrix(0, len, length(start$mu))
gamma2_out <- numeric(len)

phi <- start$phi gamma2 <- start$gamma2
mu <- start$mu Omega <- postOm(phi, mu) for(count in 1:len){ for(i in 1:n){ phi[i,] <- postPhii(phi[i,], mu, Omega, gamma2, y[[i]], times[[i]], propSd) } mu <- postmu(phi, prior$m.mu, prior$v.mu, Omega) Omega <- postOm(phi, mu) gamma2 <- postGamma2(phi) phi_out[[count]] <- phi mu_out[count,] <- mu Omega_out[count,] <- Omega gamma2_out[count] <- gamma2 if (adapt && count%%50 == 0){ propSd <- sapply(1:length(phi[1,]), function(i){ ad.propSd(sapply(phi_out[(count-50+1):count], function(mat) mat[1,i]), propSd[i], count/50) }) } } result <- list(phi = phi_out, mu = mu_out, Omega = Omega_out, gamma2 = gamma2_out) if(nMCMC > 100){ he <- matrix(0, ncol(result$mu) + 1, 2)
he[1, ] <- diagnostic(result$gamma2) for(i in 2:(ncol(result$mu)+1)) he[i, ] <- diagnostic(result$mu[,i-1]) burnIn <- max(he[, 1]) thinning <- min( max(he[, 2]), ceiling((nMCMC-burnIn)/100) ) } else { burnIn <- 0; thinning <- 1 } if(is.matrix(data)){ result <- new(Class = "est.mixedRegression", phi = result$phi, mu = result$mu, Omega = result$Omega, gamma2 = result$gamma2, model = class.to.list(model.class), t = t, Y = data, burnIn = burnIn, thinning = thinning) }else{ result <- new(Class = "est.mixedRegression", phi = result$phi, mu = result$mu, Omega = result$Omega, gamma2 = result\$gamma2,
model = class.to.list(model.class), t = t[[1]], Y = matrix(data[[1]], 1),
t.list = t, Y.list = data, burnIn = burnIn, thinning = thinning)
}

return(result)

})

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