Take the 2-minute tour ×
Cross Validated is a question and answer site for people interested in statistics, machine learning, data analysis, data mining, and data visualization. It's 100% free, no registration required.

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.

share|improve this question
1  
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
    
Thanks for your answer - the 3D indexing solved my problem. What bothers me a little is that according to what you said (JAGS preference of continuity) it should be doable via switching indexes in the Omega matrix (storing subsequent covariance matrices in subsequent D columns instead of rows), however this does not work. Nonetheless your solution is much more elegant. –  kman Nov 24 '12 at 21:46
    
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
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Browse other questions tagged or ask your own question.