# JAGS - problem with taking submatrices

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,

x[i,1:D] ~ dmnorm(mu[c[i],1:D], Omega[(D*(c[i]-1)+1):(D*c[i]),1:D])

JAGS supports 3 dimensional arrays for storing a list of matrices. Unlike WINBugs, JAGS uses column-major for storage, and like WINBugs prefers the entries to be in contiguous blocks. So, if you have $K$ groups with $\Omega_k$ the precision associated with the $k$'th group, then you would access this in JAGS via Omega[1:D, 1:D, k]. I've successfully used JAGS for mixtures of multivariate normals doing this, so it is possible. – guy Nov 24 '12 at 20:24
The matrices in your code aren't stored as contiguous blocks of memory. Suppose $\Omega$ is 100 x 100. Then Omega[1, 2] and Omega[2, 1] are very far apart, since it is stored as Omega[1, 1], Omega[2, 1], Omega[3, 1], ..., Omega[100, 1], Omega[1, 2], ... – guy Nov 25 '12 at 1:02