## Some Box-Cox demonstrations.

## Return the profile likelihood for the Box-Cox coefficient over the 
## specified grid.
boxcox_profile = function(Y, X, minL=-2, maxL=2, grid=0.1)
{
  n = length(Y)
  qq = qr(X)
  Q = qr.Q(qq)
  V = NULL
  C = sum(log(Y))

  for (L in seq(-2.05, 2.05, 0.1))
  {
    Y1 = (Y^L - 1)/L
    Yhat = Q %*% (t(Q) %*% Y1)
    s2 = sum((Y1-Yhat)^2)/length(Y)
    v = -n*log(s2)/2 + (L-1)*C
    V = rbind(V, c(L, v))
  }
  return(V)
}

## Return the optimal value of the Box-Cox coefficient over the specified
## grid.
boxcox = function(Y, X, minL=-2, maxL=2, grid=0.1)
{
  V = boxcox_profile(Y, X, minL, maxL, grid)
  ii = which.max(V[,2])
  return(V[ii,1])
}


##
## Show the profile likelihood for lambda.
##

## The true value of lambda.
lambda = 0.5

## The sample size.
n = 100

## The number of regressors.
p = 10

## The covariate data.
X = array(runif((p+1)*n), c(n,p+1))
X[,1] = 1

## The regression coefficients.
beta = runif(p+1)
beta[1] = 1

## Generate a plot shoing 10 replicated draws of the Box-Cox profile 
## likelihood.
G = array(0, c(42,10))
for (k in 1:10)
{
  ## The transformed response (make sure all values are positive).
  Y_lambda = X %*% beta
  E = rnorm(n)
  while (TRUE)
  {
    ii = (Y_lambda + E <= 0)
    if (sum(ii) == 0) { break }
    E[ii] = rnorm(sum(ii)) 
  }
  Y_lambda = Y_lambda + E

  ## The observed response.
  Y = exp(log(lambda*Y_lambda + 1) / lambda)

  ## The Box/Cox estimate of lambda.
  V = boxcox_profile(Y, X)
  F = V[,1]
  G[,k] = V[,2]
}

pdf('box_cox_1.pdf')
plot(F, G[,1], t='l', xlab='Lambda', ylab='Profile likelihood', 
          ylim=c(min(G), max(G))) 
for (k in 2:10) { lines(F, G[,k]) }
dev.off()

##
## Show a histogram of estimated Box-Cox transformation parameters.
## 

## Generate a plot shoing 10 replicated draws of the Box-Cox profile 
## likelihood.
G = array(0, 1000)
for (k in 1:1000)
{
  ## The transformed response (make sure all values are positive).
  Y_lambda = X %*% beta
  E = rnorm(n)
  while (TRUE)
  {
    ii = (Y_lambda + E <= 0)
    if (sum(ii) == 0) { break }
    E[ii] = rnorm(sum(ii)) 
  }
  Y_lambda = Y_lambda + E

  ## The observed response.
  Y = exp(log(lambda*Y_lambda + 1) / lambda)

  ## The Box/Cox estimate of lambda.
  G[k] = boxcox(Y, X)
}

pdf('box_cox_2.pdf')
hist(G, main='Estimated Box-Cox parameter (true value is 1/2)', xlab='Lambda')
dev.off()