Hello all, Given a set of observations, I would like to find lambda for a boxcox transformation so to get a symmetric/normal result as possible.
I try to use MASS::boxcox, but get different results than when using the formula from the original box-cox paper (link <http://www.jstor.org/stable/2984418?seq=1#page_scan_tab_contents>). I probably have made an error somewhere, but I can't figure out where. Here is an example in which the lambda by MASS::boxcox is 0.42424, while by the formula from the paper I get 0.40782. # Toy data ################ set.seed(13241089) x <- rnorm(1000, 10) x2 <- x**2 # we want to transform x2 to something more normal plot(density(x2)) # using MASS::boxcox ################ zpoints <- function(y) { n <- length(y) qnorm(ppoints(n))[order(order(y))] } mle <- function(BC) { with(BC, x[which.max(y)]) } a <- MASS::boxcox(x2 ~ zpoints(x2)) mle(a) # lambda: # 0.42424 # using formula from the paper ################ loglik_lambda <- function(l, y) { GM <- exp(mean(log(y))) if(l==0) x <- log(y)*GM else x <- (y^l-1)/ (l * GM^(l-1) ) # if(l==0) x <- log(y) else x <- (y^l-1)/ (l ) sd(x) } fo <- function(l) loglik_lambda(l, y = x2) V_fo <- Vectorize(fo) V_fo(2) curve(V_fo, -.5,1.5) optimize(V_fo, c(-3,3)) # lambda: # 0.40782 [[alternative HTML version deleted]] ______________________________________________ R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.