beer = read.table("beer.txt")
names(beer) = c("brand","price","calories","alcohol","type","domestic")
  # type: 1=craft lager, 2=craft ale, 3=import lager,
  #    4=regular & ice beer, 5=light/no alcohol beer
  # domestic: 1=U.S., 2=import
beer = beer[1:78,]  # remove non-alcoholic beers
attach(beer)
hist(alcohol)

n = length(alcohol)
ybar = mean(alcohol)

#priors: mu|sigma2 ~ N( 5, sigma2),  sigma2 ~ Inverse-Gamma(1, 1), thus:
m = 5
s0 = 1
alpha = 1
beta = 1

N = 1200
mu = numeric(N)
sigma2 = numeric(N)
mu[1] = 5
sigma2[1] = 1
for(i in 1:(N-1)) {
  mu[i+1] = rnorm(1, mean=( (n*ybar + s0*m)/(n + s0) ),
       sd=sqrt( sigma2[i]/(n + s0) ) )
  ysum = 0
  for(j in 1:n) ysum = ysum + (alcohol[j] - mu[i+1])^2
  sigma2[i+1] = 1/rgamma(1, (n+1)/2+alpha, 0.5*ysum + 0.5*s0*(mu[i+1]-m)^2 + beta)
}

par(mfrow=c(2,1))
plot(mu,pch=".")
plot(sigma2,pch=".")

mu = mu[201:1200]
sigma2 = sigma2[201:1200]
acf(sigma2)
acf(sigma2)

par(mfrow=c(1,2))
hist(mu)
hist(sigma2)

mean(mu)
sd(mu)

mean(sigma2)