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

#' Plot method for the Bayesian estimation results
#'
#' @description Plot method for the estimation results of the jump diffusion model.
#' @param x est.jumpDiffusion class, created with method \code{\link{estimate,jumpDiffusion-method}}
#' @param par.options list of options for function par()
#' @param style one out of "chains", "acf", "density"
#' @param par2plot logical vector, which parameters to be plotted, order: \eqn{(\phi, \theta, \gamma^2, \xi, N)}
#' @param reduced logical (1), if TRUE, the chains are thinned and burn-in phase is dropped
#' @param thinning thinning rate, if missing, the proposed one by the estimation procedure is taken
#' @param burnIn burn-in phase, if missing, the proposed one by the estimation procedure is taken
#' @param priorMeans logical(1), if TRUE (default), prior means are marked with a line
#' @param col.priorMean color of the prior mean line, default 2
#' @param lty.priorMean linetype of the prior mean line, default 1
#' @param ... optional plot parameters
#' @examples
#' 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)))
#' data <- simulate(model, t = seq(0, 1, by = 0.01), y0 = 0.5, plot.series = TRUE)
#' est <- estimate(model, t = seq(0, 1, by = 0.01), data, 1000)  # nMCMC small for example
#' plot(est)
#' plot(est, burnIn = 100, thinning = 2, reduced = TRUE)
#' plot(est, par.options = list(mar = c(5, 4.5, 4, 2) + 0.1, mfrow = c(2, 3)), xlab = "iteration")
#' # plot only for phi and xi ...
#' plot(est, style = "acf", main = "", par2plot = c(TRUE, FALSE, FALSE, TRUE, TRUE))
#' plot(est, style = "density", lwd = 2, priorMean = FALSE)
#' plot(est, style = "density", col.priorMean = 1, lty.priorMean = 2, main = "posterior")
#' plot(est, style = "acf", par.options = list(), par2plot = c(TRUE, rep(FALSE, 4)), main = "")
#' @export
#' @import stats
#' @import graphics
#' @import methods
#'
setMethod(f = "plot", signature = "est.jumpDiffusion",
definition = function(x, par.options, style = c("chains", "acf", "density"), par2plot, reduced = FALSE,
thinning, burnIn, priorMeans = TRUE, col.priorMean = 2, lty.priorMean = 1, ...) {
style <- match.arg(style)
if(reduced){
if(missing(thinning)) thinning <- x@thinning
if(missing(burnIn)) burnIn <- x@burnIn
ind <- seq(burnIn + 1, length(x@phi), by = thinning)
}else{
ind <- 1:length(x@phi)
}

p <- ncol(x@xi)
if(style == "chains"){
p1 <- c(1, 0)[c(sum(dim(x@N.est)) > 0, sum(dim(x@N.est)) == 0)]
} else {
p1 <- 0
}
if(missing(par2plot)) par2plot <- rep(TRUE, 3+p+p1)
if(missing(par.options)){
if(sum(par2plot) == 1) par.options <- list(mfrow = c(1, 1))
if(sum(par2plot) > 1) par.options <- list(mfrow = c(ceiling(sum(par2plot)/2), 2))
}

par(par.options)

if(style == "chains"){

if(par2plot[1]){
plot(x@phi[ind], type = "l", ylab = expression(phi), ...)
if(priorMeans) abline(h = x@model$phi, col = col.priorMean, lty = lty.priorMean) } if(par2plot[2]){ plot(x@theta[ind], type = "l", ylab = expression(theta), ...) if(priorMeans) abline(h = x@model$theta, col = col.priorMean, lty = lty.priorMean)
}
if(par2plot[3]){
plot(x@gamma2[ind], type = "l", ylab = expression(gamma^2), ...)
if(priorMeans)  abline(h = x@model$gamma2, col = col.priorMean, lty = lty.priorMean) } if(any(par2plot[3 + seq_len(p)])){ for(i in 1:p){ if(par2plot[3 + i]){ plot(x@xi[ind, i], type = "l", ylab = bquote(xi[.(i)]), ...) if(priorMeans) abline(h = x@model$xi[i], col = col.priorMean, lty = lty.priorMean)
}
}
}

if(sum(dim(x@N.est) > 0)){
if(par2plot[3 + p + 1]){
plot(x@t, x@N.est[,length(x@phi)], type = "l", xlab = "t", ylab = "N", ylim = range(x@N.est[,ind]), ...)
for(i in ind) lines(x@t, x@N.est[,i])
}
}
}
if(style == "acf"){

if(par2plot[1]) acf(x@phi[ind], xlab = expression(phi), ...)
if(par2plot[2]) acf(x@theta[ind], xlab = expression(theta), ...)
if(par2plot[3]) acf(x@gamma2[ind], xlab = expression(gamma^2), ...)
for(i in 1:p){
if(par2plot[3 + i]) acf(x@xi[ind, i], xlab = bquote(xi[.(i)]), ...)
}
}
if(style == "density"){

if(par2plot[1]){
plot(density(x@phi[ind]), xlab = expression(phi), ...)
if(priorMeans) abline(v = x@model$phi, col = col.priorMean, lty = lty.priorMean) } if(par2plot[2]){ plot(density(x@theta[ind]), xlab = expression(theta), ...) if(priorMeans) abline(v = x@model$theta, col = col.priorMean, lty = lty.priorMean)
}
if(par2plot[3]){
plot(density(x@gamma2[ind]), xlab = expression(gamma^2), ...)
if(priorMeans)  abline(v = x@model$gamma2, col = col.priorMean, lty = lty.priorMean) } for(i in 1:p){ if(par2plot[3 + i]){ plot(density(x@xi[ind, i]), xlab = bquote(xi[.(i)]), ...) if(priorMeans) abline(v = x@model$xi[i], col = col.priorMean, lty = lty.priorMean)
}
}
}

par(old.settings[which(names(old.settings) %in% names(par.options))])
})

