# Bayesian fitting for multinomial-dirichlet model.
graphics.off()
rm(list=ls(all=TRUE))
#install.packages("BRugs")
library(BRugs)

#------------------------------------------------------------------------------
# THE MODEL.
modelstring = "
# BUGS model specification begins here...
model {

    # needs baserate loops added at some point
    # the data will be in list format, so the model will need
    # to group by stimulus, subj, and baserate

    # Likelihood function:
    for( sb in 1 : Nsubj ) {
      for( st in 1 : Nstim ) {
        data_table[ st + Nstim*(sb-1), 1:Nresp]
          ~ dmulti( theta[st + Nstim*(sb-1), 1:Nresp], Ntr[st + Nstim*(sb-1)] )
      }
    }

    # Prior:
    for( sb in 1 : Nsubj ) {
      for( st in 1 : Nstim ) {
        theta[st + Nstim*(sb-1), 1:Nresp] ~ ddirch( alpha[st,1:Nresp] )
      }
    }

    # hyperprior:
    for( st in 1 : Nstim ) {
      for( re in 1 : Nresp ) {
        #alpha[st,re] ~ dunif(0,250)
        alpha[st,re] ~ dgamma( 3, .02)
      }
    }

}
# ... end BUGS model specification
" # close quote for modelstring
# Write model to a file:
.temp = file("model.txt","w") ; writeLines(modelstring,con=.temp) ; close(.temp)
# Load model file into BRugs and check its syntax:
modelCheck( "model.txt" )

#------------------------------------------------------------------------------
# THE DATA.

modelData(fileName="mcg_data.txt")

#------------------------------------------------------------------------------
# INTIALIZE THE CHAINS.

nchain = 3
modelCompile( numChains = nchain )
modelGenInits()

                                              
#------------------------------------------------------------------------------
# RUN THE CHAINS

# burn in
BurnInSteps = 10000 # try 20000 or so
modelUpdate( BurnInSteps )
# actual samples
samplesSet( c( "theta" , "alpha" ) )
stepsPerChain = 1100
thinStep = 200 # try 100 or so
modelUpdate( stepsPerChain , thin=thinStep )