曲線の近似, fit

References

概要

fit コマンドはデータ点列を説明 (近似) するような曲線を探索する.

f(x) = g(x, a, b, c)

fit f(x) <expression> via a, b, c

params	default	value	explanation
`f(x)`	--	(function)	近似に使う関数
`expression`	--	(string)	近似に使うデータソース

近似のモデルとする関数 \(f(x)\) を具体的に定義する. 上では関数 \(g(x,a,b,c)\) とした. ここで関数 \(g\) には自由変数 (ここでは a,b,c ) を含めることができる. fit コマンドは最も良くデータを説明する a,b,c を探索する.

内部的には Levenberg-Marquardt アルゴリズムを使って最小二乗法を実行している.

Examples

正弦波を三次関数で近似する

ノイジーな正弦曲線を描く Ruby スクリプトをデータの例に使う.

def noisy_sin(x)
  return Math.sin(x) + (rand - 0.5) * 0.2
end

100.times do |i|
  x = (i - 50) / 20.0
  puts "#{x} #{noisy_sin x}"
end

このデータを近似する関数を定義する. 例えば \(f(x) = ax^3 + bx + c\) という三次多項式で近似してみることにする. ここで \(a,b,c\) は自由変数.

自由変数は \(1.0\) の初期値が使われるが, 自分で初期化してもよい.

f(x) = a * x**3 + b * x + c

# 初期値の設定
a = 1
b = 2
c = 1

# fitting の実行
fit f(x) '<ruby -e "100.times do |i| x = (i - 50) / 20.0; puts \"#{x} #{Math.sin(x) + (rand - 0.5) * 0.2}\" end"' via a,b,c

# この時点で a,b,c には fitting した値が入っている
title = sprintf("Fitted Curve (a=%.2f, b=%.2f, c=%.2f)", a, b, c)

plot '<ruby -e "100.times do |i| x = (i - 50) / 20.0; puts \"#{x} #{Math.sin(x) + (rand - 0.5) * 0.2}\" end"' title 'Noisy Sine', \
     f(x) title title

NOTE 上の例では ruby を二度実行してるので fit に用いたデータと plot に用いたデータは異なる.

三角波を正弦波で近似する

参考
- 三角波のフーリエ級数展開

\(-\pi \leq x \leq \pi\) の範囲で

\[f(x) = \begin{cases} -\pi - x & (-\pi \leq x \lt -\pi/2) \\ x & (-\pi/2 \leq x \leq \pi/2) \\ \pi - x & (\pi/2 \lt x \leq \pi) \\ \end{cases}\]

という関数（を周期関数にしたもの）を三角波という. 参考ページに有るように

\[f(x) = \sum_{k=0}^\infty a_k \sin((2k+1)x)\]

という形式の正弦波の重ね合わせで表現できる. これを題材にやってみる.

set xrange [-pi:3.1416]
set yrange [-1.7:1.7]
set xtics ("-pi" -pi, "-pi/2" -pi/2, "0" 0, "pi/2" pi/2, "pi" pi)
set ytics ("-pi/2" -pi/2, "0" 0, "pi/2" pi/2)
set grid
set key right bottom

g(x) = a0 * sin(x) + a1 * sin(3*x) + a2 * sin(5*x)
fit g(x) '<ruby -e "(-3.1..3.1).step(0.1) { |x| f = (x <= -Math::PI/2) ? (-Math::PI - x) : (x <= Math::PI/2) ? x : (Math::PI - x); puts \"#{x} #{f}\" }"' via a0, a1, a2

title = sprintf("Fitted Curve (a0=%.2f, a1=%.2f, a2=%.2f)", a0, a1, a2)

plot '<ruby -e "(-3.1..3.1).step(0.1) { |x| f = (x <= -Math::PI/2) ? (-Math::PI - x) : (x <= Math::PI/2) ? x : (Math::PI - x); puts \"#{x} #{f}\" }"' title 'Triangle Wave', \
     g(x) title title