[R] AirPassengers 데이터를 이용한 회귀분석 비교(회귀분석,조화회귀)

위 데이터를 다루는데 있어 어떠한 모형인지부터 알고 시작하자.

data('AirPassengers')
AP <- AirPassengers
> AP
     Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
1949 112 118 132 129 121 135 148 148 136 119 104 118
1950 115 126 141 135 125 149 170 170 158 133 114 140
1951 145 150 178 163 172 178 199 199 184 162 146 166
1952 171 180 193 181 183 218 230 242 209 191 172 194
1953 196 196 236 235 229 243 264 272 237 211 180 201
1954 204 188 235 227 234 264 302 293 259 229 203 229
1955 242 233 267 269 270 315 364 347 312 274 237 278
1956 284 277 317 313 318 374 413 405 355 306 271 306
1957 315 301 356 348 355 422 465 467 404 347 305 336
1958 340 318 362 348 363 435 491 505 404 359 310 337
1959 360 342 406 396 420 472 548 559 463 407 362 405
1960 417 391 419 461 472 535 622 606 508 461 390 432

각 데이터를 일렬로 늘어뜨리고 plot를 그려보자

df <- data.frame(t=t,y=as.numeric(AP))
head(df)
plot(df$t,df$y,type='l')

위로 상승하는 트렌드가 있고 어떠한 오르락내리락하는 주기성을 가진다.

우선 주기성을 판단했으므로 동적회귀(harmonic reg)을 사용하고 일반회귀랑 비교를 해보자.

주의할 점은 훈련,검증 데이터를 주기성을 가지기에 순서대로 짤라야 함을 유의하자. 

#훈련,검증 데이터 분리 7:3
set.seed(4121)

#시계열 데이터이므로 sample()과 같은 랜덤형 방식은 X 주기성을 꼭 고려해줘야 하므로 순서대로..
train_idx <- 1:round(0.7*N)
test_idx <- (round(0.7*N)+1) : N

train_df <- df[train_idx,]
test_df <- df[test_idx,]

https://pastryofjsmath.tistory.com/90

[R] 푸리에급수를 이용한 회귀분석(harmonic regression)

여기에도 있듯이 우리는 베타값을 추정하는데 있어 cos,sin의 값을 이용해야 한다. 즉, 얘네가 계수가 됨을 인지하자.

create_fourier_feat <- function(time,K,period){
  cos_terms <- sapply(1:K, FUN = function(k) cos(2 * pi * k * time / period))
  sin_terms <- sapply(1:K, FUN = function(k) sin(2 * pi * k * time / period))
  #열 병합
  fourier_feat <- cbind(cos_terms,sin_terms)
  colnames(fourier_feat) <- c(paste0('cos',1:K),paste0('sin',1:K))
  return(fourier_feat)
}
period <- 12
K <- 2

#훈련데이터에 적용
fourier_train <- create_fourier_feat(train_df$t, K = 2, period = 12)
fourier_train_df <- cbind(y=train_df$y, fourier_train)

#검증데이터에 적용
fourier_test <- create_fourier_feat(test_df$t, K = 2,period = 12)
fourier_test_df <- cbind(y=test_df$y, fourier_test)

 

계수행렬을 모두 만들었으므로 이제 푸리에 회귀모형을 작성하자.

#푸리에 모델생성
model_fourier <- lm(y ~ ., data = data.frame(fourier_train_df))

이후 검증 모델에서도 적용시키고 MSE를 구해보자.

pred_fourier <- predict(model_fourier, newdata = as.data.frame(fourier_test))
mse_fourier <- mean((test_df$y - pred_fourier)^2)
mse_fourier
> mse_fourier
[1] 44475.54

 

이제 기본선형회귀모형도 바로 작성하자. 앞에서의 train_df, test_df를 이용하면 된다.

#기본 선형회귀모델 생성
#일반 회귀분석: 단순히 시간 t만을 사용한 선형 회귀모델 학습
model_ordinary <- lm(y ~ t, data = train_df)
pred_ordinary <- predict(model_ordinary, newdata = test_df)
mse_ordinary <- mean((test_df$y - pred_ordinary)^2)

