# copy this text and paste it into an R session,
# or just type
# source("http://www.stat.ncsu.edu/people/bloomfield/courses/st731/boxcox.r")
# at the R prompt.
boxcoxplot = function(x, lambdalo = -2, lambdahi = 2, lambdastep = 0.1,
                      n = (lambdahi - lambdalo) / lambdastep + 1,
                      useMASS = TRUE, ...)
{
# graph the Box-Cox likelihood for a single variable,
# as in Figure 4.12 of Johnson and Wichern
    lambda = seq(from = lambdalo, to = lambdahi, length = n);

    if (useMASS && require(MASS))
        return(boxcox(x ~ 1, lambda, ...));

    lik = rep(NA, n);
    for (i in 1:n)
        lik[i] = (-1/2) * boxcoxlikelihood(x, lambda[i]);
    plot(lambda, lik, type = "l");

    lmax = max(lik);
    lmin = min(lik);
    li = min((1:n)[lik == lmax]);
    lines(rep(lambda[li], 2), c(lmin, lmax), lty = 2);
    lt = format(lambda[li]);
    text(lambda[li], lmin, lt, srt = 90, c(-0.1, -0.1));

    lc = lmax + (-1/2) * qchisq(.95, 1);
    li = (1:n)[lik > lc];
    abline(h = lc, lty = 2);
    lines(rep(lambda[min(li)], 2), c(lmin, lc), lty = 2);
    lt = format(lambda[min(li)]);
    text(lambda[min(li)], lmin, lt, srt = 90, c(-0.1, -0.1));
    lines(rep(lambda[max(li)], 2), c(lmin, lc), lty = 2);
    lt = format(lambda[max(li)]);
    text(lambda[max(li)], lmin, lt, srt = 90, c(-0.1, -0.1));
}

boxcoxlikelihood = function(x, lambda, tol = 1e-6)
{
# calculate -2 x log Likelihood (equation 4-40 in Johnson and Wichern)
# x:         n x p data matrix
# lambda:    p-vector of exponents, or a single common exponent
    y = as.matrix(x);
    if (length(lambda) == 1)
        lambda = rep(lambda, ncol(y));
    for (j in 1:ncol(y))
        y[,j] = if (abs(lambda[j]) > tol)
            (y[,j]^lambda[j] - 1)/lambda[j]
            else log(y[,j]);
    ybar = apply(y, 2, mean);
    z = y - rep(1, nrow(y)) %o% ybar;
    n = nrow(z);
    sz = t(z) %*% z / n;
    n * log(det(sz)) - 2 * sum((lambda - 1) * t(log(x)));
}

boxcoxcontour = function(x, lambda1lo = -1, lambda1hi = 1,
     lambda2lo = lambda1lo, lambda2hi = lambda1hi,
     n1 = 30, n2 = n1)
{
# x            n x 2 data matrix
# lambda1lo    low end of grid for lambda1
# ...
# n1, n2       length of lambda1 and lambda2 grids
    lambda1 = seq(from = lambda1lo, to = lambda1hi, length = n1)
    lambda2 = seq(from = lambda2lo, to = lambda2hi, length = n2)
    lik = matrix(NA, n1, n2)
    for (i in 1:n1)
        for (j in 1:n2)
            lik[i, j] = boxcoxlikelihood(x, c(lambda1[i], lambda2[j]));
    lik = lik - min(lik);
    contour(lambda1, lambda2, lik);
    contour(lambda1, lambda2, lik, levels = qchisq(.95, 2),
            add = T, lty = 2);
    for (i in 1:n1)
        for (j in 1:n2)
            if (lik[i, j] <= 0) {
                points(lambda1[i], lambda2[j]);
                text(lambda1[i], lambda2[j],
                     paste("(",
                           format(round(lambda1[i], 3)), ",",
                           format(round(lambda2[j], 3)), ")"),
                     adj = 1.1);
            }
    abline(h = 1, v = 1);
}

boxcoxestimate = function(x, lambda = rep(1, ncol(x)), tol = 1e-6, ...)
{
# estimate Box-Cox transformation exponents
# x:         n x p data matrix
# lambda:    p-vector of starting values, or a single starting value
#            for a common exponent
    optim(lambda, function(l) boxcoxlikelihood(x, l, tol = tol), ...)$par;
}