Created
December 2, 2014 12:15
-
-
Save famuvie/842e9632873a16c45214 to your computer and use it in GitHub Desktop.
Plotting the marginal prior distributions for the elements of an Inverse-Wishart matrix
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
### Plotting the marginal prior distributions for the elements of | |
### an Inverse-Wishart matrix | |
library(MCMCpack) # Wishart distribution | |
library(actuar) # inverse gamma distribution | |
library(ggplot2) | |
## Hyperparameters of the inverse-Wishart prior | |
Phi <- matrix(c(1,.3,.3,2),2,2) # Scale matrix | |
nu <- 5 # Degrees of Freedom | |
## Marginal priors of the variances | |
## | |
## Use the fact that, if S ~ invW(Phi, nu), then | |
## s_ii ~ invGa(alpha = nu/2, beta = Phi_ii/2) | |
limits <- sapply(diag(Phi), | |
function(x) qinvgamma(c(.01, .99), shape = nu/2, scale = x/2)) | |
## We can choose wether to plot the marginal prior variances under | |
## the same scale or not. If that's the case, we use the maximal interval | |
limits <- c(min(limits[1, ]), max(limits[2, ])) | |
x <- seq(limits[1], limits[2], length = 501) | |
y <- do.call('c', | |
lapply(diag(Phi), | |
function(ss) dinvgamma(x, | |
shape = nu/2, | |
scale = ss/2))) | |
priordf.var <- data.frame(comp = factor(c(rep('sigma2_a', length(x)), | |
rep('sigma2_c', length(x)))), | |
x = c(x, x), | |
y) | |
## Plot prior curves and true value | |
ggplot(priordf.var, aes(x, y, col = comp)) + | |
geom_line() + | |
theme_bw() | |
## For the marginal correlation, the analytic derivation of the distribution | |
## is not very straightforward, so we just sample and plot an empirical density | |
## Use the fact that A ~ invW(Phi, nu) iff invA ~ W(invPhi, nu) | |
## and build a function to get a sample of the correlation or covariance element | |
rinvWish <- function(n, Phi, nu, corr = TRUE) { | |
invsample <- lapply(seq.int(n), | |
function(x) rwish(nu, solve(Phi))) | |
sample <- lapply(invsample, solve) | |
cov <- sapply(sample, | |
function(x) x[1,2]) | |
if( corr ) { | |
denom <- apply(sqrt(sapply(sample, diag)), 2, prod) | |
stopifnot(abs(cov) < denom) | |
cov <- cov/denom | |
} | |
return(cov) | |
} | |
N <- 1e4 | |
cor.smpl <- rinvWish(N, Phi, nu, corr = TRUE) | |
ggplot(as.data.frame(density(cor.smpl, from = -1, to = 1)[1:2]), | |
aes(x, y)) + | |
geom_line() |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment