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

# 데이터 생성: x ~ Uniform[0,1]에서 샘플, y = sin(2*pi*x) + 노이즈
set.seed(123)
n <- 100
x <- runif(n, 0, 1)
y <- sin(2 * pi * x) + 0.5 * rnorm(n)

plot(x,y,pch=16,col='blue')
f <- function(x) {sin(2 * pi * x)} 
curve(f,add=T,col='green',lwd=3)

y=sin(2pix)와 그의 산점도

K=4일 때를 보겠다. 

K가 뭐든간에 cos과 sin의 특징(값)을 추출해내야 한다.

이를 어떠한 행렬형태로 나타낼 것이다.

Cos(2pi*1*x)	Cos(2pi*2*x)	Cos(2pi*3*x)	Cos(2pi*4*x)	Cos(2pi*5*x)	Cos(2pi*6*x)
data11	data21	..	..	..	..
data12	data22	..	..	..	..

 이런 형태는 R에서 sapply를 이용해서 만들면 좋다. k=1,2,3,4... 일때 만들기에 특화되어있기 때문이다.

#K=4
K <- 4
cos_feat <- sapply(1:K,FUN =  function(k) cos(2*pi*k*x))
colnames(cos_feat) <- paste0('cos',1:K)
head(cos_feat)
> head(cos_feat)
           cos1       cos2       cos3        cos4
[1,] -0.2339190 -0.8905638  0.6505585  0.58620794
[2,]  0.2383614 -0.8863677 -0.6609131  0.57129536
[3,] -0.8408661  0.4141116  0.1444412 -0.65702312
[4,]  0.7418151  0.1005793 -0.5925926 -0.97976760
[5,]  0.9308533  0.7329757  0.4337325  0.07450688
[6,]  0.9593123  0.8405602  0.6534071  0.41308283

sin도 마저 작업하고 열을 기준으로 합치자.

sin_feat <- sapply(1:K, FUN =  function(k) sin(2*pi*k*x))
colnames(sin_feat) <- paste0('sin',1:K)

#col로 합치기.
fourier_feat <- cbind(cos_feat,sin_feat)
head(fourier_feat)
> head(fourier_feat)
           cos1       cos2       cos3        cos4       sin1       sin2       sin3       sin4
[1,] -0.2339190 -0.8905638  0.6505585  0.58620794  0.9722561 -0.4548583 -0.7594562  0.8101606
[2,]  0.2383614 -0.8863677 -0.6609131  0.57129536 -0.9711765 -0.4629820  0.7504625  0.8207445
[3,] -0.8408661  0.4141116  0.1444412 -0.65702312  0.5412432 -0.9102261  0.9895134 -0.7538704
[4,]  0.7418151  0.1005793 -0.5925926 -0.97976760 -0.6706045 -0.9949290 -0.8055023 -0.2001386
[5,]  0.9308533  0.7329757  0.4337325  0.07450688 -0.3653931 -0.6802548 -0.9010417 -0.9972205
[6,]  0.9593123  0.8405602  0.6534071  0.41308283  0.2823471  0.5417182  0.7570067  0.9106935

 

이제 회귀모형을 작성하자. 앞서 배운 행렬 A가 열공간에 속할 때 열벡터가 행렬 t(A)의 영공간에 속한다는 개념을 쓰지말고 이번엔 그냥 R에서 제공하는 패키지 lm()을 바로 이용해서 코드를 간략화 해보자.

> data <- data.frame(y,fourier_feat)
> head(data)
            y       cos1       cos2       cos3        cos4       sin1       sin2       sin3       sin4
1  1.09891536 -0.2339190 -0.8905638  0.6505585  0.58620794  0.9722561 -0.4548583 -0.7594562  0.8101606
2 -0.98544990  0.2383614 -0.8863677 -0.6609131  0.57129536 -0.9711765 -0.4629820  0.7504625  0.8207445
3  0.51980796 -0.8408661  0.4141116  0.1444412 -0.65702312  0.5412432 -0.9102261  0.9895134 -0.7538704
4  0.01369668  0.7418151  0.1005793 -0.5925926 -0.97976760 -0.6706045 -0.9949290 -0.8055023 -0.2001386
5 -0.47827861  0.9308533  0.7329757  0.4337325  0.07450688 -0.3653931 -0.6802548 -0.9010417 -0.9972205
6  1.04058245  0.9593123  0.8405602  0.6534071  0.41308283  0.2823471  0.5417182  0.7570067  0.9106935
> model <- lm(y~., data=data)
> summary(model)

Call:
lm(formula = y ~ ., data = data)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.98382 -0.30025 -0.02694  0.29285  0.95893 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -0.03987    0.04783  -0.834   0.4067    
cos1         0.07344    0.06922   1.061   0.2915    
cos2         0.01882    0.06967   0.270   0.7877    
cos3        -0.10207    0.06828  -1.495   0.1384    
cos4         0.07067    0.06720   1.052   0.2958    
sin1         1.07358    0.06559  16.368   <2e-16 ***
sin2        -0.10595    0.06533  -1.622   0.1083    
sin3        -0.04625    0.06676  -0.693   0.4902    
sin4        -0.15149    0.06821  -2.221   0.0288 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.469 on 91 degrees of freedom
Multiple R-squared:  0.7572,	Adjusted R-squared:  0.7358 
F-statistic: 35.47 on 8 and 91 DF,  p-value: < 2.2e-16

sin1일 때, 즉, sin(2pi*x) 일때 매~~~~우 유효하다고 나와있다. 

이제 모델지도도 끝났으므로 새 데이터로 예측을 잘하는지 확인해보자.

#Predict new data
x_new <- seq(0,1,length.out=100)

cos_feat_new <- sapply(1:K,FUN =  function(k) cos(2*pi*k*x_new))
colnames(cos_feat_new) <- paste0('cos',1:K)


sin_feat_new <- sapply(1:K, FUN =  function(k) sin(2*pi*k*x_new))
colnames(sin_feat_new) <- paste0('sin',1:K)

fourier_feat_new <- cbind(cos_feat_new,sin_feat_new)

#View(fourier_feat_new)

pred <- predict(model, newdata = data.frame(fourier_feat_new))
plot(x,y,pch=16,col='blue')
lines(x_new,pred,col='red',lwd=3)

f <- function(x) {sin(2 * pi * x)} 
curve(f,add=T,col='green',lwd=3)

 

 

 

원 함수에 어떻게 근사되고 있긴하다. MSE를 비교해보자면

#생성한 예측값에 대한 MSE
> (sum(y-pred))^2
[1] 9.97877
#기존 원함수에 대한 MSE
> (sum(y-f(x)))^2
[1] 7.221603

아니 그러면 기존 함수가 더 좋지 않냐고 할 수 있다.

하지만, 우리가 원함수를 예측 할 수 있는가? 보통 힘든일이 아닐 것이다. 

 

다음은 연습문제를 하나 남겨두겠다. 이에대한 정답코드도 같이 첨부해두겠다.

 

푸리에 회귀.pdf

0.02MB