I'm trying to use JAGS for clustering mixtures of multivariate normal distributions. In order to model different covariance structures in each cluster I wanted to use one big (C*D)xD matrix (C - number of clusters, D - number of dimensions) which would represent all of the covariance matrices stacked one upon another. The problem is, that while JAGS is perfectly happy with sampling such submatrices from Wishart distribution, it will not allow me to subsequently use them as parameters for the normal distribution (or at least so I suspect). Here is my (not so) minimal working example:
library(rjags)
library(MASS)
clusters.bug <-
"model {
# N - observations
# D - dimensions
# C - clusters
for (i in 1:N) {
# JAGS has some problem with the following indexing here:
x[i,1:D] ~ dmnorm(mu[c[i],1:D], Omega[(D*(c[i]-1)+1):(D*c[i]),1:D])
# x[i,1:D] ~ dmnorm(mu[c[i],1:D], Omega0)
c[i] ~ dcat(p)
}
for (j in 1:C) {
mu[j,1:D] ~ dmnorm(mu0, Omega0)
Omega[(D*(j-1)+1):(D*j),1:D] ~ dwish(Omega0, D)
}
for (k in 1:D) {
mu0[k] <- 0
for (l in 1:D) {
Omega0[k,l] <- ifelse(k == l, 1, 0)
}
}
p ~ ddirch(alpha)
}"
## Data
N <- 1000
D <- 2
C <- 3
p <- c(0.2, 0.3, 0.5)
mu <- matrix(c(1,2, 3,4, 5,6), ncol = 2, byrow = TRUE)
cv <- matrix(c(1,0, 0,1, 1,.5, .5,2, 2,.3, .3,1), ncol = 2, byrow = TRUE) / 10
x <- NULL
for(i in 1:N) {
k <- sample(1:C, 1, prob = p)
x <- rbind(x, mvrnorm(1, mu[k,], cv[(D*(k-1)+1):(D*k),]))
}
contour(kde2d(x[,1], x[,2]))
alpha <- rep(1, C)
Omega <- NULL
for(i in 1:C) {
Omega <- rbind(Omega, diag(D))
}
## JAGS
jags <- jags.model(textConnection(clusters.bug),
data = list('x' = x,
'N' = N,
'D' = D,
'C' = C,
'alpha' = alpha),
n.chains = 1,
n.adapt = 100)
If you comment out:
x[i,1:D] ~ dmnorm(mu[c[i],1:D], Omega[(D*(c[i]-1)+1):(D*c[i]),1:D])
and uncomment the next line, everything is ok.