# METROPOLIS ALGORITHM FOR VARIANCE COMPONENTS BAYESIAN MODEL
#   WITH "DYESTUFFS" DATA


# LOG OF THE VARIANCE COMPONENTS POSTERIOR DENSITY FUNCTION
logvarcompdens = function(V, W, mu, theta, Y) {
  # theta = vector of length K, and Y = vector of length K*J
  if ( (V <= 0) || (W <= 0) ) {
    return(-Inf);
  } else {
    KK = length(theta);
    JJ = length(Y) / length(theta);
    thetamat = t( matrix( rep(theta,JJ), nrow=JJ ) );
    return(
      -b1/V + log(V)*(-a1-1) + (-b2/W) + log(W)*(-a2-1) -
      (mu-a3)^2 / (2*b3) + log(V)*(-KK/2) + log(W)*(-JJ*KK/2) -
      sum( (theta - mu)^2 ) / (2*V) - sum( (Y - thetamat)^2 ) / (2*W)
    );
  }
}

# SAME THING, BUT DEFINED AS A FUNCTION OF A SINGLE VECTOR
thelogdens = function(argvec) {
  logvarcompdens(argvec[1], argvec[2], argvec[3], argvec[4:(K+3)], Ydye)
}

# MATRIX giving the famous "dyestuff" batch data from Davies (1967).
# Defined so Ydye[i,j] equals yield (in grams) from j'th sample of i'th batch.
# Valid for i=1,2,3,4,5,6, and j=1,2,3,4,5, i.e. K=6 and J=5.
Ydye = t( matrix(
	c(1545, 1440, 1440, 1520, 1580,
         1540, 1555, 1490, 1560, 1495,
         1595, 1550, 1605, 1510, 1560,
         1445, 1440, 1595, 1465, 1545,
         1595, 1630, 1515, 1635, 1625,
         1520, 1455, 1450, 1480, 1445), nrow=5) )
K=6; J=5

# PRIOR DISTRIBUTION CONSTANTS
a1 = a2 = b1 = b2 = 100;
a3 = 1500;
b3 = (200)^2;
# b1 = 5

# MCMC RUN VARIABLES
M = 50000
B = 20000
printlength = 100
sigma = 0.5
numacc = 0
Vlist = Wlist = mulist = theta1list = theta2list = biglist = rep(0,M)

# INITIAL DISTRIBUTION:
currvec = c(rnorm(2,100,100), rnorm(K+1,1500,200))
for (i in 1:length(currvec))
  currvec[i] = max(1.0, currvec[i]);

# RUN THE METROPOLIS ALGORITHM
for (i in 1:M) {
  if (i == printlength*floor(i/printlength))
    cat(i, "/", M, " ... ", sep="")
  newvec = currvec + sigma * rnorm(K+3)
  U = runif(1)
  if (log(U) < thelogdens(newvec) - thelogdens(currvec)) {
    # accept the proposal (newvec)
    currvec = newvec
    numacc = numacc + 1
  }
  Vlist[i] = currvec[1]
  Wlist[i] = currvec[2]
  mulist[i] = currvec[3]
  theta1list[i] = currvec[4]
  theta2list[i] = currvec[5]
  if (Vlist[i] < Wlist[i])
    biglist[i] = 1;
}

cat("\n\nacceptance rate =", numacc/M, "\n")
cat("mean of V =", mean(Vlist[(B+1):M]), "\n")
cat("mean of W =", mean(Wlist[(B+1):M]), "\n")
cat("mean of W/V =", mean(Wlist[(B+1):M]/Vlist[(B+1):M]), "\n")
cat("mean of mu =", mean(mulist[(B+1):M]), "\n")
cat("mean of theta1 =", mean(theta1list[(B+1):M]), "\n")
cat("mean of theta2 =", mean(theta2list[(B+1):M]), "\n")
# cat("P(V<W)  =", mean(biglist[(B+1):M]), "\n\n")

plot(Wlist, type='l')
# plot(Vlist, type='l')
# plot(Vlist[(B+1):M],type='l')
# plot(mulist, type='l')
# plot(theta1list, type='l')
# plot(theta2list, type='l')