#' Plot method for the Bayesian estimation results
#'
#' @description Plot method for the estimation results of the Merton model.
#' @param x est.Merton class, created with method \code{\link{estimate,Merton-method}}
#' @param par.options list of options for function par()
#' @param style one out of "chains", "acf", "density"
#' @param par2plot logical vector, which parameters to be plotted, order: \eqn{(\phi, \widetilde{\theta}, \gamma^2, \xi, N)}
#' @param reduced logical (1), if TRUE, the chains are thinned and burn-in phase is dropped
#' @param thinning thinning rate, if missing, the proposed one by the estimation procedure is taken
#' @param burnIn burn-in phase, if missing, the proposed one by the estimation procedure is taken
#' @param priorMeans logical(1), if TRUE (default), prior means are marked with a line
#' @param col.priorMean color of the prior mean line, default 2
#' @param lty.priorMean linetype of the prior mean line, default 1
#' @param ... optional plot parameters
#' @examples
#' model <- set.to.class("Merton", Lambda = function(t, xi) (t/xi[2])^xi[1],
#' parameter = list(thetaT = 0.1, phi = 0.05, gamma2 = 0.1, xi = c(3, 1/4)))
#' data <- simulate(model, t = seq(0, 1, by = 0.01), y0 = 0.5, plot.series = TRUE)
#' est <- estimate(model, t = seq(0, 1, by = 0.01), data, 1000)  # nMCMC small for example
#' plot(est)
#' plot(est, burnIn = 100, thinning = 2, reduced = TRUE)
#' plot(est, par.options = list(mar = c(5, 4.5, 4, 2) + 0.1, mfrow = c(2, 3)), xlab = "iteration")
#' # plot only for phi and xi ...
#' plot(est, style = "acf", main = "", par2plot = c(TRUE, FALSE, FALSE, TRUE, TRUE))
#' plot(est, style = "density", lwd = 2, priorMean = FALSE)
#' plot(est, style = "density", col.priorMean = 1, lty.priorMean = 2, main = "posterior")
#' plot(est, style = "acf", par.options = list(), par2plot = c(TRUE, rep(FALSE, 4)), main = "")
#' @export
setMethod(f = "plot", signature = "est.Merton",
definition = function(x, par.options, style = c("chains", "acf", "density"), par2plot, reduced = FALSE,
thinning, burnIn, priorMeans = TRUE, col.priorMean = 2, lty.priorMean = 1, ...) {
style <- match.arg(style)
if(reduced){
if(missing(thinning)) thinning <- x@thinning
if(missing(burnIn)) burnIn <- x@burnIn
ind <- seq(burnIn + 1, length(x@phi), by = thinning)
}else{
ind <- 1:length(x@phi)
}

p <- ncol(x@xi)
if(style == "chains"){
p1 <- c(1, 0)[c(sum(dim(x@N.est)) > 0, sum(dim(x@N.est)) == 0)]
} else {
p1 <- 0
}
if(missing(par2plot)) par2plot <- rep(TRUE, 3+p+p1)
if(missing(par.options)){
if(sum(par2plot) == 1) par.options <- list(mfrow = c(1, 1))
if(sum(par2plot) > 1) par.options <- list(mfrow = c(ceiling(sum(par2plot)/2), 2))
}

par(par.options)

if(style == "chains"){

if(par2plot[1]){
plot(x@phi[ind], type = "l", ylab = expression(phi), ...)
if(priorMeans) abline(h = x@model$phi, col = col.priorMean, lty = lty.priorMean) } if(par2plot[2]){ plot(x@thetaT[ind], type = "l", ylab = expression(tilde(theta)), ...) if(priorMeans) abline(h = x@model$thetaT, col = col.priorMean, lty = lty.priorMean)
}
if(par2plot[3]){
plot(x@gamma2[ind], type = "l", ylab = expression(gamma^2), ...)
if(priorMeans)  abline(h = x@model$gamma2, col = col.priorMean, lty = lty.priorMean) } if(any(par2plot[3 + seq_len(p)])){ for(i in 1:p){ if(par2plot[3 + i]){ plot(x@xi[ind, i], type = "l", ylab = bquote(xi[.(i)]), ...) if(priorMeans) abline(h = x@model$xi[i], col = col.priorMean, lty = lty.priorMean)
}
}
}

if(sum(dim(x@N.est) > 0)){
if(par2plot[3 + p + 1]){
plot(x@t, x@N.est[,length(x@phi)], type = "l", xlab = "t", ylab = "N", ylim = range(x@N.est[,ind]), ...)
for(i in ind) lines(x@t, x@N.est[,i])
}
}
}
if(style == "acf"){

if(par2plot[1]) acf(x@phi[ind], xlab = expression(phi), ...)
if(par2plot[2]) acf(x@thetaT[ind], xlab = expression(tilde(theta)), ...)
if(par2plot[3]) acf(x@gamma2[ind], xlab = expression(gamma^2), ...)
for(i in 1:p){
if(par2plot[3 + i]) acf(x@xi[ind, i], xlab = bquote(xi[.(i)]), ...)
}
}
if(style == "density"){

if(par2plot[1]){
plot(density(x@phi[ind]), xlab = expression(phi), ...)
if(priorMeans) abline(v = x@model$phi, col = col.priorMean, lty = lty.priorMean) } if(par2plot[2]){ plot(density(x@thetaT[ind]), xlab = expression(tilde(theta)), ...) if(priorMeans) abline(v = x@model$theta, col = col.priorMean, lty = lty.priorMean)
}
if(par2plot[3]){
plot(density(x@gamma2[ind]), xlab = expression(gamma^2), ...)
if(priorMeans)  abline(v = x@model$gamma2, col = col.priorMean, lty = lty.priorMean) } for(i in 1:p){ if(par2plot[3 + i]){ plot(density(x@xi[ind, i]), xlab = bquote(xi[.(i)]), ...) if(priorMeans) abline(v = x@model$xi[i], col = col.priorMean, lty = lty.priorMean)
}
}
}

par(old.settings[which(names(old.settings) %in% names(par.options))])
})

