R/fitted.bsam.R

In bsamGP: Bayesian Spectral Analysis Models using Gaussian Process Priors

"fitted.bsam" <- function(object, alpha = 0.05, HPD = TRUE, ...) {
  smcmc <- object$mcmc$smcmc
  n <- object$n
  nint <- object$nint + 1
  nfun <- object$nfun

  fxobsg <- object$fit.draws$fxobs
  fxgridg <- object$fit.draws$fxgrid
  wbg <- object$fit.draws$wbeta
  if (object$model != "gbsar") {
    yhatg <- object$fit.draws$yhat
  } else {
    yhatg <- object$fit.draws$muhat
  }
  fxobs <- list()
  fxobsm <- apply(fxobsg, c(1, 2), mean)
  fxobs$mean <- fxobsm

  fxgrid <- list()
  fxgridm <- apply(fxgridg, c(1, 2), mean)
  fxgrid$mean <- fxgridm

  wbeta <- list()
  wbm <- apply(wbg, 2, mean)
  wbeta$mean <- wbm

  yhat <- list()
  ym <- apply(yhatg, 2, mean)
  yhat$mean <- ym

  if (HPD) {
    prob <- 1 - alpha

    fx.l <- fx.u <- matrix(0, n, nfun)
    fxg.l <- fxg.u <- matrix(0, nint, nfun)
    for (i in 1:nfun) {
      fxobsg.o <- apply(fxobsg[, i, ], 1, sort)
      gap <- max(1, min(smcmc - 1, round(smcmc * prob)))
      init <- 1:(smcmc - gap)
      inds <- apply(fxobsg.o[init + gap, , drop = FALSE] - fxobsg.o[init, , drop = FALSE], 2, which.min)
      fx.l[, i] <- fxobsg.o[cbind(inds, 1:n)]
      fx.u[, i] <- fxobsg.o[cbind(inds + gap, 1:n)]

      fxgridg.o <- apply(fxgridg[, i, ], 1, sort)
      gap <- max(1, min(smcmc - 1, round(smcmc * prob)))
      init <- 1:(smcmc - gap)
      inds <- apply(fxgridg.o[init + gap, , drop = FALSE] - fxgridg.o[init, , drop = FALSE], 2, which.min)
      fxg.l[, i] <- fxgridg.o[cbind(inds, 1:nint)]
      fxg.u[, i] <- fxgridg.o[cbind(inds + gap, 1:nint)]
    }
    fxobs$lower <- fx.l
    fxobs$upper <- fx.u
    fxgrid$lower <- fxg.l
    fxgrid$upper <- fxg.u

    wbg.o <- apply(wbg, 2, sort)
    gap <- max(1, min(smcmc - 1, round(smcmc * prob)))
    init <- 1:(smcmc - gap)
    inds <- apply(wbg.o[init + gap, , drop = FALSE] - wbg.o[init, , drop = FALSE], 2, which.min)
    wbeta$lower <- wbg.o[cbind(inds, 1:n)]
    wbeta$upper <- wbg.o[cbind(inds + gap, 1:n)]

    yhatg.o <- apply(yhatg, 2, sort)
    gap <- max(1, min(smcmc - 1, round(smcmc * prob)))
    init <- 1:(smcmc - gap)
    inds <- apply(yhatg.o[init + gap, , drop = FALSE] - yhatg.o[init, , drop = FALSE], 2, which.min)
    yhat$lower <- yhatg.o[cbind(inds, 1:n)]
    yhat$upper <- yhatg.o[cbind(inds + gap, 1:n)]
  } else {
    fxobs$lower <- apply(fxobsg, c(1, 2), function(x) quantile(x, prob = alpha/2))
    fxobs$upper <- apply(fxobsg, c(1, 2), function(x) quantile(x, prob = 1 - alpha/2))

    fxgrid$lower <- apply(fxgridg, c(1, 2), function(x) quantile(x, prob = alpha/2))
    fxgrid$upper <- apply(fxgridg, c(1, 2), function(x) quantile(x, prob = 1 - alpha/2))

    wbeta$lower <- apply(wbg, 2, function(x) quantile(x, prob = alpha/2))
    wbeta$upper <- apply(wbg, 2, function(x) quantile(x, prob = 1 - alpha/2))

    yhat$lower <- apply(yhatg, 2, function(x) quantile(x, prob = alpha/2))
    yhat$upper <- apply(yhatg, 2, function(x) quantile(x, prob = 1 - alpha/2))
  }

  out <- object
  out$alpha <- alpha
  out$HPD <- HPD
  out$yhat <- yhat
  out$wbeta <- wbeta
  out$fxobs <- fxobs
  out$fxgrid <- fxgrid
  out$xgrid <- object$fit.draws$xgrid
  class(out) <- "fitted.bsam"
  out
}