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 ~ chi-sq( 2, 1) = exponential (1), thus:
mu0 = 5
kappa = 1
v = 2
w = 1

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 + kappa*mu0)/(n + kappa) ),
       sd=1/sqrt( (n + kappa)*tau[i]) )
  ysum = 0
  for(j in 1:n) ysum = ysum + (alcohol[j] - mu[i+1])^2
  tau[i+1] = rgamma(1, (n + v + 1)/2, (ysum + kappa*(mu[i+1]-mu0) + v/w)/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))