R/sample_omega.R

In josbop/BayesBiSurv:

sample_omega <- function(omega, t, y, x, theta, gamma, partition, lambda) {

  bound.h <- cum_h(t, partition, lambda)

  len.out <- length(omega)
  omega.out <- vector(mode = "numeric", length = len.out)

  for (i in 1:len.out) {
    bound <- bound.h[i] * exp(x[[i]] %*% theta)
    ran.gam <- rgamma(n = 1, shape = 1 + y[i], scale = 1 + gamma)
    proposal <- ran.gam + bound

    # if the previous observation is not greater than the bound,
    # the Metropolis step always accepts

    if (omega[i] <= bound) {

      omega.out[i] <- proposal

    } else {

      log.alpha <- y[i] * (log(proposal) - log(ran.gam) - log(omega[i]) + log(omega[i] - bound))
      prob <- min(exp(log.alpha), 1)
      u <- runif(n = 1)

      omega.out[i] <- proposal * (u <= prob) + omega[i] * (u > prob)

    }


  }

  return(omega.out)

}