# Foundations Lab, Day 3: Some Bayesian Packages
# by Stephen R. Haptonstahl
# EITM 2008, Washington University in St Louis

# Check out the Bayesian Task View on CRAN for a complete
# list of R packages for Bayesian Inference:
# http://cran.wustl.edu/web/views/Bayesian.html

###  Create a dataset  ###
n <- 100
mydata <- data.frame(x1=rnorm(n, sd=3))
mydata$y <- mydata$x1 / 3 + .2 + rnorm(n, sd=.2)
mydata$z <- pnorm(mydata$y)
mydata$y2 <- sapply(1:n, function(i) 
  sample(0:1, 1, replace=T, prob=c(1-mydata$z[i],mydata$z[i]))
)

par(mfrow=c(1,2))
plot(y ~ x1, data=mydata)
plot(y2 ~ x1, data=mydata, pch="|")
points(z ~ x1, data=mydata)

###  Frequentist analyses  ###

lm.freq <- lm(y ~ x1, data=mydata)
summary(lm.freq)

probit.freq <- glm(y2 ~ x1, data=mydata, family=gaussian())
summary(probit.freq)

###  MCMCpack  ###
library(MCMCpack)

lm.posterior <- MCMCregress(y ~ x1, data=mydata)
summary(lm.posterior)
plot(lm.posterior)

probit.posterior <- MCMCprobit(y2 ~ x1, data=mydata)
summary(probit.posterior)
plot(probit.posterior)

# compare to frequentist results
summary(lm.freq)$coefficients[,1:2]; summary(lm.posterior)$statistics[1:2,1:2]