#title 변수 변환
[[TableOfContents]]

==== power 변환 ====
아래와 같은 지수 분포가 있다.
{{{
set.seed(1)
x <- rexp(1000,1)
hist(x)
}}}

아래와 같이 해서 변환한다. 뭔가 병신같지만 쓸만하다. 
{{{
pwr <- function(x){
    x <- sample(x, 500)
    pval <- 0
    for(p in seq(-5, 5, 0.01)){
        if(p != 0) {
            tmp <- shapiro.test(x^p)$p.value
        }
        
        if(pval <= tmp) {
            pval <- tmp
        }
    }
    return(pval)
}
hist(x^pwr(x))
}}}

==== 정규분포 확인 ====
shapiro.test로 확인 한다. p-value가 0.05이면 정규분포라고 본다.
{{{
data(trees)

shapiro.test(trees$Girth)
shapiro.test(trees$Volume)
}}}

결과
{{{
> shapiro.test(trees$Height)

	Shapiro-Wilk normality test

data:  trees$Height
W = 0.9655, p-value = 0.4034

> shapiro.test(trees$Volume)

	Shapiro-Wilk normality test

data:  trees$Volume
W = 0.8876, p-value = 0.003579
}}}

Heigth는 정규분포이고, Volume은 정규분포가 아니다. log변환을 해보자.

==== log 변환 ====
{{{
shapiro.test(log(trees$Volume))
hist(log(trees$Volume))
}}}

결과
{{{
> shapiro.test(log(trees$Volume))

	Shapiro-Wilk normality test

data:  log(trees$Volume)
W = 0.9643, p-value = 0.3766
}}}
log 변환을 했더니 정규분포스러워졌다. boxcox 변환을 해볼까?

==== boxcox 변환 ====
그냥 만들어 봤다. shapiro.test를 이용하므로 5000개 이하의 변수값만 된다.

공식 = http://onlinestatbook.com/2/transformations/box-cox.html
{{{
boxcox.f <- function(x)
{
  p <- c()
  l <- c()
  for(lambda in seq(-2,2,0.1))
  {
    if (lambda != 0)
    {
      p.tmp <- shapiro.test((x^lambda-1) / lambda)$p.value
      if (p.tmp >= 0.05) 
      {
        l <- c(l, lambda)
        p <- c(p, p.tmp)
      }
    }
  }
  d <- data.frame(p,l)
  #names(d)[names(d)=="p"] <- "pvalue"
  #names(d)[names(d)=="l"] <- "lambda"
  #d[d$p==min(d$p)|d$p==max(d$p), c("pvalue", "lambda")]
  max.lambda <- d[d$p==max(d$p), c("l")]
  rs <- (x^max.lambda-1) / max.lambda
  return (rs)
}

x <- boxcox.f(trees$Volume)
}}}

==== boxcox 변환:회귀분석 ====
회귀진단으로 잔차에 대한 분석을 하는데, 잔차는
 * 정규분포를 따라야하고,  --> shapiro.test(resid(model))
 * 분산이 일정하고, --> bptest(model)
 * 특별한 추세를 보이지 않아야 한다. --> dwtest(model)
즉, 잔차는 정규성, 등분산성, 독립성을 만족해야 한다. 참고로 bptest, dwtest는 lmtest library를 사용한다.

먼저 변수 변환을 하지 않은 상태와 로그변환을 한 상태에서 회귀분석을 해보자. 
 * model <- lm(Volume~Height+Girth,data=trees)
 * model <- lm(I(log(Volume))~Height+Girth,data=trees)
둘다 독립성을 만족시키지 못한다. 그래서 box-cox변환을 실시해보자. 

{{{
library("MASS") 
boxcox(lm(Volume~Height+Girth,data=trees),lambda=seq(-1,1,by=.1))
}}}
attachment:변수변환/boxcox.png

그림을 보면 알겠지만 대략 0.3에서 lambda가 결정된다는 것을 볼 수 있다. 최적의 lambda 값을 찾아보자.
{{{
bc <- boxcox(lm(Volume~Height+Girth,data=trees),lambda=seq(-1,1,by=.1))
lambda <- bc$x[which.max(bc$y)]
}}}

결과
{{{
> bc <- boxcox(lm(Volume~Height+Girth,data=trees),lambda=seq(-1,1,by=.1))
> lambda <- bc$x[which.max(bc$y)]
> lambda
[1] 0.3030303
}}}

{{{
lambda <- 0.3
model <- lm(I((Volume^lambda - 1)/lambda)~Height+Girth,data=trees)
shapiro.test(resid(model))

#install.packages("lmtest")
#library("lmtest")
bptest(model)
dwtest(model)
}}}

결과
{{{
> shapiro.test(resid(model))

	Shapiro-Wilk normality test

data:  resid(model)
W = 0.96822, p-value = 0.4714

> library("lmtest")
필요한 패키지를 로딩중입니다: zoo

다음의 패키지를 부착합니다: ‘zoo’

The following objects are masked from ‘package:base’:

    as.Date, as.Date.numeric

> dwtest(model)

	Durbin-Watson test

data:  model
DW = 2.069, p-value = 0.4811
alternative hypothesis: true autocorrelation is greater than 0

> dwtest(model)

	Durbin-Watson test

data:  model
DW = 2.069, p-value = 0.4811
alternative hypothesis: true autocorrelation is greater than 0

> 
}}}
잔차의 정규성, 등분산성, 독립성을 모두 만족한다.

참고: 변수변환한 것으로 다시 원래대로..
{{{
b <- 100
lambda <- 0.2
y <- (b^lambda - 1)/lambda
((y*lambda) + 1)^(1/lambda)
}}}

==== 푸리에 변환 ====
 * [푸리에 변환]

==== 웨이블릿 변환 ====
 * attachment:변수변환/wavelet_transform_in_r.pdf
 * https://atmos.washington.edu/~breth/classes/AS552/matlab/lect/html/wavelet_turbulence.html

==== 공간 변환 ====
{{{
x <- rnorm(100, 80, 10)
y <- rnorm(100, 80, 10)

mydata <- data.frame(x, y, grp="a")
mydata <- rbind(mydata, data.frame(x = rnorm(10, 30, 5), y = rnorm(10, 30, 5), grp = "b"))
head(mydata)
plot(y~x, data=mydata, col=grp)

mydata$x <- scale(mydata$x)
mydata$y <- scale(mydata$y)

mydata$x1 <- mydata$x / sqrt(mydata$x^2 + mydata$y^2)
mydata$y1 <- mydata$y / sqrt(mydata$x^2 + mydata$y^2)
plot(mydata$x1)
plot(mydata$y1)

par(mfrow=c(1,2))
plot(mydata$x, mydata$y, main="이상치가 있는원래 데이터", col=mydata$grp)
plot(mydata$x1, mydata$y1, main="공간 변형", col=mydata$grp)
par(mfrow=c(1,1))
}}}
attachment:변수변환/t1_2.png