#
#Simple routines for conjugate prior analysis of multinomial likelihood
#problems, for the R statistical computing system.
#
#by Carter Butts <ctb@andrew.cmu.edu>
#
#v0.0.1, last updated August 23, 2000
#
#Contents:
#
#bmulti.compare - Compare two multinomial models
#bmulti.post - Fit a Bayesian multinomial model with conjugate priors
#plot.bayes.multi.post - Plot a multinomial model object
#print.bayes.multi.post - Print a multinomial model object
#print.bmulti.compare - Print a multinomial model comparison object
#rdirichlet - Draw from the Dirichlet distribution
#summary.bayes.multi.post - Summarize a multinomial model object
#
#NOTES:
#
#License: This software is made available under the terms of the GNU
#Public License (GPL), version 2.0 or later.  Please visit the Free
#Software Foundation (www.fsf.org) for more details.
#
#Lack of support: It must be noted that this code is provided AS IS,
#without support or warranty (express or implied).  The author will
#generally try to fix bugs were possible, but the end user is
#generally on his/her own with respect to this software; the author
#apologizes for this fact, but would like to point out that he has
#numerous professional and other obligations which prevent his
#working full-time on this project.
#
#Implementation: This code has been written for the R implementation of
#the S language, and is currently thought to work with version 1.0 and
#higher.  Systematic testing, however, has been limited: caveat emptor. 
#

rdirichlet<-function(ndraw,params){
   n<-length(params)
   draws<-matrix(ncol=n,nrow=ndraw)
   for(i in 1:ndraw){
      draws[i,]<-rgamma(n,params)
      draws[i,]<-draws[i,]/sum(draws[i,])
   }
   draws
}

bmulti.post<-function(counts,priors=rep(1,length(counts)),catnames=c(1:length(counts)),draws=10000){
   out<-data.frame()
   n<-length(priors) 
   out$draws<-rdirichlet(draws,priors+counts)
   out$catnames<-catnames
   out$parameters<-priors+counts
   class(out)<-"bayes.multi.post"
   out
}

print.bayes.multi.post<-function(bmp){
   n<-dim(bmp$draws)[2] 
   cat("\nConjugate Bayesian Analysis - Multinomial Likelihood\n\n")
   cat("Joint Posterior Distribution:\n\n")
   cat("p(theta|y) ~ Dirichlet w/ parameters\n\n")
   pparam<-matrix(nrow=1,ncol=n)
   pparam[1,]<-bmp$parameters
   rownames(pparam)<-""
   colnames(pparam)<-bmp$catnames
   print(pparam)
   cat("\n")

   ptab<-matrix(nrow=8,ncol=n)
   rownames(ptab)<-c("Min","5%","25%","50%","75%","95%","Max","Mean")
   colnames(ptab)<-bmp$catnames
   for(i in 1:n){
      ptab[1:7,i]<-quantile(bmp$draws[,i],c(0,0.05,0.25,0.5,0.75,0.95,1))
      ptab[8,i]<-mean(bmp$draws[,i])   
   }
   cat("Marginal Posterior Parameter Distributions:\n\n")
   print(ptab)
   cat("\n")
}

summary.bayes.multi.post<-function(bmp){
   n<-dim(bmp$draws)[2] 
   cat("\nConjugate Bayesian Analysis - Multinomial Likelihood\n\n")
   cat("Joint Posterior Distribution:\n\n")
   cat("p(theta|y) ~ Dirichlet w/ parameters\n\n")
   pparam<-matrix(nrow=1,ncol=n)
   pparam[1,]<-bmp$parameters
   rownames(pparam)<-""
   colnames(pparam)<-bmp$catnames
   print(pparam)
   cat("\n")

   ptab<-matrix(nrow=8,ncol=n)
   rownames(ptab)<-c("Min","5%","25%","50%","75%","95%","Max","Mean")
   colnames(ptab)<-bmp$catnames
   for(i in 1:n){
      ptab[1:7,i]<-quantile(bmp$draws[,i],c(0,0.05,0.25,0.5,0.75,0.95,1))
      ptab[8,i]<-mean(bmp$draws[,i])   
   }
   cat("Marginal Posterior Parameter Distributions:\n\n")
   print(ptab)
   cat("\n")

   pagrb<-matrix(nrow=n,ncol=n)
   rownames(pagrb)<-bmp$catnames
   colnames(pagrb)<-bmp$catnames
   for(i in 1:n)
      for(j in 1:n)
         pagrb[i,j]<-mean(as.numeric(bmp$draws[,i]>bmp$draws[,j]))
   cat("Marginal Posterior Probability of Rowparam>Colparam:\n\n")
   print(pagrb)
   cat("\n")
}

plot.bayes.multi.post<-function(bmp,mode="density",...){
   n<-dim(bmp$draws)[2]
   par(mfrow=c(ceiling(sqrt(n)),ceiling(sqrt(n)))) 
   if(mode=="density"){
      for(i in 1:n){
         plot(density(bmp$draws[,i]),main=paste("Marginal Posterior of",bmp$catnames[i]),xlim=c(0,1),...)
         abline(v=median(bmp$draws[,i]),lty=2)
         abline(v=quantile(bmp$draws[,i],0.05),lty=3)
         abline(v=quantile(bmp$draws[,i],0.95),lty=3)
      }
      #mtext(side=3,line=-2,outer=TRUE,"Densities of Marginal Posterior Draws")
   }else{
      for(i in 1:n){
         hist(bmp$draws[,i],main=paste("Marginal Posterior of",bmp$catnames[i]),probability=TRUE,freq=FALSE,xlim=c(0,1),xlab=bmp$catnames[i],...)
         abline(v=median(bmp$draws[,i]),lty=2)
         abline(v=quantile(bmp$draws[,i],0.05),lty=3)
         abline(v=quantile(bmp$draws[,i],0.95),lty=3)
      }
      #mtext(side=3,line=-2,outer=TRUE,"Histograms of Marginal Posterior Draws")
   }
}

bmulti.compare<-function(bmp1,bmp2,lab){
   out<-data.frame()
   n<-dim(bmp1$draws)[2]
   m<-matrix(nrow=2,ncol=n)
   for(i in 1:n){
      m[1,i]<-mean(as.numeric(bmp1$draws[,i]>=bmp2$draws[,i]))
      m[2,i]<-mean(as.numeric(bmp1$draws[,i]<=bmp2$draws[,i]))
   }
   out$probs<-m
   out$catnames<-bmp1$catnames
   out$setnames<-lab
   class(out)<-"bmulti.compare"
   out
}

print.bmulti.compare<-function(bmc){
   cat("Posterior Comparison of ",bmc$setnames[1],"and",bmc$setnames[2],"\n\n")
   m<-bmc$probs
   rownames(m)<-c(paste("p(",bmc$setnames[1],">=",bmc$setnames[2],")"),paste("p(",bmc$setnames[1],"<=",bmc$setnames[2],")"))
   colnames(m)<-bmc$catnames
   print(m)
   cat("\n")
}