## Julia program for Grouped Multinomial Data (Ex. 7.2.3)
# sample from Dirichlet distributions
function gmulti(T, x, a, b, alpha1 = 0.5, alpha2 = 0.5, alpha3 = 0.5)
z = ones(T+1, size(x, 1)-1) # initial z satisfy `z <= x`
# sample from g_1(theta | y)
dir = Dirichlet([z[t, 1] + z[t, 2] + alpha1, z[t, 3] + z[t, 4] + alpha2, x[5] + alpha3])
# sample from g_2(z | x, theta)
bi = Binomial(x[i], a[i]*mu[t+1]/(a[i]*mu[t+1]+b[i]))
z[t+1, i] = rand(bi, 1)[1]
bi = Binomial(x[i], a[i]*eta[t+1]/(a[i]*eta[t+1]+b[i]))
z[t+1, i] = rand(bi, 1)[1]
a = [0.06, 0.14, 0.11, 0.09];
b = [0.17, 0.24, 0.19, 0.20];