inst/sbc/lkj.R

In OncoBayes2: Bayesian Logistic Regression for Oncology Dose-Escalation Trials

## code from Ben Goodrich posted here
## https://groups.google.com/forum/#!msg/stan-users/3gDvAs_qwN8/Xpgi2rPlx68J
##

rgbeta <- function(num, shape) {
  if (shape == Inf) {
    rep(0, num)
  } else if (shape > 0) {
    -1 + 2 * rbeta(num, shape, shape)
  } else if (shape == 0) {
    -1 + 2 * rbinom(num, 1, 0.5)
  } else {
    stop("shape must be non-negative")
  }
}

## LKJ
rcorvine <- function(n, eta = 1, cholesky = FALSE, permute = !cholesky) {
  alpha <- eta + (n - 2) / 2
  L <- matrix(0, n, n)
  L[1, 1] <- 1
  L[-1, 1] <- partials <- rgbeta(n - 1, alpha)
  if (n == 2) {
    L[2, 2] <- sqrt(1 - L[2, 1]^2)
    if (cholesky) {
      return(L)
    }
    Sigma <- tcrossprod(L)
    if (permute) {
      ord <- sample(n)
      Sigma <- Sigma[ord, ord]
    }
    return(Sigma)
  }
  W <- log(1 - partials^2)
  for (i in 2:(n - 1)) {
    gap <- (i + 1):n
    gap1 <- i:(n - 1)
    alpha <- alpha - 0.5
    partials <- rgbeta(n - i, alpha)
    L[i, i] <- exp(0.5 * W[i - 1])
    L[gap, i] <- partials * exp(0.5 * W[gap1])
    W[gap1] <- W[gap1] + log(1 - partials^2)
  }
  L[n, n] <- exp(0.5 * W[n - 1])
  if (cholesky) {
    return(L)
  }
  Sigma <- tcrossprod(L)
  if (permute) {
    ord <- sample(n)
    Sigma <- Sigma[ord, ord]
  }
  return(Sigma)
}