BoxCox <- function(y, lambda)
#
#  Carry out Box-Cox transformation of variable y.
#
	if(lambda == 0) log(y) else (y^lambda - 1)/lambda

BoxCoxLike <- function(Lambda=0, GLMfit, OrigY) {
#
#   Compute the loglikelihood of data when the Box-Cox 
#   transformed data have a normal(mu, sigma^2) distribution.
#
  RSS <- deviance(GLMfit)
  n <- length(GLMfit$y)
  sigma.hat <- deviance(GLMfit)/GLMfit$df.residual
  Logli <- (Lambda-1)*sum(log(OrigY)) - n*log(sqrt(sigma.hat)) 
  -2*Logli
}
BoxCoxPL <- function(Lambda=0, GLMfit, OrigY) {
#
#   Compute the profile likelihood of data when the Box-Cox 
#   transformed data have a normal(mu, sigma^2) distribution.
#
  RSS <- deviance(GLMfit)
  n <- length(GLMfit$y)
  pl <- -((n/2)*log(RSS)) + ((Lambda-1)*sum(log(OrigY)))
  ProfLike <- pl
  -2*ProfLike
}

#
#   Example computing the loglikelihood, etc. with 
#   Box-Cox transformations carried out on a grid.
#
lambda <- seq(from=-2, to=2, by=0.1)
pl.vec <- rep(NA, length(lambda))
logli.vec <- pl.vec
likeli.vec <- pl.vec
y <- energy
for (i in 1:length(lambda)) {
  l <- lambda[i]
  y.l <- BoxCox(y,l)
  fit <- glm(y.l ~ weight, family=gaussian, link=identity)
  pl.vec[i] <- BoxCoxPL(Lambda=l, GLMfit=fit, OrigY=y)
  logli.vec[i] <- BoxCoxLike(Lambda=l, GLMfit=fit, OrigY=y)
  likeli.vec[i] <- exp(-logli.vec[i]/2)
}
#
#   Plot -2*loglikelihood as a function of the Box-Cox 
#   transformation parameter "lambda."
#
plot(lambda, logli.vec, type="b", xaxt="n",
     xlab="Lambda", ylab="-2*loglikelihood")
axis(1, at=(-20:20)/10, labels=F, tck=-0.02)
axis(1, at=(-4:4)/2, tck=-0.04)

plot(lambda, likeli.vec, type="b", xaxt="n",
     xlab="Lambda", ylab="likelihood")
axis(1, at=(-20:20)/10, labels=F, tck=-0.02)
axis(1, at=(-4:4)/2, tck=-0.04)