UPDATE: 2022-10-23 11:48:53

はじめに

ここでは、タイトル通りフーリエ級数による曲線あてはめをRで行う。ちょっと調べる必要があって、フーリエ級数とか細かい話忘れていたので、おさらい。数式は手書きです… LaTeXで書くのが面倒だったので。ごめんなさい。

フーリエ級数

フーリエ級数は下記の三角関数の式のこと。赤枠がフーリエ係数。フーリエ級数は、複雑な周期をもつ関数であっても、単純な周期関数の無限の和によって表現できるとかいうスグレモノ。画像にある通り区間[a,b]で、データが2m+1個であれば、画像の2個目の式で近似できる。画像はiとnを書き間違えている。

ということで入門はじめての時系列分析のp52のサンプルデータを使ってRで再現する。この書籍はEXCELで再現する方法はのっているが、Rはないので、Rでやっておく。

df <- data.frame(y = c(2,3,5,6,7,7,6),
                 x = c(1,2,3,4,5,6,7))

plot(df$x, df$y)

このサンプルデータの場合、サイズは7、つまり2*3+1となりm=3で近似できる。また、フーリエ係数は、下記の画像の通り、行列の形で表現すると、線形回帰と同様に、逆行列を用いて、フーリエ係数を計算できる。

このサンプルデータから計算すると、このような式が得られる。画像はiとnを書き間違えている。

Rでも推定しておく。同じ結果が得られる。

n <- 1/8
fit <- lm(y ~ 
            cos(2*pi*n*x) + sin(2*pi*n*x) + 
            cos(4*pi*n*x) + sin(4*pi*n*x) +  
            cos(6*pi*n*x) + sin(6*pi*n*x),
          data = df)

fit
## 
## Call:
## lm(formula = y ~ cos(2 * pi * n * x) + sin(2 * pi * n * x) + 
##     cos(4 * pi * n * x) + sin(4 * pi * n * x) + cos(6 * pi * 
##     n * x) + sin(6 * pi * n * x), data = df)
## 
## Coefficients:
##         (Intercept)  cos(2 * pi * n * x)  sin(2 * pi * n * x)  
##           5.000e+00           -1.207e+00           -2.061e+00  
## cos(4 * pi * n * x)  sin(4 * pi * n * x)  cos(6 * pi * n * x)  
##          -1.321e-15           -5.000e-01            2.071e-01  
## sin(6 * pi * n * x)  
##          -6.066e-02

この係数を使って、さきほどの散布図に曲線を当てはめる。

pred <- function(x){
  res <- fit$coefficients[[1]] + #a
    fit$coefficients[[2]] * cos(2*pi*n*x) +
    fit$coefficients[[3]] * sin(2*pi*n*x) + 
    fit$coefficients[[4]] * cos(4*pi*n*x) + 
    fit$coefficients[[5]] * sin(4*pi*n*x) + 
    fit$coefficients[[6]] * cos(6*pi*n*x) + 
    fit$coefficients[[7]] * sin(6*pi*n*x)

  return(res)  
}

xx <- seq(0,8,0.1)
plot(df$x, df$y)
lines(x = xx, y = pred(x = xx), col = "red")

以上でおしまい。