> # 결과
> cat("Fourier 회귀분석 MSE:", mse_fourier)
Fourier 회귀분석 MSE: 44475.54
> cat("일반 선형 회귀 MSE:", mse_ordinary)
일반 선형 회귀 MSE: 5783.176

t만을 독립변수로 한 회귀모형보다 성능이 매우 하찮다...

이를 plot으로 그려보자.

# 예측 결과 시각화
plot(df$t, df$y, type = "l", main = "AirPassengers 회귀(동적,일반선형) 예측 비교",
     xlab = "시간 (t)", ylab = "승객 수", col = "grey")
lines(test_df$t, pred_fourier, col = "blue", lwd = 2, lty = 2)
lines(test_df$t, pred_ordinary, col = "red", lwd = 2, lty = 2)
legend("topleft", legend = c("실제값", "Fourier 회귀 예측", "일반 선형회귀 예측"),
       col = c("grey", "blue", "red"), lty = c(1,2,2), lwd = 2)

구릴수밖에 없다. 우상향하는 트렌드를 반영하지 못했기 때문이다. 그렇다면 어떻게 하면 좋을까?

우상향하는 트렌드를 납작하게 짓눌러버리면 이쁜 모양이 나올 것 같다.

납작하게 눌러버리는 그러한 방식이 로그를 씌우면 해결된다. 그럼 실제값하고 차이가 난다고 할 수 있지만 나중에 지수를 씌워주면 문제가 없다.이를 트렌드와 주기성을 고려한 모델이라고 하자.

 

이는 y를 logy로 변환해서 위의 코드대로 작성하면 된다.

data("AirPassengers")
AP <- AirPassengers
N <- length(AP)
t <- 1:N
log_y <- log(as.numeric(AP))  # 로그 변환

# 데이터프레임
df <- data.frame(
  t = t,
  log_y = log_y
)

# 훈련/테스트 분할
train_idx <- 1:round(0.8 * N)
test_idx <- (round(0.8 * N) + 1):N
train_df <- df[train_idx, ]
test_df  <- df[test_idx, ]

# Fourier 기저 생성
period <- 12
K <- 2  # 필요에 따라 조정
create_fourier_feat <- function(time, K, period) {
  cos_terms <- sapply(1:K, function(k) cos(2 * pi * k * time / period))
  sin_terms <- sapply(1:K, function(k) sin(2 * pi * k * time / period))
  cbind(cos_terms, sin_terms)
}
fourier_train <- create_fourier_feat(train_df$t, K, period)
fourier_test  <- create_fourier_feat(test_df$t,  K, period)

# 데이터프레임으로 변환
fourier_train_df <- data.frame(
  log_y = train_df$log_y,
  t = train_df$t,
  fourier_train
)
fourier_test_df <- data.frame(
  log_y = test_df$log_y,
  t = test_df$t,
  fourier_test
)

# 회귀 모델 (로그 변환된 y, 추세 t, Fourier)
model_fourier_trend <- lm(log_y ~ t + ., data = fourier_train_df)

# 예측 (로그 스케일)
pred_log <- predict(model_fourier_trend, newdata = fourier_test_df)
# 지수화하여 원래 스케일로 복원
pred <- exp(pred_log)

# 실제 y값
actual_y <- exp(test_fourier_df$log_y)

# MSE 계산
mse_fourier_trend <- mean((actual_y - pred)^2)
cat("fourier + trend + log scale MSE:", mse_fourier_trend)

# 시각화
plot(df$t, exp(df$log_y), type = "l", col = "grey",
     main = "AirPassengers (log scale + trend + fourier)",
     xlab = "Time", ylab = "Passengers")
lines(test_df$t, pred, col = "blue", lwd = 2,lty=2)
legend("topleft", legend = c("실제값", "예측값"), 
       col = c("grey", "blue"), lty = c(1,2), lwd=2)

트렌드까지 고려한 푸리에 회귀모형