varcompdens = function(V, W, mu, theta, Y) {
  # theta = vector of length K, and Y = vector of length K*J
  if ( (V <= 0) || (W <= 0) ) {
    return(0)
  } else {
    KK = length(theta)
    JJ = length(Y) / length(theta)
    thetarep = rep(theta,JJ)
    return(
      exp(-b1/V) * V^(-a1-1) * exp(-b2/W) * W^(-a2-1) *
      exp(-(mu-a3)^2 / (2*b3)) * V^(-KK/2) * W^(-JJ*KK/2) *
      exp( - sum( (theta - mu)^2 ) / (2*V) - sum( (Y - thetarep)^2 ) / (2*W) )
    )
  }
}

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)
    thetarep = rep(theta,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 - thetarep)^2 ) / (2*W)
    )
  }
}

varcompquick = function(V, W, mu, theta, Y) {
    return( exp( logvarcompdens(V,W,mu,theta,Y) ) )
}

# SOME PRIOR DISTRIBUTION CONSTANTS
a1 = b1 = a3 = 0
a2 = 0.5
b2 = 1
b3 = 1e+12

# A TEST COMPUTATION
theta = c(3,2,1)
Y = c(3,4,5,6,7,8)
print( varcompdens(1,2,1,theta,Y) )
print( logvarcompdens(1,2,1,theta,Y) )
print( varcompquick(1,2,1,theta,Y) )

# THE FAMOUS "DYESTUFF DATA", FIRST PRESENTED IN DAVIES (1967):

# K=6, J=5
Ydye = 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)

thetatst = rep(1500,6)
mutst = 1550
Vtst=Wtst=60^2
print( varcompdens(Vtst,Wtst,mutst,thetatst,Ydye) )
print( varcompquick(Vtst,Wtst,mutst,thetatst,Ydye) )