R/final.resamp9.R

In milokmilo/Stranded: Calculate Biological Parameters from Stranded Cetaceans

#' Final resampling
#'
#' Final resampling. Modified from final.resamp.R of HPbayes_0.1 package
#' @param K
#' @param B1
#' @param H.new
#' @param H.k
#' @param log.like
#' @param d.keep
#' @param prior
#' @param h.mu
#' @param h.sig
#' @param nrisk
#' @param ndeath
#' @param B. Defalt = 400
#' @param theta.dim. Default = 9
#' @param age. Default = c(1e-05, 1, seq(5, 100, 5))
#' @keywords resampling
#' @export
#' @examples
#' final.resamp9() 

final.resamp9 <- function (K, B1, H.new, H.k, log.like, d.keep, prior, h.mu, h.sig, 
                           nrisk, ndeath, B = 400, theta.dim = 9, age = c(1e-05, 1, 
                                                                          seq(5, 100, 5))) {
  lx <- nrisk
  dx <- ndeath
  mp8.mle <- function(theta, x.fit = age) {
    p.hat <- mod9p(theta = theta, x = x.fit)
    ll <- ll.binom(x = dx, n = lx, p = p.hat)
    return(ll)
  }
  B0 <- theta.dim * 1000
  prior.cov <- cov(prior)
  imp.keep <- theta.dim * 100
  vwts <- NULL
  ewts <- NULL
  mwts <- NULL
  nup <- NULL
  nck <- NULL
  mwt.case <- NULL
  frac.nup <- NULL
  prilo <- apply(prior, 2, min)
  prihi <- apply(prior, 2, max)
  lo <- prilo
  hi <- prihi
  for (k in 1:K) {
    if (k == 1) {
      log.like.k <- rep(NA, B1)
      for (i in 1:B1) {
        log.like.k[i] <- mp8.mle(H.new[i, ])
        if (i%%100 == 0) {
        }
      }
    }
    if (k > 1) {
      log.like.k <- rep(NA, B)
      for (i in 1:B) {
        log.like.k[i] <- mp8.mle(H.new[i, ])
        if (i%%100 == 0) {
        }
      }
      H.k <- rbind(H.k, H.new)
    }
    log.like <- c(log.like, log.like.k)
    n.k <- B * (k - 1) + B0 + B1
    ncol.qk <- 1 + d.keep + (k - 1)
    DENS.MAT <- matrix(NA, nrow = n.k, ncol = ncol.qk)
    DENS.MAT[, 1] <- (B0/n.k) * dens.prior9(H.k, pri.lo = prilo, pri.hi = prihi)
    for (i in 1:(ncol.qk - 1)) {
      DENS.MAT[, (i + 1)] <- B/n.k * dmvnorm(x = H.k, mean = h.mu[i, 
                                                                  ], sigma = h.sig[, , i])
    }
    qk <- rowSums(DENS.MAT)
    a <- -max(log.like, na.rm = TRUE)
    like.k <- exp(log.like)
    wts.k.0 <- like.k * dens.prior9(H.k, pri.lo = prilo, pri.hi = prihi)
    wts.k <- wts.k.0/qk
    wts.k <- wts.k/sum(wts.k, na.rm = TRUE)
    wts.k[is.na(wts.k)] <- 0
    maxw.k <- max(wts.k)
    thresh.k <- min(maxw.k, 0.05)
    imp.theta <- H.k[wts.k == maxw.k, ]
    dist.to.mean <- mahalanobis(x = H.k, center = imp.theta, 
                                cov = prior.cov)
    keep.4.cov <- sort(dist.to.mean, decreasing = FALSE)[imp.keep]
    inputs.4.cov <- H.k[dist.to.mean <= keep.4.cov, ]
    wts.4.cov <- wts.k[dist.to.mean <= keep.4.cov]
    wts.4.cov <- ((wts.4.cov + 1/length(wts.k)) * 0.5)/sum(((wts.4.cov + 
                                                               1/length(wts.k)) * 0.5))
    Sigma.k <- cov.wt(inputs.4.cov, wt = wts.4.cov)$cov
    h.mu <- rbind(h.mu, imp.theta)
    h.sig[, , (d.keep + k)] <- Sigma.k
    H.new <- mvrnorm(n = B, mu = imp.theta, Sigma = Sigma.k)
    for (i in 1:B) {
      lt0 <- table(H.new[i, ] < 0)
      while (lt0[1] != ncol(prior) | 
               H.new[i, 1] < lo[1] | H.new[i, 1] > hi[1] | 
               H.new[i, 2] < lo[2] | H.new[i, 2] > hi[2] | 
               H.new[i, 3] < lo[3] | H.new[i, 3] > hi[3] | 
               H.new[i, 4] < lo[4] | H.new[i, 4] > hi[4] | 
               H.new[i, 5] < lo[5] | H.new[i, 5] > hi[5] | 
               H.new[i, 6] < lo[6] | H.new[i, 6] > hi[6] | 
               H.new[i, 7] < lo[7] | H.new[i, 7] > hi[7] | 
               H.new[i, 8] < lo[8] | H.new[i, 8] > hi[8] | 
               H.new[i, 9] < lo[9] | H.new[i, 9] > hi[9]) {
        H.new[i, ] <- mvrnorm(n = 1, mu = imp.theta, 
                              Sigma = Sigma.k)
        lt0 <- table(H.new[i, ] < 0)
      }
      if (i%%100 == 0) {
      }
    }
    H.new <- cbind(H.new[, 1:9])
    vwts <- c(vwts, var.rwts(wts.k))
    ewts <- c(ewts, entropy.wts(wts.k))
    mwts <- c(mwts, maxw.k)
    nup <- c(nup, expt.upts(wts.k, m = B))
    frac.nup <- c(frac.nup, expt.upts(wts.k, m = B)/B)
    N.k <- K * B
    case.n <- 1:N.k
    mwt.case <- c(mwt.case, case.n[wts.k == maxw.k])
    assign(paste("imp.theta.k", k, sep = ""), imp.theta)
    assign(paste("Sigma.k", k, sep = ""), Sigma.k)
    if (k%%5 == 0) {
      assign(paste("log.like.k", k, sep = ""), log.like.k)
      assign(paste("wts.k", k, sep = ""), wts.k)
      assign(paste("maxw.k", k, sep = ""), maxw.k)
      assign(paste("imp.theta", k, sep = ""), imp.theta)
      assign(paste("H.k", k, sep = ""), H.k)
      assign(paste("qk", k, sep = ""), qk)
      save.vectors <- c(paste("log.like.k", k, sep = ""), 
                        paste("wts.k", k, sep = ""), paste("maxw.k", 
                                                           k, sep = ""), paste("imp.theta", k, sep = ""), 
                        paste("Sigma.k", k, sep = ""), paste("H.k", k, 
                                                             sep = ""), "qk", "log.like")
      rm.vectors <- c(paste("log.like.k", k, sep = ""), 
                      paste("wts.k", k, sep = ""), paste("maxw.k", 
                                                         k, sep = ""))
      rm(list = rm.vectors)
    }
    if (frac.nup[k] >= 0.632) {
      break
    }
  }
  return(list(H.new = H.new, vwts = vwts, ewts = ewts, mwts = mwts, 
              nup = nup, frac.nup = frac.nup, wts.k = wts.k, mwt.case = mwt.case))
}