model {
for (i in 1:n) {
velocity[i] ~ dnorm(mu[i], tau[i])
mu[i] <- mu.mix[group[i]]
tau[i] <- tau.mix[group[i]]
group[i] ~ dcat(pi[])
for (j in 1:C) {
gind[i,j] <- equals(j, group[i])
}
}
p[1] <- q[1]
for (j in 2:C) {
p[j] <- q[j]*(1 - q[j-1])*p[j-1]/q[j-1]
}
for (j in 1:C) {
q[j] ~ dbeta(1, alpha)
pi[j] <- p[j]/sum(p[])
mu.mix[j] ~ dnorm(amu, mu.prec[j])
mu.prec[j] <- bmu*tau.mix[j]
tau.mix[j] ~ dgamma(aprec, bprec)
}
alpha <- 1
# Could replace constant alpha with a prior
# alpha ~ dgamma(2, 4)
# alpha ~ dunif(0.3, 10)
amu ~ dnorm(0, 0.001)
bmu ~ dgamma(0.5, 50)
aprec <- 2
bprec ~ dgamma(2, 1)
K <- sum(cl[])
for (j in 1:C) {
sumind[j] <- sum(gind[,j])
cl[j] <- step(sumind[j]-1+0.001) # cluster j used in this
# iteration
}
for (j in 1:ndens) {
for (i in 1:C) {
dens.cpt[i,j] <- pi[i]*
sqrt(tau.mix[i] / (2*3.141592654))*
exp(-0.5*tau.mix[i]*(mu.mix[i] - dens.x[j])
*(mu.mix[i] - dens.x[j]))
}
dens[j] <- sum(dens.cpt[,j])
}
}


Data:
list(velocity=c(9.172, 9.35, 9.483, 9.558, 9.775, 10.227, 10.406, 16.084, 16.17, 18.419, 18.552, 18.6, 18.927, 19.052, 19.07, 19.33, 19.343, 19.349, 19.44, 19.473, 19.529, 19.541, 19.547, 19.663, 19.846, 19.856, 19.863, 19.914, 19.918, 19.973, 19.989, 20.166, 20.175, 20.179, 20.196, 20.215, 20.221, 20.415, 20.629, 20.795, 20.821, 20.846, 20.875, 20.986, 21.137, 21.492, 21.701, 21.814, 21.921, 21.96, 22.185, 22.209, 22.242, 22.249, 22.314, 22.374, 22.495, 22.746, 22.747, 22.888, 22.914, 23.206, 23.241, 23.263, 23.484, 23.538, 23.542, 23.666, 23.706, 23.711, 24.129, 24.285, 24.289, 24.368, 24.717, 24.99, 25.633, 26.96, 26.995, 32.065, 32.789, 34.279), dens.x=c(8, 8.1, 8.2, 8.3, 8.4, 8.5, 8.6, 8.7, 8.8, 8.9, 9, 9.1, 9.2, 9.3, 9.4, 9.5, 9.6, 9.7, 9.8, 9.9, 10, 10.1, 10.2, 10.3, 10.4, 10.5, 10.6, 10.7, 10.8, 10.9, 11, 11.1, 11.2, 11.3, 11.4, 11.5, 11.6, 11.7, 11.8, 11.9, 12, 12.1, 12.2, 12.3, 12.4, 12.5, 12.6, 12.7, 12.8, 12.9, 13, 13.1, 13.2, 13.3, 13.4, 13.5, 13.6, 13.7, 13.8, 13.9, 14, 14.1, 14.2, 14.3, 14.4, 14.5, 14.6, 14.7, 14.8, 14.9, 15, 15.1, 15.2, 15.3, 15.4, 15.5, 15.6, 15.7, 15.8, 15.9, 16, 16.1, 16.2, 16.3, 16.4, 16.5, 16.6, 16.7, 16.8, 16.9, 17, 17.1, 17.2, 17.3, 17.4, 17.5, 17.6, 17.7, 17.8, 17.9, 18, 18.1, 18.2, 18.3, 18.4, 18.5, 18.6, 18.7, 18.8, 18.9, 19, 19.1, 19.2, 19.3, 19.4, 19.5, 19.6, 19.7, 19.8, 19.9, 20, 20.1, 20.2, 20.3, 20.4, 20.5, 20.6, 20.7, 20.8, 20.9, 21, 21.1, 21.2, 21.3, 21.4, 21.5, 21.6, 21.7, 21.8, 21.9, 22, 22.1, 22.2, 22.3, 22.4, 22.5, 22.6, 22.7, 22.8, 22.9, 23, 23.1, 23.2, 23.3, 23.4, 23.5, 23.6, 23.7, 23.8, 23.9, 24, 24.1, 24.2, 24.3, 24.4, 24.5, 24.6, 24.7, 24.8, 24.9, 25, 25.1, 25.2, 25.3, 25.4, 25.5, 25.6, 25.7, 25.8, 25.9, 26, 26.1, 26.2, 26.3, 26.4, 26.5, 26.6, 26.7, 26.8, 26.9, 27, 27.1, 27.2, 27.3, 27.4, 27.5, 27.6, 27.7, 27.8, 27.9, 28, 28.1, 28.2, 28.3, 28.4, 28.5, 28.6, 28.7, 28.8, 28.9, 29, 29.1, 29.2, 29.3, 29.4, 29.5, 29.6, 29.7, 29.8, 29.9, 30, 30.1, 30.2, 30.3, 30.4, 30.5, 30.6, 30.7, 30.8, 30.9, 31, 31.1, 31.2, 31.3, 31.4, 31.5, 31.6, 31.7, 31.8, 31.9, 32, 32.1, 32.2, 32.3, 32.4, 32.5, 32.6, 32.7, 32.8, 32.9, 33, 33.1, 33.2, 33.3, 33.4, 33.5, 33.6, 33.7, 33.8, 33.9, 34, 34.1, 34.2, 34.3, 34.4, 34.5, 34.6, 34.7, 34.8, 34.9, 35), ndens=271, C=20, n=82)


Inits:
list(amu=0, bmu=1)

Inits (if there is a prior on alpha):
list(amu=0, bmu=1, alpha=1)

Click for statistics with fixed alpha