#' Plot method for the Bayesian estimation results
#'
#' @description Plot method for the estimation results of the diffusion model.
#' @param x est.Diffusion class, created with method \code{\link{estimate,Diffusion-method}}
#' @param par.options list of options for function par()
#' @param style one out of "chains", "acf", "density"
#' @param par2plot logical vector, which parameters to be plotted, order: \eqn{(\phi, \gamma^2)}
#' @param reduced logical (1), if TRUE, the chains are thinned and burn-in phase is dropped
#' @param thinning thinning rate, if missing, the proposed one by the estimation procedure is taken
#' @param burnIn burn-in phase, if missing, the proposed one by the estimation procedure is taken
#' @param priorMeans logical(1), if TRUE (default), prior means are marked with a line
#' @param col.priorMean color of the prior mean line, default 2
#' @param lty.priorMean linetype of the prior mean line, default 1
#' @param ... optional plot parameters
#' @examples
#' model <- set.to.class("Diffusion", b.fun = function(phi, t, y) phi[1]-phi[2]*y,
#'     parameter = list(phi = c(10, 1), gamma2 = 0.1))
#' data <- simulate(model, t = seq(0, 1, by = 0.01), y0 = 0.5, plot.series = TRUE)
#' est <- estimate(model, t = seq(0, 1, by = 0.01), data, 1000)  # nMCMC small for example
#' plot(est)
#' plot(est, burnIn = 100, thinning = 2, reduced = TRUE)
#' plot(est, par.options = list(mar = c(5, 4.5, 4, 2) + 0.1, mfrow = c(3,1)), xlab = "iteration")
#' plot(est, style = "acf", main = "", par2plot = c(TRUE, TRUE, FALSE))
#' plot(est, style = "density", lwd = 2, priorMean = FALSE)
#' plot(est, style = "density", col.priorMean = 1, lty.priorMean = 2, main = "posterior")
#' plot(est, style = "acf", par.options = list(), main = "", par2plot = c(FALSE, FALSE, TRUE))
#' @export
setMethod(f = "plot", signature = "est.Diffusion",
definition = function(x, par.options, style = c("chains", "acf", "density"), par2plot, reduced = FALSE,
thinning, burnIn, priorMeans = TRUE, col.priorMean = 2, lty.priorMean = 1,...) {

style <- match.arg(style)
if(reduced){
if(missing(thinning)) thinning <- x@thinning
if(missing(burnIn)) burnIn <- x@burnIn
ind <- seq(burnIn + 1, length(x@gamma2), by = thinning)
}else{
ind <- 1:length(x@gamma2)
}

p <- ncol(x@phi)

if(missing(par2plot)) par2plot <- rep(TRUE, 1+p)
if(missing(par.options)){
if(sum(par2plot) == 1) par.options <- list(mfrow = c(1, 1))
if(sum(par2plot) > 1) par.options <- list(mfrow = c(ceiling(sum(par2plot)/2), 2))
}

par(par.options)

if(style == "chains"){
for(i in 1:p) {
if(par2plot[i]){
plot(x@phi[ind, i], type = "l", ylab = bquote(phi[.(i)]), ...)
if(priorMeans) abline(h = x@model$phi[i], col = col.priorMean, lty = lty.priorMean) } } if(par2plot[p+1]){ plot(x@gamma2[ind], type = "l", ylab = expression(gamma^2), ...) if(priorMeans) abline(h = x@model$gamma2, col = col.priorMean, lty = lty.priorMean)
}
}
if(style == "acf"){
for(i in 1:p) {
if(par2plot[i]){
acf(x@phi[ind, i], xlab = bquote(phi[.(i)]), ...)
}
}
if(par2plot[p+1]) acf(x@gamma2[ind], xlab = expression(gamma^2), ...)
}
if(style == "density"){
for(i in 1:p) {
if(par2plot[i]){
plot(density(x@phi[ind, i]), xlab = bquote(phi[.(i)]), ...)
if(priorMeans) abline(v = x@model$phi[i], col = col.priorMean, lty = lty.priorMean) } } if(par2plot[p+1]){ plot(density(x@gamma2[ind]), xlab = expression(gamma^2), ...) if(priorMeans) abline(v = x@model$gamma2, col = col.priorMean, lty = lty.priorMean)
}
}

par(old.settings[which(names(old.settings) %in% names(par.options))])
})

#' Plot method for the Bayesian estimation results
#'
#' @description Plot method for the estimation results of the hierarchical (mixed) diffusion model.
#' @param x est.mixedDiffusion class, created with method \code{\link{estimate,mixedDiffusion-method}}
#' @param par.options list of options for function par()
#' @param style one out of "chains", "acf", "density", "int.phi"
#' @param par2plot logical vector, which parameters to be plotted, order: \eqn{(\mu, \Omega, \gamma^2)}
#' @param reduced logical (1), if TRUE, the chains are thinned and burn-in phase is dropped
#' @param thinning thinning rate, if missing, the proposed one by the estimation procedure is taken
#' @param burnIn burn-in phase, if missing, the proposed one by the estimation procedure is taken
#' @param priorMeans logical(1), if TRUE (default), prior means are marked with a line
#' @param col.priorMean color of the prior mean line, default 2
#' @param lty.priorMean linetype of the prior mean line, default 1
#' @param level level for style = "int.phi"
#' @param phi in the case of simulation study: known values for phi
#' @param ... optional plot parameters
#' @examples
#' mu <- c(10, 3, 1); Omega = c(1, 0.4, 0.01)
#' phi <- sapply(1:3, function(i) rnorm(20, mu[i], sqrt(Omega[i])))
#' model <- set.to.class("mixedDiffusion", b.fun = function(phi, t, y) phi[1]-phi[2]*y,
#'     parameter = list(mu = mu, Omega = Omega, phi = phi, gamma2 = 0.1),
#'     y0 = function(phi, t) phi[3], sT.fun = function(t, x) sqrt(abs(x)))
#' data <- simulate(model, t = seq(0, 1, by = 0.02), plot.series = TRUE)
#' est <- estimate(model, t = seq(0, 1, by = 0.02), data, 100)  # nMCMC small for example
#' plot(est, burnIn = 10, thinning = 2, reduced = TRUE)
#' plot(est, par.options = list(mar = c(5, 4.5, 4, 2) + 0.1, mfrow = c(2,1)), xlab = "iteration")
#' plot(est, style = "acf", main = "")
#' plot(est, style = "density", lwd = 2, priorMean = FALSE)
#' plot(est, style = "density", col.priorMean = 1, lty.priorMean = 2, main = "posterior")
#' plot(est, style = "acf", par.options = list(), main = "", par2plot = c(rep(FALSE, 6), TRUE))
#' plot(est, style = "int.phi", phi = phi, par2plot = c(TRUE, FALSE, FALSE))
#' @export
setMethod(f = "plot", signature = "est.mixedDiffusion",
definition = function(x, par.options, style = c("chains", "acf", "density", "int.phi"), par2plot, reduced = FALSE,
thinning, burnIn, priorMeans = TRUE, col.priorMean = 2, lty.priorMean = 1, level = 0.05, phi, ...) {

style <- match.arg(style)
if(reduced){
if(missing(thinning)) thinning <- x@thinning
if(missing(burnIn)) burnIn <- x@burnIn
ind <- seq(burnIn + 1, length(x@gamma2), by = thinning)
}else{
ind <- 1:length(x@gamma2)
}

p <- ncol(x@mu)
k <- nrow(x@phi[[1]])

if(missing(par2plot)) par2plot <- rep(TRUE, 1+2*p)
if(missing(par.options)){
if(sum(par2plot) == 1) par.options <- list(mfrow = c(1, 1))
if(sum(par2plot) > 1) par.options <- list(mfrow = c(ceiling(sum(par2plot)/2), 2))
if(style == "int.phi") par.options <- list(mfrow = c(sum(par2plot[1:p]), 1))
}

par(par.options)

if(style == "chains"){
for(i in 1:p) {
if(par2plot[i]){
plot(x@mu[ind, i], type = "l", ylab = bquote(mu[.(i)]), ...)
if(priorMeans) abline(h = x@model$mu[i], col = col.priorMean, lty = lty.priorMean) } } for(i in 1:p) { if(par2plot[p+i]){ plot(x@Omega[ind, i], type = "l", ylab = bquote(omega[.(i)]), ...) if(priorMeans) abline(h = x@model$Omega[i], col = col.priorMean, lty = lty.priorMean)
}
}
if(par2plot[2*p+1]){
plot(x@gamma2[ind], type = "l", ylab = expression(gamma^2), ...)
if(priorMeans) abline(h = x@model$gamma2, col = col.priorMean, lty = lty.priorMean) } } if(style == "acf"){ for(i in 1:p) { if(par2plot[i]) acf(x@mu[ind, i], xlab = bquote(mu[.(i)]), ...) } for(i in 1:p){ if(par2plot[p+i]) acf(x@Omega[ind, i], xlab = bquote(omega[.(i)]), ...) } if(par2plot[2*p+1]) acf(x@gamma2[ind], xlab = expression(gamma^2), ...) } if(style == "density"){ for(i in 1:p) { if(par2plot[i]){ plot(density(x@mu[ind, i]), xlab = bquote(mu[.(i)]), ...) if(priorMeans) abline(v = x@model$mu[i], col = col.priorMean, lty = lty.priorMean)
}
}
for(i in 1:p) {
if(par2plot[p+i]){
plot(density(x@Omega[ind, i]), xlab = bquote(omega[.(i)]), ...)
if(priorMeans) abline(v = x@model$Omega[i], col = col.priorMean, lty = lty.priorMean) } } if(par2plot[2*p+1]){ plot(density(x@gamma2[ind]), xlab = expression(gamma^2), ...) if(priorMeans) abline(v = x@model$gamma2, col = col.priorMean, lty = lty.priorMean)
}
}
if(style == "int.phi"){

qu <- lapply(1:p, function(i) {
apply( sapply(x@phi[ind], function(m) m[, i]), 1, function(ve) quantile(ve, c(level/2, 0.5, 1-level/2)) )
})
for(i in 1:p){
if(par2plot[i]){
plot(qu[[i]][2, ], pch = 20, ylim = range(qu[[i]]), ylab = bquote(phi[.(i)]), ...)
segments(1:k, qu[[i]][1, ], 1:k, qu[[i]][3, ], ...)
segments(1:k-0.1, qu[[i]][1, ], 1:k+0.1, qu[[i]][1, ], ...)
segments(1:k-0.1, qu[[i]][3, ], 1:k+0.1, qu[[i]][3, ], ...)
if(!missing(phi)) points(phi[, i], pch = 20, col = col.priorMean)
}
}
par(old.settings[which(names(old.settings) %in% names(par.options))])
return(invisible(qu))
}

par(old.settings[which(names(old.settings) %in% names(par.options))])
})

