## simulation from the multinomial-Dirichlet model
## ADM 10/10/2004 (with help from KQ)

library(MCMCpack) # provides Dirichlet density and generator 
G <- 10000
pi.post <- rdirichlet(G, c(278,82,35,10,1,26))

# posterior summary
cat("# posterior mean\n")
print(apply(pi.post, 2, mean))
cat("# posterior standard deviation\n")
print(apply(pi.post, 2, sd))
cat("# Pr(Dean > Other)\n")
prob.Dean.gt.Other <- mean(pi.post[,3] > pi.post[,6])
print(prob.Dean.gt.Other)

# posterior histograms
pdf("multinomial.pdf", height=6.5, width=6.5)
par(mfrow=c(3,2))
hist(pi.post[,1], xlim=c(0,.8), main="Kerry")
box()
hist(pi.post[,2], xlim=c(0,.8), main="Edwards")
box()
hist(pi.post[,3], xlim=c(0,.8), main="Dean")
box()
hist(pi.post[,4], xlim=c(0,.8), main="Kucinich")
box()
hist(pi.post[,5], xlim=c(0,.8), main="Sharpton")
box()
hist(pi.post[,6], xlim=c(0,.8), main="Other")
box()
dev.off()