[Normal-Gamma sampling]

normalgamma.R

library(assert)
source("multinorm.R")

rnormalgamma <- (function(k, b, Sigma, nu, kappa) {
  # Sample from a Normal-Gamma distribution Attaches both `rnormalgamma` and 
  # `rnormalgamma.fast` to the current environment.
  # 
  # Usage:
  #   rnormalgamma(b, Sigma, nu, kappa)
  #   rnormalgamma.fast(b, Sigma, nu, kappa)
  #   
  #   Returns pairs (beta, phi) such that 
  #     beta | phi  ~ N(b, Sigma/phi)
  #     phi         ~ Gamma(nu/2, kappa/2)
  #   where phi is a vector of sampled values and beta is a matrix where each 
  #   column represents a sample.
  
  rnormalgamma.fast <- function(k, b, Sigma, nu, kappa) {
    phi = rgamma(k, nu/2, kappa/2)
    beta = sapply(phi, function(phi) {
      rmultinorm.fast(1, b, Sigma/phi)
    })
    return(list(phi=phi, beta=beta))
  }
  rnormalgamma.fast(6, c(0,1), matrix(c(1,0,0,1), nrow=2), 1, 1)
  
  rnormalgamma <- function(k, b, Sigma, nu, kappa) {
    assert(is.numeric(k),
           is.numeric(b),
           is.matrix(as.matrix(Sigma)),
           is.numeric(nu),
           is.numeric(kappa),
           )
    assert(length(k) == 1,
           k > 0,
           k %% 1 == 0,
           length(b) == nrow(Sigma),
           length(b) == ncol(Sigma),
           length(nu) == 1,
           nu > 0,
           length(kappa) == 1,
           kappa > 0
    )
    return(rnormalgamma.fast(k, b, Sigma, nu, kappa))
  }
  
  assign("rnormalgamma.fast", rnormalgamma.fast, parent.env(environment()))
  return(rnormalgamma)
})()