#comparing models for flipping a coin
n=10
y=seq(0,10,by=1)
plot(y, 1/(1+ 2^n*gamma(y+1)*gamma(n-y+1)/gamma(n+2)) )
abline(h=0.5)

n=100
y=seq(0,100,by=1)
plot(y, 1/(1+ 2^n*gamma(y+1)*gamma(n-y+1)/gamma(n+2)) )
abline(h=0.5)

#comparing regression models with Bayes factors and BIC
beer = read.table("beer.txt")
names(beer) = c("brand","price","calories","alcohol","type","domestic")

library(MCMCpack)
model1 = MCMCregress(price ~ domestic + calories + alcohol, data=beer, B0=.1,
  marginal.likelihood="Laplace")
summary(model1)

model2 = MCMCregress(price ~ domestic + calories, data=beer, B0=.1,
  marginal.likelihood="Laplace")
summary(model2)

model3 = MCMCregress(price ~ domestic, data=beer, B0=.1,
  marginal.likelihood="Laplace")
summary(model3)

model4 = MCMCregress(price ~  calories, data=beer, B0=.1,
  marginal.likelihood="Laplace")
summary(model4)

bf = BayesFactor(model1, model2, model3, model4)
summary(bf)

model1.mle = lm(price ~ domestic + calories + alcohol, data=beer)
summary(model1.mle)
model2.mle = lm(price ~ domestic + calories, data=beer)
summary(model2.mle)
model3.mle = lm(price ~ domestic, data=beer)
summary(model3.mle)
model4.mle = lm(price ~ calories, data=beer)
summary(model4.mle)

library(nlme)   # contains BIC function
BIC(model1.mle) # this is -2*BIC from class notes
BIC(model2.mle) # so best model minimizes this BIC
BIC(model3.mle)
BIC(model4.mle)