bayes.factor.example.R

Jim — Nov 6, 2013, 8:12 AM

# computing bayes factors

# write a new function to faciliate the comparison of two models

bayes.factor <- function(model1, model2){
    require(LearnBayes)
    fit1 <- laplace(model1$logpost, model1$mode, model1$stuff)
    fit2 <- laplace(model2$logpost, model2$mode, model2$stuff)
    log.bayes.factor <- fit1$int - fit2$int
    print(paste('The log Bayes factor in support of model 1 over model 2 is',
          round(log.bayes.factor, 2)))
}

# example
# compare two models
# M1:  y is binomial(20, p), p is beta(3, 10)
# M2:  y is binomial(20, p), p is beta(10, 3)
# observe y = 5

# write a log posterior
# (don't forget normalizing constants)

lbinomial.beta <- function(p, stuff){
  y <- stuff$data[1]
  n <- stuff$data[2]
  a <- stuff$prior[1]
  b <- stuff$prior[2]
  dbinom(y, size=n, prob=p, log=TRUE) +
    dbeta(p, a, b, log=TRUE)
}

# use bayes.factor function

model1 <- list(logpost=lbinomial.beta, 
               mode=0.5, 
               stuff=list(data=c(5, 20), prior=c(3, 10)))

model2 <- list(logpost=lbinomial.beta, 
               mode=0.5, 
               stuff=list(data=c(5, 20), prior=c(10, 3)))

bayes.factor(model1, model2)

Loading required package: LearnBayes

[1] "The log Bayes factor in support of model 1 over model 2 is 4.61"

# here we can compare approximate Bayes factor with exact value

log.pred1 <- log(pbetap(c(3, 10), 20, 5)) -
             log(pbetap(c(10, 3), 20, 5))
round(log.pred1, 2)

[1] 4.61