この記事は
Juliaでカーブフィッティングの三角関数をフィッティングする内容をCurveFitパッケージで実施しました.
実施環境
- macOS 10.14.6
- julia 1.4.1
データを作る
フィッティングされる方の曲線を作ります.今回は三角関数に誤差項が乗っているデータを扱います.
using Plots
using CurveFit
X = 0.01:0.01:1.0
Y = @. 0.5 + 0.4 * sin(2pi*X) + 0.05randn()
@.
は定義した数式をベクトルの全ての要素に適用するシンタックスです.
関数定義は非常に簡単にできます.また、2pi
のように数字を変数の前につけた場合は2*pi
と解釈されるので可読性が高くなります(さすがに2piX
はエラーになりました).randやsinなどの数学関数が組み込みなのが便利です.
このデータを見てみます.
gr()
scatter(X, Y)
サンプル
フィッティングに使用するサンプルを採取します.今回は等間隔に10点とります.
X0 = X[1:10:91]
Y0 = Y[1:10:91]
JuliaのインデックスはMatlabと同様に1から始まります.他の言語に慣れていると気持ち悪く感じるかもしれません.
フィッティング
CurveFitにはいくつかのフィッティングモデルがあるのですが,今回は多項式モデルのPolyを使用します.
fit1 = curve_fit(Poly, X0, Y0, 1)
Y1 = fit1.(X)
非常に簡単です.curve_fit
の第1引数にモデルのPoly
を指定し,第2, 3引数にフィッティングに使用するサンプルを指定します.第4引数は多項式の最大次数です.ここでは1次式でのフィッティングをしています.
第2行はフィッティングしたモデルを使って,近似線の値を計算しています.それではプロットしてみましょう.
plot(X, Y, seriestype = :scatter, label="data")
plot!(X, Y1, label="Poly deg1")
2行目のplot
に!
がついているのは,1行目のプロットに重ねて描画するためです(引数を変更する関数呼び出しを破壊的メソッドと呼ぶらしい).もし!
をつけないで2行目をplot
と書いた場合,1行目のプロット内容はリフレッシュされてしまい,2行目のプロット内容しか描画されません.
1次近似では無理がありますが,頑張ってますね.
3次近似、6次近似、9次近似
次数を増やしてみましょう.
fit3 = curve_fit(Poly, X0, Y0, 3)
Y3 = fit3.(X)
plot(X, Y, seriestype = :scatter, label="data")
plot!(X, Y3, label="Poly deg3")
3次近似だとこんな感じ.いい感じですね.
次は6次近似を試してみます.
fit6 = curve_fit(Poly, X0, Y0, 6)
Y6 = fit6.(X)
plot(X, Y, seriestype = :scatter, label="data")
plot!(X, Y6, label="Poly deg6")
ややオーバーフィッティングし始めましたね.
さらに次数を増やし9次にしてみます.
fit9 = curve_fit(Poly, X0, Y0, 9)
Y9 = fit9.(X)
plot(X, Y, seriestype = :scatter, label="data")
plot!(X, Y9, label="Poly deg9")
9次近似では見事なまでにオーバーフィッティングしてあらぬ方向に行っています.
以上,JuliaのCurveFit
を使ったカーブフィッティングを紹介しました.CurveFit
を使用すると,非常に簡単にフィッティングができて便利ですね.
参考にしたJuliaでカーブフィッティングでは10次近似をしているですが,curve_fit(Poly, X0, Y0, 10)
はエラーとなってしまい,フィッティングできませんでした.CurveFit
ではモデル式の変数がサンプル数より多い場合にはエラーを返す実装になっているようです.