#' Plot method for the Bayesian estimation results
#'
#' @description Plot method for the estimation results of the hidden diffusion model.
#' @param x est.hiddenDiffusion class, created with method \code{\link{estimate,hiddenDiffusion-method}}
#' @param par.options list of options for function par()
#' @param style one out of "chains", "acf", "density"
#' @param par2plot logical vector, which parameters to be plotted, order: \eqn{(\phi, \gamma^2, \sigma^2, Y)}
#' @param reduced logical (1), if TRUE, the chains are thinned and burn-in phase is dropped
#' @param thinning thinning rate, if missing, the proposed one by the estimation procedure is taken
#' @param burnIn burn-in phase, if missing, the proposed one by the estimation procedure is taken
#' @param priorMeans logical(1), if TRUE (default), prior means are marked with a line
#' @param col.priorMean color of the prior mean line, default 2
#' @param lty.priorMean linetype of the prior mean line, default 1
#' @param ... optional plot parameters
#' @examples
#' model <- set.to.class("hiddenDiffusion", b.fun = function(phi, t, y) phi[1]-phi[2]*y,
#'     parameter = list(phi = c(10, 1), gamma2 = 1, sigma2 = 0.1),
#'     y0 = function(phi, t) 0.5)
#' data <- simulate(model, t = seq(0, 1, by = 0.01), plot.series = TRUE)
#' est <- estimate(model, t = seq(0, 1, by = 0.01), data$Y, 100) # nMCMC small for example #' plot(est) #' plot(est, par2plot = c(rep(FALSE, 3), TRUE, FALSE), ylim = c(0.001, 0.1), par.options = list()) #' plot(est, burnIn = 10, thinning = 2, reduced = TRUE) #' plot(est, par.options = list(mar = c(5, 4.5, 4, 2) + 0.1, mfrow = c(3,1)), xlab = "iteration") #' plot(est, style = "acf", main = "", par2plot = c(TRUE, TRUE, FALSE, FALSE)) #' plot(est, style = "density", lwd = 2, priorMean = FALSE) #' plot(est, style = "density", col.priorMean = 1, lty.priorMean = 2, main = "posterior") #' plot(est, style = "acf", par.options = list(), main = "", par2plot = c(FALSE, FALSE, TRUE, TRUE)) #' @export setMethod(f = "plot", signature = "est.hiddenDiffusion", definition = function(x, par.options, style = c("chains", "acf", "density"), par2plot, reduced = FALSE, thinning, burnIn, priorMeans = TRUE, col.priorMean = 2, lty.priorMean = 1, ...) { old.settings <- par(no.readonly = TRUE) style <- match.arg(style) if(reduced){ if(missing(thinning)) thinning <- x@thinning if(missing(burnIn)) burnIn <- x@burnIn ind <- seq(burnIn + 1, length(x@gamma2), by = thinning) }else{ ind <- 1:length(x@gamma2) } p <- ncol(x@phi) if(missing(par2plot)){ if(style == "chains") par2plot <- rep(TRUE, 3+p) if(style %in% c("acf", "density")) par2plot <- rep(TRUE, 2+p) } if(missing(par.options)){ if(sum(par2plot) == 1) par.options <- list(mfrow = c(1, 1)) if(sum(par2plot) > 1) par.options <- list(mfrow = c(ceiling(sum(par2plot)/2), 2)) } par(par.options) if(style == "chains"){ for(i in 1:p) { if(par2plot[i]){ plot(x@phi[ind, i], type = "l", ylab = bquote(phi[.(i)]), ...) if(priorMeans) abline(h = x@model$phi[i], col = col.priorMean, lty = lty.priorMean)
}
}
if(par2plot[p+1]){
plot(x@gamma2[ind], type = "l", ylab = expression(gamma^2), ...)
if(priorMeans) abline(h = x@model$gamma2, col = col.priorMean, lty = lty.priorMean) } if(par2plot[p+2]){ plot(x@sigma2[ind], type = "l", ylab = expression(sigma^2), ...) if(priorMeans) abline(h = x@model$sigma2, col = col.priorMean, lty = lty.priorMean)
}
if(par2plot[p+3]){
plot(x@t, x@Z, pch = 20, xlab = "t", ylab = "data", col = col.priorMean, ylim = range(c(x@Z, range(x@Y.est[ind, ]))))
for(i in ind) lines(x@t, x@Y.est[i,])
points(x@t, x@Z, pch = 20, col = col.priorMean)
}
}
if(style == "acf"){
for(i in 1:p) {
if(par2plot[i]){
acf(x@phi[ind, i], xlab = bquote(phi[.(i)]), ...)
}
}
if(par2plot[p+1]) acf(x@gamma2[ind], xlab = expression(gamma^2), ...)
if(par2plot[p+2]) acf(x@sigma2[ind], xlab = expression(sigma^2), ...)
}
if(style == "density"){
for(i in 1:p) {
if(par2plot[i]){
plot(density(x@phi[ind, i]), xlab = bquote(phi[.(i)]), ...)
if(priorMeans) abline(v = x@model$phi[i], col = col.priorMean, lty = lty.priorMean) } } if(par2plot[p+1]){ plot(density(x@gamma2[ind]), xlab = expression(gamma^2), ...) if(priorMeans) abline(v = x@model$gamma2, col = col.priorMean, lty = lty.priorMean)
}
if(par2plot[p+2]){
plot(density(x@sigma2[ind]), xlab = expression(sigma^2), ...)
if(priorMeans) abline(v = x@model$sigma2, col = col.priorMean, lty = lty.priorMean) } } par(old.settings[which(names(old.settings) %in% names(par.options))]) }) #' Plot method for the Bayesian estimation results #' #' @description Plot method for the estimation results of the hidden hierarchical diffusion model. #' @param x est.hiddenmixedDiffusion class, created with method \code{\link{estimate,hiddenmixedDiffusion-method}} #' @param par.options list of options for function par() #' @param style one out of "chains", "acf", "density", "int.phi" #' @param par2plot logical vector, which parameters to be plotted, order: \eqn{(\mu, \Omega, \gamma^2, \sigma^2, Y)} #' @param reduced logical (1), if TRUE, the chains are thinned and burn-in phase is dropped #' @param thinning thinning rate, if missing, the proposed one by the estimation procedure is taken #' @param burnIn burn-in phase, if missing, the proposed one by the estimation procedure is taken #' @param priorMeans logical(1), if TRUE (default), prior means are marked with a line #' @param col.priorMean color of the prior mean line, default 2 #' @param lty.priorMean linetype of the prior mean line, default 1 #' @param level level for style = "int.phi" #' @param phi in the case of simulation study: known values for phi #' @param ... optional plot parameters #' @examples #' \dontrun{ #' mu <- c(10, 3, 1); Omega = c(1, 0.4, 0.01) #' phi <- sapply(1:3, function(i) rnorm(20, mu[i], sqrt(Omega[i]))) #' model <- set.to.class("hiddenmixedDiffusion", b.fun = function(phi, t, y) phi[1]-phi[2]*y, #' parameter = list(mu = mu, Omega = Omega, phi = phi, gamma2 = 1, sigma2 = 0.1), #' y0 = function(phi, t) phi[3]) #' data <- simulate(model, t = seq(0, 1, by = 0.02), plot.series = TRUE) #' est <- estimate(model, t = seq(0, 1, by = 0.02), data$Z, 1000)
#' plot(est, burnIn = 10, thinning = 2, reduced = TRUE)
#' plot(est, par.options = list(mar = c(5, 4.5, 4, 2) + 0.1, mfrow = c(2,1)), xlab = "iteration")
#' plot(est, style = "acf", main = "", par2plot = c(TRUE, TRUE, rep(FALSE, 7)))
#' plot(est, style = "density", lwd = 2, priorMean = FALSE,
#'    par2plot = c(rep(FALSE, 6), TRUE, TRUE, FALSE))
#' plot(est, style = "density", col.priorMean = 1, lty.priorMean = 2, main = "posterior")
#' plot(est, style = "acf", par.options = list(), main = "", par2plot = c(rep(FALSE, 6), TRUE, TRUE))
#' plot(est, style = "int.phi", phi = phi, par2plot = c(TRUE, FALSE, FALSE))
#' }
#' @export
setMethod(f = "plot", signature = "est.hiddenmixedDiffusion",
definition = function(x, par.options, style = c("chains", "acf", "density", "int.phi"), par2plot, reduced = FALSE,
thinning, burnIn, priorMeans = TRUE, col.priorMean = 2, lty.priorMean = 1, level = 0.05, phi, ...) {

style <- match.arg(style)
if(reduced){
if(missing(thinning)) thinning <- x@thinning
if(missing(burnIn)) burnIn <- x@burnIn
ind <- seq(burnIn + 1, length(x@gamma2), by = thinning)
}else{
ind <- 1:length(x@gamma2)
}

p <- ncol(x@mu)
k <- nrow(x@phi[[1]])

if(missing(par2plot)){
if(style == "chains")  par2plot <- rep(TRUE, 3+2*p)
if(style %in% c("acf", "density"))  par2plot <- rep(TRUE, 2+2*p)
if(style == "int.phi") par2plot <- rep(TRUE, p)
}
if(missing(par.options)){
if(sum(par2plot) == 1) par.options <- list(mfrow = c(1, 1))
if(sum(par2plot) > 1) par.options <- list(mfrow = c(ceiling(sum(par2plot)/2), 2))

if(style == "int.phi") par.options <- list(mfrow = c(sum(par2plot[1:p]), 1))
}

par(par.options)

if(style == "chains"){
for(i in 1:p) {
if(par2plot[i]){
plot(x@mu[ind, i], type = "l", ylab = bquote(mu[.(i)]), ...)
if(priorMeans) abline(h = x@model$mu[i], col = col.priorMean, lty = lty.priorMean) } } for(i in 1:p) { if(par2plot[p+i]){ plot(x@Omega[ind, i], type = "l", ylab = bquote(omega[.(i)]), ...) if(priorMeans) abline(h = x@model$Omega[i], col = col.priorMean, lty = lty.priorMean)
}
}
if(par2plot[2*p+1]){
plot(x@gamma2[ind], type = "l", ylab = expression(gamma^2), ...)
if(priorMeans) abline(h = x@model$gamma2, col = col.priorMean, lty = lty.priorMean) } if(par2plot[2*p+2]){ plot(x@sigma2[ind], type = "l", ylab = expression(sigma^2), ...) if(priorMeans) abline(h = x@model$sigma2, col = col.priorMean, lty = lty.priorMean)
}
if(par2plot[2*p+3]){
plot(x@t, x@Z[1,], pch = 20, xlab = "t", ylab = "first data series", col = col.priorMean,
ylim = range(c(x@Z[1,], range(sapply(ind, function(i) x@Y.est[[i]][[1]])))))
for(i in ind) lines(x@t, x@Y.est[[i]][[1]])
points(x@t, x@Z[1,], pch = 20, col = col.priorMean)
}
}
if(style == "acf"){
for(i in 1:p) {
if(par2plot[i]) acf(x@mu[ind, i], xlab = bquote(mu[.(i)]), ...)
}
for(i in 1:p){
if(par2plot[p+i]) acf(x@Omega[ind, i], xlab = bquote(omega[.(i)]), ...)
}
if(par2plot[2*p+1]) acf(x@gamma2[ind], xlab = expression(gamma^2), ...)
if(par2plot[2*p+2]) acf(x@sigma2[ind], xlab = expression(sigma^2), ...)
}
if(style == "density"){
for(i in 1:p) {
if(par2plot[i]){
plot(density(x@mu[ind, i]), xlab = bquote(mu[.(i)]), ...)
if(priorMeans) abline(v = x@model$mu[i], col = col.priorMean, lty = lty.priorMean) } } for(i in 1:p) { if(par2plot[p+i]){ plot(density(x@Omega[ind, i]), xlab = bquote(omega[.(i)]), ...) if(priorMeans) abline(v = x@model$Omega[i], col = col.priorMean, lty = lty.priorMean)
}
}
if(par2plot[2*p+1]){
plot(density(x@gamma2[ind]), xlab = expression(gamma^2), ...)
if(priorMeans) abline(v = x@model$gamma2, col = col.priorMean, lty = lty.priorMean) } if(par2plot[2*p+2]){ plot(density(x@sigma2[ind]), xlab = expression(sigma^2), ...) if(priorMeans) abline(v = x@model$sigma2, col = col.priorMean, lty = lty.priorMean)
}
}
if(style == "int.phi"){

qu <- lapply(1:p, function(i) {
apply( sapply(x@phi[ind], function(m) m[, i]), 1, function(ve) quantile(ve, c(level/2, 0.5, 1-level/2)) )
})
for(i in 1:p){
if(par2plot[i]){
plot(qu[[i]][2, ], pch = 20, ylim = range(qu[[i]]), ylab = bquote(phi[.(i)]), ...)
segments(1:k, qu[[i]][1, ], 1:k, qu[[i]][3, ], ...)
segments(1:k-0.1, qu[[i]][1, ], 1:k+0.1, qu[[i]][1, ], ...)
segments(1:k-0.1, qu[[i]][3, ], 1:k+0.1, qu[[i]][3, ], ...)
if(!missing(phi)) points(phi[, i], pch = 20, col = col.priorMean)
}
}
par(old.settings[which(names(old.settings) %in% names(par.options))])
return(invisible(qu))
}

par(old.settings[which(names(old.settings) %in% names(par.options))])
})

