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 ~ N( 5, 1),  tau ~ exponential (1) = Gamma(2/2, 2/2), thus:
mu0 = 5
k = 1
alpha = 2
beta = 2

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

par(mfrow=c(2,1))
plot(mu)
plot(tau)

mu = mu[201:1200]
tau = tau[201:1200]
par(mfrow=c(1,2))
hist(mu)
hist(tau)

mean(mu)
sd(mu)

mean(tau)
1/sqrt(mean(tau))