Sampling from a bivariate Normal-Inverse Gamma Distribution

We load dthe MCMCpack, which we need to sample froma an inverse gamma univariate distribution which we will need to do

suppressPackageStartupMessages(library(MCMCpack))
suppressPackageStartupMessages(library(ggplot2))

To sample a bivariate Normal-Inverse Gamma \( (X,Y) \sim NIG(\mu, \sigma, shape, scale) \) we follow three steps.

1. Set all the parameters of the distribution and the size of the sample

n <- 1000  #size of the sample
mu <- 3
sigma <- 5
shape <- 4
scale <- 6

2. Sample from \( Y \sim IG(shape, scale) \)

Y <- rinvgamma(n, shape, scale)

3. Sample from a \( X|Y \sim N(\mu, Y\sigma^2) \)

X <- rnorm(n, mu, sd = sqrt(Y * sigma^2))

Lets graph the bivariate sample we have just drawn

smpl.NIG <- data.frame(N = X, IG = Y)

We graph \( X \sim N(\mu, Y\sigma^2) \)

hist(smpl.NIG$N, 100, prob = TRUE)
lines(x = sort(smpl.NIG$N), y = dnorm(sort(smpl.NIG$N), mu, sigma), col = "red")

We graph \( Y \sim IG(shape, scale) \)

hist(smpl.NIG$IG, 100, prob = TRUE)
lines(x = sort(smpl.NIG$IG), y = dinvgamma(sort(smpl.NIG$IG), shape, scale), 
    col = "red")

We graph the bivariate \( (X,Y) \sim NIG(\mu, \sigma, shape,scale) \)

p <- ggplot(data = smpl.NIG)
p <- p + geom_jitter(aes(x = N, y = IG), alpha = 0.9, colour = "grey")
p <- p + geom_density2d(aes(x = N, y = IG))
p