y.bar = 5.494  # y_i ~ N(mu, .106^2)
n = 40
#priors: mu ~ N(5.5, 1/tau), tau ~ exp(1)

N = 5100
mu = numeric(N)
tau = numeric(N)
mu[1] = 5.5
tau[1] = 1
for(i in 1:(N-1)) {
  mu[i+1] = rnorm(1, (n*y.bar/(.106^2)+5.5*tau[i])/(n/(.106^2)+tau[i]),
      sd=1/sqrt(n/(.106^2)+tau[i]))
  tau[i+1] = rgamma(1, 1.5, 0.5*(mu[i+1]-5.5)^2+1)
}
par(mfrow=c(2,1))
plot(mu,pch=".")
plot(tau,pch=".")

mu = mu[101:5100]
tau = tau[101:5100]
par(mfrow=c(1,2))
hist(mu)
hist(tau)
lines(seq(from=0.5,to=9.5,by=1),5000*dexp(seq(from=0.5,to=9.5,by=1)))

par(mfrow=c(1,1))
qqplot(sort(tau),qexp(1:5000/5001))

mean(mu); sd(mu); mean(tau); sd(tau)

par(mfrow=c(1,2))
acf(mu); acf(tau)


#convergence diagnostics -- back to Metropolis-Hastings
accept.mu = 0
accept.tau = 0
for(i in 1:(N-1)) {
  #sample mu using a normal proposal
  mu.star = rnorm(1, mu[i])
  alpha = exp(-0.5 * (n/(.106^2)+tau[i]) *
    ( (mu.star - (n*y.bar/(.106^2)+5.5*tau[i])/(n/(.106^2)+tau[i]))^2 -
        (mu[i] - (n*y.bar/(.106^2)+5.5*tau[i])/(n/(.106^2)+tau[i]))^2 ) )
  if(runif(1) < alpha) {
    mu[i+1] = mu.star
    accept.mu = accept.mu + 1
  }
  else mu[i+1] = mu[i]
  #sample tau using a normal proposal
  tau.star = rnorm(1, tau[i])
  if(tau.star > 0) {
    alpha = sqrt(tau.star/tau[i])*
      exp( (tau[i]-tau.star) * (0.5*(mu[i+1]-5.5)^2 + 1) )
    if(runif(1) < alpha) {
      tau[i+1] = tau.star
      accept.tau = accept.tau + 1
    }
    else tau[i+1] = tau[i]
  }
  else tau[i+1] = tau[i]
}
par(mfrow=c(2,1))  # note - no burn-in
plot(mu,pch=".")
plot(tau,pch=".")
accept.mu/N
accept.tau/N

acf(mu)
acf(tau)

#adjust both proposals
  mu.star = rnorm(1, mu[i], sd=0.1)
  tau.star = rnorm(1, tau[i], sd=2)

#overcompensate
  mu.star = rnorm(1, mu[i], sd=0.001)
  tau.star = rnorm(1, tau[i], sd=10)

#tau steps too small
  tau.star = rnorm(1, tau[i], sd=0.1)