통계학

Box-Cox 변환

써머23 2025. 5. 19. 19:43
728x90

 

-----------------------------

 

 

# 타구속도 데이터를 정규화 하기 위한 최적의 람다값 및 그 변환된 값을 구해보자

-------------------------

# 패키지 로드
library(readxl)
library(MASS)
library(ggplot2)

# 1. 데이터 불러오기
data <- read_excel("d:/mydata/PCA용.xlsx", sheet = "Sheet1")
y <- data$타구속도

# 2. Shift (음수나 0 방지)
if (any(y <= 0)) {
  shift <- abs(min(y)) + 1e-6
  y <- y + shift
  cat("Shift 적용됨\n")
} else {
  shift <- 0
}

# 3. Box-Cox 변환
model <- lm(y ~ 1)
bc <- boxcox(model, lambda = seq(-2, 2, 0.01))
lambda_max <- bc$x[which.max(bc$y)]
cat("최적 λ:", lambda_max, "\n")

# 4. 실제 변환 적용
if (lambda_max == 0) {
  y_trans <- log(y)
} else {
  y_trans <- (y^lambda_max - 1) / lambda_max
}

# ------------------------------
# 📊 A. 히스토그램 + 정규곡선
# ------------------------------
# 정규분포 곡선을 위해 평균과 sd 계산
mu <- mean(y_trans)
sigma <- sd(y_trans)

# 데이터 프레임화
df_trans <- data.frame(y_trans = y_trans)

# 히스토그램 + 정규곡선 그래프
ggplot(df_trans, aes(x = y_trans)) +
  geom_histogram(aes(y = ..density..), bins = 15, fill = "skyblue", colour = "black") +
  stat_function(fun = dnorm, args = list(mean = mu, sd = sigma), 
                colour = "red", size = 1.2) +
  labs(title = "Box-Cox 변환된 데이터의 히스토그램과 정규곡선",
       x = "Box-Cox 변환 값", y = "밀도") +
  theme_minimal()

# ------------------------------
# 📈 B. Q-Q Plot (정규분포 비교)
# ------------------------------
qqnorm(y_trans, main = "Box-Cox 변환 데이터의 Q-Q Plot")
qqline(y_trans, col = "red", lwd = 2)

 

 

 

# 다변량 벡터의 경우

-----------------------------------------------------------------------------------------------------------------------------

# 필요한 패키지 로드
library(readxl)
library(MASS)
library(ggplot2)
library(gridExtra)  # 여러 그래프 배치
library(tidyr)
library(dplyr)
library(nortest) 

# 1. 데이터 불러오기
data <- read_excel("d:/mydata/PCA용.xlsx", sheet = "Sheet1")
vars <- c("컨택률", "스윙률", "타구속도")

# 결과 저장 리스트 초기화
results <- list()
plots <- list()

# 반복: 각 변수마다 처리
for (var in vars) {
  cat("====", var, "====\n")
  
  y <- data[[var]]
  
  # 음수/0 방지 shift
  if (any(y <= 0)) {
    shift <- abs(min(y)) + 1e-6
    y_shifted <- y + shift
  } else {
    shift <- 0
    y_shifted <- y
  }
  
  # 평균모형으로 Box-Cox 수행
  model <- lm(y_shifted ~ 1)
  bc <- boxcox(model, lambda = seq(-2, 2, 0.01))
  lambda_max <- bc$x[which.max(bc$y)]
  cat("최적 λ:", round(lambda_max, 4), "\n")
  
  # 변환 수행
  if (lambda_max == 0) {
    y_trans <- log(y_shifted)
  } else {
    y_trans <- (y_shifted^lambda_max - 1) / lambda_max
  }
  
  # 저장
  results[[var]] <- list(original = y, transformed = y_trans, lambda = lambda_max)
  
  # 히스토그램: 원본
  mu_orig <- mean(y)
  sd_orig <- sd(y)
  p1 <- ggplot(data.frame(y), aes(x = y)) +
    geom_histogram(aes(y = ..density..), bins = 15, fill = "gray", color = "black") +
    stat_function(fun = dnorm, args = list(mean = mu_orig, sd = sd_orig), color = "red", size = 1) +
    labs(title = paste0(var, " (원본)"), x = var, y = "밀도") +
    theme_minimal()
  
  # 히스토그램: 변환 후
  mu_tr <- mean(y_trans)
  sd_tr <- sd(y_trans)
  p2 <- ggplot(data.frame(y_trans), aes(x = y_trans)) +
    geom_histogram(aes(y = ..density..), bins = 15, fill = "skyblue", color = "black") +
    stat_function(fun = dnorm, args = list(mean = mu_tr, sd = sd_tr), color = "red", size = 1) +
    labs(title = paste0(var, " (Box-Cox 변환 후, λ=", round(lambda_max, 2), ")"), x = "변환값", y = "밀도") +
    theme_minimal()
  
  # Q-Q plot: 원본
  qq_original <- function() {
    qqnorm(y, main = paste0(var, " 원본 Q-Q Plot"))
    qqline(y, col = "red", lwd = 2)
  }
  
  # Q-Q plot: 변환 후
  qq_transformed <- function() {
    qqnorm(y_trans, main = paste0(var, " Box-Cox 변환 Q-Q Plot"))
    qqline(y_trans, col = "red", lwd = 2)
  }
  
  # 저장
  plots[[var]] <- list(hist_before = p1, hist_after = p2,
                       qq_before = qq_original, qq_after = qq_transformed)
}

# =========================
# 📊 히스토그램 결과 보기
# =========================
# 예시로 컨택률부터 출력
grid.arrange(plots[["컨택률"]]$hist_before,
             plots[["컨택률"]]$hist_after,
             ncol = 2)

# 예시로 스윙률 출력
grid.arrange(plots[["스윙률"]]$hist_before,
             plots[["스윙률"]]$hist_after,
             ncol = 2)

# 예시로 타구속도 출력
grid.arrange(plots[["타구속도"]]$hist_before,
             plots[["타구속도"]]$hist_after,
             ncol = 2)

# =========================
# 📈 Q-Q Plot 출력 (Base R)
# =========================
# 원하는 변수 Q-Q Plot 확인
par(mfrow = c(1, 2))  # 한 줄에 두 개씩 출력
plots[["컨택률"]]$qq_before()
plots[["컨택률"]]$qq_after()

par(mfrow = c(1, 2))
plots[["스윙률"]]$qq_before()
plots[["스윙률"]]$qq_after()

par(mfrow = c(1, 2))
plots[["타구속도"]]$qq_before()
plots[["타구속도"]]$qq_after()

# 5. Shapiro-Wilk Test
shapiro_result <- shapiro.test(y_trans)
cat("\n[Shapiro-Wilk Test]\n")
print(shapiro_result)

# 6. Anderson-Darling Test
ad_result <- ad.test(y_trans)
cat("\n[Anderson-Darling Test]\n")
print(ad_result)

 

728x90