keebler = scan("cookies.txt")
chips = scan("chips.txt")
safeway = scan("safeway.txt")

y.sum = numeric(3)
y.sum[1] = sum(keebler) 
y.sum[2] = sum(chips) 
y.sum[3] = sum(safeway) 
n = numeric(3)
n[1] = length(keebler)
n[2] = length(chips)
n[3] = length(safeway)

N.iter = 5000
N.burnin = 100
alpha = numeric(N.iter+N.burnin)
lambda = matrix(0,ncol=3,nrow=N.iter+N.burnin)
alpha[1] = 250                           # start at prior means
lambda[1,] = c(25,25,25)
accept.alpha = 0
for(i in 1:(N.iter+N.burnin-1)) {
  alpha.star = rnorm(1,alpha[i],sd=15)   # may need to tune sd 
  if(alpha.star > 0) {
    prob = exp( 3*alpha.star*log(10)+624*log(alpha.star)-3*lgamma(alpha.star)-2.5*alpha.star + alpha.star*sum(log(lambda[i,])) - (3*alpha[i]*log(10)+624*log(alpha[i])-3*lgamma(alpha[i])-2.5*alpha[i] + alpha[i]*sum(log(lambda[i,]))) )
    if(runif(1) < prob) {
      alpha[i+1] = alpha.star
      accept.alpha = accept.alpha + 1
    }
    else alpha[i+1] = alpha[i]
  }
  else alpha[i+1] = alpha[i]
  #sample lambda.j ~ Gamma(alpha+y.sum.j, 10 + n.j)
  for(j in 1:3)   lambda[i+1,j] = rgamma(1,alpha[i+1]+y.sum[j], 10+n[j])
}
accept.alpha/(N.iter+N.burnin)
par(mfrow=c(2,2))
plot(alpha,pch=".")
for(j in 1:3) plot(lambda[,j],pch=".",ylab=paste("lambda[",j,"]"))

alpha = alpha[(N.burnin+1):(N.iter+N.burnin)]
lambda = lambda[(N.burnin+1):(N.iter+N.burnin),]
hist(alpha)
for(j in 1:3) hist(lambda[,j],pch=".",main=paste("Histogram of theta[",j,"]"))

mean(alpha)
var(alpha)
mean(lambda)
var(lambda)

apply(lambda,2,mean)   #separate means by column (for "by row", 2nd arg is 1)
y.sum/n                #compare to MLE's

acf(alpha)
for(j in 1:3) acf(lambda[,j],main=paste("Series theta[",j,"]"))