Jim — Nov 6, 2013, 7:44 PM
# Two Bayes Factor Examples
###############################################
# Link Barker Example
# Described in Albert Notes
###############################################
# Poisson log posterior
# y is vector of observations
# prior is uniform on [0, T]
lb.lpost1 <- function(lambda, datapar){
y <- datapar$y
T <- datapar$T
sum(dpois(y, lambda, log=TRUE)) - log(T)
}
# geometric log posterior
# y is vector of observations
# prior is 1/T 1/p^2
lb.lpost2 <- function(p, datapar){
y <- datapar$y
T <- datapar$T
sum(dgeom(y, p, log=TRUE)) - log(T) - 2 * log(p)
}
# enter in data
y <- c(0, 1, 2, 3, 8)
# the new bayes.factor function facilitates 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',
model1$name, 'over', model2$name, 'is',
round(log.bayes.factor, 2)))
}
# for each model, create a list
# - name: name of model
# - logpost: name of function containg the log posterior
# - mode: rough guess at mode of the posterior
# - stuff: list of the inputs needed in the function
model1 <- list(name="Poisson Model",
logpost=lb.lpost1,
mode=2,
stuff=list(T=20, y=y))
model2 <- list(name="Geometric Model",
logpost=lb.lpost2,
mode=0.2,
stuff=list(T=20, y=y))
# compare the geometric and Poisson models -- output is log Bayes factor
bayes.factor(model2, model1)
Loading required package: LearnBayes
[1] "The log Bayes factor in support of Geometric Model over Poisson Model is 2.66"
# to get the bayes factor, take the exponential of the value
exp(2.66)
[1] 14.3
############################################
# Compare Two Binomial Proportions
############################################
# yj distributed binomial(nj, pj), j=1, 2
# want to test p1 = p2
# Model 1: two proportions, each given a uniform prior
# Model 2: p1=p2, common proportion given a uniform prior
# define the two log posteriors
log.2p <- function(theta, stuff){
p1 <- theta[1]
p2 <- theta[2]
y1 <- stuff[1]
n1 <- stuff[2]
y2 <- stuff[3]
n2 <- stuff[4]
dbinom(y1, size=n1, prob=p1, log=TRUE) +
dbinom(y2, size=n2, prob=p2, log=TRUE)
}
log.1p <- function(p, stuff){
y1 <- stuff[1]
n1 <- stuff[2]
y2 <- stuff[3]
n2 <- stuff[4]
dbinom(y1, size=n1, prob=p, log=TRUE) +
dbinom(y2, size=n2, prob=p, log=TRUE)
}
# data
y1 <- 4; n1 <- 30; y2 <- 12; n2 <- 30
# define each model
model1 <- list(name="One Proportion Model",
logpost=log.1p,
mode=0.5,
stuff=c(y1, n1, y2, n2))
model2 <- list(name="Two Proportions Model",
logpost=log.2p,
mode=c(0.5, 0.5),
stuff=c(y1, n1, y2, n2))
# compute the log Bayes factor in support of model 2
bayes.factor(model2, model1)
[1] "The log Bayes factor in support of Two Proportions Model over One Proportion Model is 1.41"
# contrast with the usual frequentist test of p1 = p2
prop.test(c(y1, y2), c(n1, n2))
2-sample test for equality of proportions with continuity
correction
data: c(y1, y2) out of c(n1, n2)
X-squared = 4.176, df = 1, p-value = 0.041
alternative hypothesis: two.sided
95 percent confidence interval:
-0.51337 -0.01996
sample estimates:
prop 1 prop 2
0.1333 0.4000
#######################################################