#' Plot method for the Bayesian estimation results
#'
#' @description Plot method for the estimation results of the jump regression model.
#' @param x est.jumpRegression class, created with method \code{\link{estimate,jumpRegression-method}}
#' @param par.options list of options for function par()
#' @param style one out of "chains", "acf", "density"
#' @param par2plot logical vector, which parameters to be plotted, order: \eqn{(\phi, \theta, \gamma^2, \xi, N)}
#' @param reduced logical (1), if TRUE, the chains are thinned and burn-in phase is dropped
#' @param thinning thinning rate, if missing, the proposed one by the estimation procedure is taken
#' @param burnIn burn-in phase, if missing, the proposed one by the estimation procedure is taken
#' @param priorMeans logical(1), if TRUE (default), prior means are marked with a line
#' @param col.priorMean color of the prior mean line, default 2
#' @param lty.priorMean linetype of the prior mean line, default 1
#' @param ... optional plot parameters
#' @examples
#' 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 = seq(0, 1, by = 0.01), plot.series = TRUE)
#' est <- estimate(model, t = seq(0, 1, by = 0.01), data, 1000)  # nMCMC small for example
#' plot(est)
#' plot(est, burnIn = 100, thinning = 2, reduced = TRUE)
#' plot(est, par.options = list(mar = c(5, 4.5, 4, 2) + 0.1, mfrow = c(2, 3)), xlab = "iteration")
#' plot(est, style = "acf", main = "", par2plot = c(TRUE, FALSE, FALSE, TRUE, TRUE))
#' plot(est, style = "density", lwd = 2, priorMean = FALSE)
#' plot(est, style = "density", col.priorMean = 1, lty.priorMean = 2, main = "posterior")
#' plot(est, style = "acf", par.options = list(), par2plot = c(TRUE, rep(FALSE, 4)), main = "")
#' @export
setMethod(f = "plot", signature = "est.jumpRegression",
definition = function(x, par.options, style = c("chains", "acf", "density"), par2plot, reduced = FALSE,
thinning, burnIn, priorMeans = TRUE, col.priorMean = 2, lty.priorMean = 1, ...) {

style <- match.arg(style)
if(reduced){
if(missing(thinning)) thinning <- x@thinning
if(missing(burnIn)) burnIn <- x@burnIn
ind <- seq(burnIn + 1, length(x@gamma2), by = thinning)
}else{
ind <- 1:length(x@gamma2)
}

p1 <- ncol(x@theta)
p2 <- ncol(x@xi)
if(style == "chains"){
p3 <- c(1, 0)[c(sum(dim(x@N.est)) > 0, sum(dim(x@N.est)) == 0)]
} else {
p3 <- 0
}
if(missing(par2plot)) par2plot <- rep(TRUE, p1+p2+p3+1)
if(missing(par.options)){
if(sum(par2plot) == 1) par.options <- list(mfrow = c(1, 1))
if(sum(par2plot) > 1) par.options <- list(mfrow = c(ceiling(sum(par2plot)/2), 2))
}

par(par.options)

if(style == "chains"){

for(i in 1:p1) {
if(par2plot[i]){
plot(x@theta[ind, i], type = "l", ylab = bquote(theta[.(i)]), ...)
if(priorMeans) abline(h = x@model$theta[i], col = col.priorMean, lty = lty.priorMean) } } if(par2plot[p1 + 1]){ plot(x@gamma2[ind], type = "l", ylab = expression(gamma^2), ...) if(priorMeans) abline(h = x@model$gamma2, col = col.priorMean, lty = lty.priorMean)
}
for(i in 1:p2){
if(par2plot[p1 + 1 + i]){
plot(x@xi[ind, i], type = "l", ylab = bquote(xi[.(i)]), ...)
if(priorMeans)  abline(h = x@model$xi[i], col = col.priorMean, lty = lty.priorMean) } } if(sum(dim(x@N.est) > 0)){ if(par2plot[p1 +p2 + 1]){ plot(x@t, x@N.est[,length(x@gamma2)], type = "l", xlab = "t", ylab = "N", ylim = range(x@N.est[,ind]), ...) for(i in ind) lines(x@t, x@N.est[,i]) } } } if(style == "acf"){ for(i in 1:p1) { if(par2plot[i]) acf(x@theta[ind, i], xlab = bquote(theta[.(i)]), ...) } if(par2plot[p1 + 1]) acf(x@gamma2[ind], xlab = expression(gamma^2), ...) for(i in 1:p2){ if(par2plot[p1 + 1 + + i]) acf(x@xi[ind, i], xlab = bquote(xi[.(i)]), ...) } } if(style == "density"){ for(i in 1:p1) { if(par2plot[i]){ plot(density(x@theta[ind, i]), xlab = bquote(theta[.(i)]), ...) if(priorMeans) abline(v = x@model$theta[i], col = col.priorMean, lty = lty.priorMean)
}
}
if(par2plot[p1 + 1]){
plot(density(x@gamma2[ind]), xlab = expression(gamma^2), ...)
if(priorMeans)  abline(v = x@model$gamma2, col = col.priorMean, lty = lty.priorMean) } for(i in 1:p2){ if(par2plot[p1 + 1 + i]){ plot(density(x@xi[ind, i]), xlab = bquote(xi[.(i)]), ...) if(priorMeans) abline(v = x@model$xi[i], col = col.priorMean, lty = lty.priorMean)
}
}

}

par(old.settings[which(names(old.settings) %in% names(par.options))])
})

