model {
for (i in 1:M) {
	y[i] ~ dnorm(mu[z[i]], tauinv[z[i]])
	z[i] ~ dcat(p[])
}

# stick-breaking prior
p[1] <- r[1]
for (j in 2:(N-1)) {p[j] <- r[j]*(1-r[j-1])*p[j-1]/r[j-1]}
for (k in 1:(N-1)) {r[k] ~ dbeta(1, alpha)}

#  sum p
ps <- sum(p[1:(N-1)])
p[N] <- 1 - ps

# baseline
for (k in 1:N) {
    mu[k] ~ dnorm(basemu, basetau)
    tauinv[k] ~ dgamma(3, b)
}
basemu ~ dnorm(0, .01)
basetau <- pow(sigmaF0, -2)
sigmaF0 ~ dunif(0, 10)
b ~ dgamma(0.03, 0.03)

# DPP parameter prior
alpha ~ dunif(0.3, 10)
# alpha  <- 1
}

library('rjags')
library('coda)
jags <- jags.model('example1.bug',
	data = list('x' = x,
	'N' = N),
	n.chains = 4,
	n.adapt = 100)

update(jags, 1000)

jags.samples(jags,
         c('mu', 'tau'),
         1000)

output <- coda.samples(model=jags,variable.names=c("basemu", "p","mu", "basetau"), n.iter=10000, thin=1)

print(summary(output))
plot(output)

data <- list(x = x, y = y)
init <- list(beta = 0, tau = 1)
jags <- jags.model(..., inits = init, data = data,
	n.chain = 3, n.adapt = 100)

par = c('alpha', 'beta')
samp = coda.sample(jags, par, n.iter = 20000)
plot(samp)
summary(samp)

burn.in = 1000
summary(window(samp, start = burn.in))

codamenu()

gelman.diag(samp)
gelman.plot(samp)

round(summary(window(samps[, paste("y[", 1:25, "]",
sep = "")], start = burn.in))$quantiles[, c(3, 1, 5)], 2)

y ~ dcat(p[]) # 1 ~ n

y[i, 1:3] ~ dmulti(q[], n[i]) # (1,2,n-3)
n[i] = sum(y[i,])

q[1:3] ~ ddirch(alpha[])
for (k in 1:3){alpha[k] <- 1}

T(,): a priori restriction
I(,): a posteriori restriction

Y[i] ~ dnormmix(mu[], tau[], p[])