###############################################
#Multinomial Data Example
#Pre-Election Poll (see pg. 83 in the text)
#1/19/06
###############################################

#Survey Data
y1 <- 727            #number of Bush supporters
y2 <- 583            #number of Dukakis supporters
y3 <- 137            #number who support other candidates
n.individs <- y1+y2+y3

#Parameters of the posterior distribution of theta
alpha.post <- c(y1,y2,y3)

#Function to sample from a Dirichlet distribution -- see algorithm on pg. 582 of the textbook
rdirichlet <- function(n=1,alpha){
                k <- length(alpha)
                x <- rep(0,k)
                theta <- matrix(0,nrow=n,ncol=k)
                for(i in 1:n){
                    x <- apply(as.matrix(alpha),1,function(shape.x) return(rgamma(1,shape=shape.x,scale=1)))
                    theta[i,] <- x/sum(x)
                    }
                return(theta)
                }

#Draw a random sample from the posterior
rdirichlet(1,alpha.post)

#Draw 1000 random samples from the posterior
post.samples <- rdirichlet(1000,alpha.post)

#Histograms of the posterior samples of the theta's
hist.labels <- c(expression(theta[1]),expression(theta[2]),expression(theta[3]))
hist.titles <- c("Posterior Samples","","")
par(mfrow=c(3,1))
for(j in 1:3){
    hist(post.samples[,j],main=hist.titles[j],xlab=hist.labels[j],ylab="",xlim=c(0,1))
    }

#Histogram of samples from the posterior distribution of theta_1-theta_2 -- INTERPRETATION?
par(mfrow=c(1,1))
hist(post.samples[,1]-post.samples[,2],main="Posterior Samples",xlab=expression(theta[1]-theta[2]),ylab="")