［R］周期変位を含むデータの回帰モデルを作成する

線形回帰であれば関数lmを使用すれば簡単にできる。以下の式で表される周期変動を含む回帰モデルを作成してみた例を示す。

y = a + b×x + c sin(2πx) + d cos(2πx)

a,b,c,d：定数

この例では周期変動の周期は1としている。以下に、回帰モデルの作成例を示す。

> # 独立変数(横軸方向)のデータを作成(特に意味はないが不等間隔)
> x <- c(0.0, 0.1, 0.5, 0.7, 1.0, 1.2, 1.4, 1.9, 2.0)
> # 従属変数(縦軸方向)のデータを適当に作成
> y <- c(0.1, 1.2, -0.4, -0.5, 0.0, 1.4, 0.4, 0.2, 1.5)
> # 関数lmを使用して回帰モデルを作成する
> r <- lm(y ~ x + sin(2 * pi * x) + cos(2 * pi * x))
> # 作成した回帰モデルを描画するために、描画のための横軸方向のデータを作成
> xest <- seq(0, 2, 0.05)
> # 描画のために作成した独立変数に基づく従属変数の値を作成
> yest <- predict(r, data.frame(x = xest))
> # 独立変数と従属変数のプロット
> plot(x,y)
> # 回帰モデルの描画
> lines(xest, yest)
> # 最後に作成した回帰モデルのモデルパラメーターを描画
> coef(r)
    (Intercept)               x sin(2 * pi * x) cos(2 * pi * x)
     -0.1074583       0.3660297       0.8998182       0.3977064
> # 切片
> coef(r)[1]
(Intercept)
-0.1074583
> # 三角関数の合成により、求めた式の振幅を計算する
> sqrt(coef(r)[3] ^ 2 + coef(r)[4] ^ 2)
sin(2 * pi * x)
      0.9837902

関数coefにより表示された値が、左から順にそれぞれa、b、c、dの値である。

sinとcosにそれぞれ「2 * pi * x」としたところがポイント。これは関数sinも関数cosも引数はラジアン単位で与えなければならないことから、xの値に2πをかけることで周期を1とした値に変換している。