#' Plot method for the Bayesian estimation results
#'
#' @description Plot method for the estimation results of the NHPP.
#' @param x est.NHPP class, created with method \code{\link{estimate,NHPP-method}}
#' @param par.options list of options for function par()
#' @param style one out of "chains", "acf", "density"
#' @param par2plot logical vector, which parameters to be plotted, order: \eqn{(\phi, \theta, \gamma^2, \xi, N)}
#' @param reduced logical (1), if TRUE, the chains are thinned and burn-in phase is dropped
#' @param thinning thinning rate, if missing, the proposed one by the estimation procedure is taken
#' @param burnIn burn-in phase, if missing, the proposed one by the estimation procedure is taken
#' @param priorMeans logical(1), if TRUE (default), prior means are marked with a line
#' @param col.priorMean color of the prior mean line, default 2
#' @param lty.priorMean linetype of the prior mean line, default 1
#' @param ... optional plot parameters
#' @examples
#' model <- set.to.class("NHPP", parameter = list(xi = c(5, 1/2)),
#'   Lambda = function(t, xi) (t/xi[2])^xi[1])
#' data <- simulate(model, t = seq(0, 1, by = 0.01), plot.series = TRUE)
#' est <- estimate(model, t = seq(0, 1, by = 0.01), data$Times, 10000) # nMCMC small for example #' plot(est) #' plot(est, burnIn = 1000, thinning = 20, reduced = TRUE) #' plot(est, xlab = "iteration") #' plot(est, style = "acf", main = "", par2plot = c(TRUE, FALSE), par.options = list(mfrow = c(1, 1))) #' plot(est, style = "density", lwd = 2, priorMean = FALSE) #' plot(est, style = "density", col.priorMean = 1, lty.priorMean = 2, main = "posterior") #' plot(est, style = "acf", par.options = list(), par2plot = c(FALSE, TRUE), main = "") #' @export setMethod(f = "plot", signature = "est.NHPP", definition = function(x, par.options, style = c("chains", "acf", "density"), par2plot, reduced = FALSE, thinning, burnIn, priorMeans = TRUE, col.priorMean = 2, lty.priorMean = 1, ...) { old.settings <- par(no.readonly = TRUE) style <- match.arg(style) if(reduced){ if(missing(thinning)) thinning <- x@thinning if(missing(burnIn)) burnIn <- x@burnIn ind <- seq(burnIn + 1, nrow(x@xi), by = thinning) }else{ ind <- 1:nrow(x@xi) } p <- ncol(x@xi) if(missing(par2plot)) par2plot <- rep(TRUE, p) if(missing(par.options)){ if(sum(par2plot) == 1) par.options <- list(mfrow = c(1, 1)) if(sum(par2plot) > 1) par.options <- list(mfrow = c(ceiling(sum(par2plot)/2), 2)) } par(par.options) if(style == "chains"){ for(i in 1:p){ if(par2plot[i]){ plot(x@xi[ind, i], type = "l", ylab = bquote(xi[.(i)]), ...) if(priorMeans) abline(h = x@model$xi[i], col = col.priorMean, lty = lty.priorMean)
}
}

}
if(style == "acf"){

for(i in 1:p){
if(par2plot[i]) acf(x@xi[ind, i], xlab = bquote(xi[.(i)]), ...)
}
}
if(style == "density"){

for(i in 1:p){
if(par2plot[i]){
plot(density(x@xi[ind, i]), xlab = bquote(xi[.(i)]), ...)
if(priorMeans)  abline(v = x@model$xi[i], col = col.priorMean, lty = lty.priorMean) } } } par(old.settings[which(names(old.settings) %in% names(par.options))]) }) #' Plot method for the Bayesian estimation results #' #' @description Plot method for the estimation results of the regression model. #' @param x est.Regression class, created with method \code{\link{estimate,Regression-method}} #' @param par.options list of options for function par() #' @param style one out of "chains", "acf", "density" #' @param par2plot logical vector, which parameters to be plotted, order: \eqn{(\phi, \gamma^2)} #' @param reduced logical (1), if TRUE, the chains are thinned and burn-in phase is dropped #' @param thinning thinning rate, if missing, the proposed one by the estimation procedure is taken #' @param burnIn burn-in phase, if missing, the proposed one by the estimation procedure is taken #' @param priorMeans logical(1), if TRUE (default), prior means are marked with a line #' @param col.priorMean color of the prior mean line, default 2 #' @param lty.priorMean linetype of the prior mean line, default 1 #' @param ... optional plot parameters #' @examples #' 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 = seq(0, 1, by = 0.01), plot.series = TRUE) #' est <- estimate(model, t = seq(0, 1, by = 0.01), data, 1000) # nMCMC small for example #' plot(est) #' plot(est, burnIn = 100, thinning = 2, reduced = TRUE) #' plot(est, par.options = list(mar = c(5, 4.5, 4, 2) + 0.1, mfrow = c(3,1)), xlab = "iteration") #' plot(est, style = "acf", main = "", par2plot = c(TRUE, TRUE, FALSE)) #' plot(est, style = "density", lwd = 2, priorMean = FALSE) #' plot(est, style = "density", col.priorMean = 1, lty.priorMean = 2, main = "posterior") #' plot(est, style = "acf", par.options = list(), main = "", par2plot = c(FALSE, FALSE, TRUE)) #' @export setMethod(f = "plot", signature = "est.Regression", definition = function(x, par.options, style = c("chains", "acf", "density"), par2plot, reduced = FALSE, thinning, burnIn, priorMeans = TRUE, col.priorMean = 2, lty.priorMean = 1, ...) { old.settings <- par(no.readonly = TRUE) style <- match.arg(style) if(reduced){ if(missing(thinning)) thinning <- x@thinning if(missing(burnIn)) burnIn <- x@burnIn ind <- seq(burnIn + 1, length(x@gamma2), by = thinning) }else{ ind <- 1:length(x@gamma2) } p <- ncol(x@phi) if(missing(par2plot)) par2plot <- rep(TRUE, 1+p) if(missing(par.options)){ if(sum(par2plot) == 1) par.options <- list(mfrow = c(1, 1)) if(sum(par2plot) > 1) par.options <- list(mfrow = c(ceiling(sum(par2plot)/2), 2)) } par(par.options) if(style == "chains"){ for(i in 1:p) { if(par2plot[i]){ plot(x@phi[ind, i], type = "l", ylab = bquote(phi[.(i)]), ...) if(priorMeans) abline(h = x@model$phi[i], col = col.priorMean, lty = lty.priorMean)
}
}
if(par2plot[p+1]){
plot(x@gamma2[ind], type = "l", ylab = expression(gamma^2), ...)
if(priorMeans) abline(h = x@model$gamma2, col = col.priorMean, lty = lty.priorMean) } } if(style == "acf"){ for(i in 1:p) { if(par2plot[i]){ acf(x@phi[ind, i], xlab = bquote(phi[.(i)]), ...) } } if(par2plot[p+1]) acf(x@gamma2[ind], xlab = expression(gamma^2), ...) } if(style == "density"){ for(i in 1:p) { if(par2plot[i]){ plot(density(x@phi[ind, i]), xlab = bquote(phi[.(i)]), ...) if(priorMeans) abline(v = x@model$phi[i], col = col.priorMean, lty = lty.priorMean)
}
}
if(par2plot[p+1]){
plot(density(x@gamma2[ind]), xlab = expression(gamma^2), ...)
if(priorMeans) abline(v = x@model$gamma2, col = col.priorMean, lty = lty.priorMean) } } par(old.settings[which(names(old.settings) %in% names(par.options))]) }) #' Plot method for the Bayesian estimation results #' #' @description Plot method for the estimation results of the hierarchical (mixed) regression model. #' @param x est.mixedRegression class, created with method \code{\link{estimate,mixedRegression-method}} #' @param par.options list of options for function par() #' @param style one out of "chains", "acf", "density", "int.phi" #' @param par2plot logical vector, which parameters to be plotted, order: \eqn{(\mu, \Omega, \gamma^2)} #' @param reduced logical (1), if TRUE, the chains are thinned and burn-in phase is dropped #' @param thinning thinning rate, if missing, the proposed one by the estimation procedure is taken #' @param burnIn burn-in phase, if missing, the proposed one by the estimation procedure is taken #' @param priorMeans logical(1), if TRUE (default), prior means are marked with a line #' @param col.priorMean color of the prior mean line, default 2 #' @param lty.priorMean linetype of the prior mean line, default 1 #' @param level level for style = "int.phi" #' @param phi in the case of simulation study: known values for phi #' @param ... optional plot parameters #' @examples #' mu <- c(1, 3); Omega = c(0.4, 0.01) #' phi <- sapply(1:2, function(i) rnorm(20, mu[i], sqrt(Omega[i]))) #' model <- set.to.class("mixedRegression", fun = function(phi, t) phi[1]*t + phi[2], #' parameter = list(mu = mu, Omega = Omega, phi = phi, gamma2 = 0.1)) #' data <- simulate(model, t = seq(0, 1, by = 0.02), plot.series = TRUE) #' est <- estimate(model, t = seq(0, 1, by = 0.02), data, 100) # nMCMC small for example #' plot(est, burnIn = 10, thinning = 2, reduced = TRUE) #' plot(est, par.options = list(mar = c(5, 4.5, 4, 2) + 0.1, mfrow = c(2,1)), xlab = "iteration") #' plot(est, style = "acf", main = "") #' plot(est, style = "density", lwd = 2, priorMean = FALSE) #' plot(est, style = "density", col.priorMean = 1, lty.priorMean = 2, main = "posterior") #' plot(est, style = "acf", par.options = list(), main = "", par2plot = c(rep(FALSE, 4), TRUE)) #' plot(est, style = "int.phi", phi = phi, par2plot = c(TRUE, FALSE)) #' @export setMethod(f = "plot", signature = "est.mixedRegression", definition = function(x, par.options, style = c("chains", "acf", "density", "int.phi"), par2plot, reduced = FALSE, thinning, burnIn, priorMeans = TRUE, col.priorMean = 2, lty.priorMean = 1, level = 0.05, phi, ...) { old.settings <- par(no.readonly = TRUE) style <- match.arg(style) if(reduced){ if(missing(thinning)) thinning <- x@thinning if(missing(burnIn)) burnIn <- x@burnIn ind <- seq(burnIn + 1, length(x@gamma2), by = thinning) }else{ ind <- 1:length(x@gamma2) } p <- ncol(x@mu) k <- nrow(x@phi[[1]]) if(missing(par2plot)) par2plot <- rep(TRUE, 1+2*p) if(missing(par.options)){ if(sum(par2plot) == 1) par.options <- list(mfrow = c(1, 1)) if(sum(par2plot) > 1) par.options <- list(mfrow = c(ceiling(sum(par2plot)/2), 2)) if(style == "int.phi") par.options <- list(mfrow = c(sum(par2plot[1:p]), 1)) } par(par.options) if(style == "chains"){ for(i in 1:p) { if(par2plot[i]){ plot(x@mu[ind, i], type = "l", ylab = bquote(mu[.(i)]), ...) if(priorMeans) abline(h = x@model$mu[i], col = col.priorMean, lty = lty.priorMean)
}
}
for(i in 1:p) {
if(par2plot[p+i]){
plot(x@Omega[ind, i], type = "l", ylab = bquote(omega[.(i)]), ...)
if(priorMeans) abline(h = x@model$Omega[i], col = col.priorMean, lty = lty.priorMean) } } if(par2plot[2*p+1]){ plot(x@gamma2[ind], type = "l", ylab = expression(gamma^2), ...) if(priorMeans) abline(h = x@model$gamma2, col = col.priorMean, lty = lty.priorMean)
}
}
if(style == "acf"){
for(i in 1:p) {
if(par2plot[i]) acf(x@mu[ind, i], xlab = bquote(mu[.(i)]), ...)
}
for(i in 1:p){
if(par2plot[p+i]) acf(x@Omega[ind, i], xlab = bquote(omega[.(i)]), ...)
}
if(par2plot[2*p+1]) acf(x@gamma2[ind], xlab = expression(gamma^2), ...)
}
if(style == "density"){
for(i in 1:p) {
if(par2plot[i]){
plot(density(x@mu[ind, i]), xlab = bquote(mu[.(i)]), ...)
if(priorMeans) abline(v = x@model$mu[i], col = col.priorMean, lty = lty.priorMean) } } for(i in 1:p) { if(par2plot[p+i]){ plot(density(x@Omega[ind, i]), xlab = bquote(omega[.(i)]), ...) if(priorMeans) abline(v = x@model$Omega[i], col = col.priorMean, lty = lty.priorMean)
}
}
if(par2plot[2*p+1]){
plot(density(x@gamma2[ind]), xlab = expression(gamma^2), ...)
if(priorMeans) abline(v = x@model\$gamma2, col = col.priorMean, lty = lty.priorMean)
}
}
if(style == "int.phi"){

qu <- lapply(1:p, function(i) {
apply( sapply(x@phi[ind], function(m) m[, i]), 1, function(ve) quantile(ve, c(level/2, 0.5, 1-level/2)) )
})
for(i in 1:p){
if(par2plot[i]){
plot(qu[[i]][2, ], pch = 20, ylim = range(qu[[i]]), ylab = bquote(phi[.(i)]), ...)
segments(1:k, qu[[i]][1, ], 1:k, qu[[i]][3, ], ...)
segments(1:k-0.1, qu[[i]][1, ], 1:k+0.1, qu[[i]][1, ], ...)
segments(1:k-0.1, qu[[i]][3, ], 1:k+0.1, qu[[i]][3, ], ...)
if(!missing(phi)) points(phi[, i], pch = 20, col = col.priorMean)
}
}
par(old.settings[which(names(old.settings) %in% names(par.options))])
return(invisible(qu))
}

par(old.settings[which(names(old.settings) %in% names(par.options))])
})

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