# To get basically the same data set as described on pages 10-13 of textbook:

library(MASS)

set.seed(1203)

mymat1 <- mvrnorm(200, mu=c(0,0), Sigma=matrix(c(1,0.5,0.5,1),ncol=2))

mymat2 <- log(abs(mymat1)) + 12

# mymat2 can be thought of as the observed data matrix which is not multivariate normal.

chisplot(mymat2)

boxcox(lm(mymat2[,1]~1),lambda=seq(-10,10,1/5))
boxcox(lm(mymat2[,2]~1),lambda=seq(-10,10,1/5))


# To see the actual lambda values:
box.x <- boxcox(lm(mymat2[,1]~1),lambda=seq(-10,10,1/5),plotit=F)$x
box.y <- boxcox(lm(mymat2[,1]~1),lambda=seq(-10,10,1/5),plotit=F)$y
cbind(box.x, box.y)

box.x <- boxcox(lm(mymat2[,2]~1),lambda=seq(-10,10,1/5),plotit=F)$x
box.y <- boxcox(lm(mymat2[,2]~1),lambda=seq(-10,10,1/5),plotit=F)$y
cbind(box.x, box.y)


# Good choices for the lambdas:

lambda=c(5.5, 6)

mymat3 <- cbind( ((mymat2[,1])^5.5 - 1)/5.5,  ((mymat2[,2])^6 - 1)/6 )

# mymat3 is the transformed data matrix

#### A more efficient way to do the transformations of each column of mymat2.
#### This assumes all the lambda_j's are in a vector called "lambda", as above.

newmat <- NULL
for (j in 1:length(lambda)){
if (lambda[j]!=0){
newcol<- (airpol.data.sub[,j]^lambda[j] - 1)/lambda[j]
} else {
newcol<-log(airpol.data.sub[,j])
}
newmat <- cbind(newmat,newcol)
}
mymat3 <- newmat
colnames(mymat3) <-colnames(mymat2)

# mymat3 is the transformed data matrix
####
####

chisplot(mymat3)