2013年11月13日水曜日

Gnuplotによる区分関数の最小二乗法

不等号で定義域を分けて異なる関数を組み合わせた作った関数を区分関数（Piecewise Function）という。 例えば絶対値を用いた関数$f(x)=|x|$も$x=0$を境に$x>=0$における$f(x)=x$と、$x<0$における$f(x)=-x$の組み合わせと考えられる。 誤差を含むデータ点から区分関数のパラメータを求める最適化はGnuplotのfitコマンドで意外と簡単にできる。 面白いのは区分関数の区切り位置（2つの関数の境目）自体もパラメータとすることでその位置も求められること。

V字関数のフィッティング

データの作成が容易なので、次の簡単な区分関数のフィッティング方法を紹介する。 正確な名前は分からないがここではV字関数と仮名を付けておく。 この関数どういう形かというと、$a<0,\ b>0$のとき$x$軸上の$x=c$で$y=0$となり、折り返すまさにV字の関数となる。

\[ f(x)= \left\{ \begin{array}{ll} a(x-c) & (x < c)\\ b(x-c) & (x \geq c)\\ \end{array} \right. \]

上記関数に対して誤差を含むデータ点が与えられたとき、 Gnuplotのfitコマンドでパラメータ$a,b,c$を同時に非線形最適化して求めてみる。 使用したGnuplotのコマンドとデータは以下のとおり。

Gnuplotコマンド

reset

#関数の定義
f(x,a,b,c)= (x<c) ? a*(x-c) : b*(x-c)
 
#初期値の設定
a=0
b=0
c=0

#関数のフィッティング
fit f(x,a,b,c) 'sample.txt' via a,b,c
 
#表示設定
set xr[-10:10]
set samples 100
set terminal wxt

#関数とデータの描画
plot f(x,a,b,c) 
replot 'sample.txt' lc 3

データ "sample.txt"

-4 11
-3 7.5
-2 6
-1 3
-0 2
1 0
2 5
3 12
4 16
5 20
6 27

このデータは、$a=-2, b=5, c=1$の時の関数の値に対してわずかな誤差を人手により与えている。 最適化の初期値は意図的に正解値から遠ざけて$a=0, b=0, c=0$としたが、おおよそ上の正解値に近づいた結果$a=-2.0544, b=5.21874, c=0.936731$が得られている。

データとフィッティングされた関数

フィッティング結果

gnuplot> load 'fitting_sample.txt'
 
 Iteration 0
 WSSR        : 1780.25           delta(WSSR)/WSSR   : 0
 delta(WSSR) : 0                 limit for stopping : 1e-005
 lambda   : 1.91485
 
initial set of free parameter values
 
a               = 1e-030
b               = 1e-030
c               = 1e-030
/
　　　：
　　（省略）
　　　：
 
 Iteration 6
 WSSR        : 5.75979           delta(WSSR)/WSSR   : -9.82704e-009
 delta(WSSR) : -5.66017e-008     limit for stopping : 1e-005
 lambda   : 1.91485e-005
 
resultant parameter values
 
a               = -2.0544
b               = 5.21874
c               = 0.936731
 
After 6 iterations the fit converged.
final sum of squares of residuals : 5.75979
rel. change during last iteration : -9.82704e-009
 
degrees of freedom    (FIT_NDF)                        : 8
rms of residuals      (FIT_STDFIT) = sqrt(WSSR/ndf)    : 0.848513
variance of residuals (reduced chisquare) = WSSR/ndf   : 0.719974
 
Final set of parameters            Asymptotic Standard Error
=======================            ==========================
 
a               = -2.0544          +/- 0.1336       (6.501%)
b               = 5.21874          +/- 0.1976       (3.787%)
c               = 0.936731         +/- 0.1152       (12.3%)
 
correlation matrix of the fit parameters:
 
               a      b      c      
a               1.000 
b               0.403  1.000 
c               0.490  0.822  1.000・

参考

投稿者 時刻: