綾小路龍之介の素人思考

[数値計算] 最小二乗法フィルタ

あるデータファイルに記載されたxとyの関係(実験の測定値)が次の関数でを満たしていると仮定するんだな。y=a*x+b。これから最小二乗法でaとbを定める事を考えるんだな。例えばデータファイルは「xの値、スペース、yの値、改行」のような書式となっているとするんだな。

とりあえずデータファイルを次のようにして生成したんだな。

C:\>perl -le "sub f(){$a=1+rand(0.1);$b=1+rand(0.1);return $a*$_+$b;} for(1..20){$x=$_+rand(0.1);$f=&f($x);print \"$x $f\";}">a.dat
C:\>

これにより生成したa.datの内容を表示すると下のようになるんだな。

C:\>perl -pe "" a.dat
1.01757202148438 2.11520385742188
2.05708923339844 3.05591735839844
3.03453979492188 4.1234130859375
4.03664245605469 5.36311340332031
5.02697448730469 6.1601806640625
6.01442565917969 7.09348449707031
7.04530944824219 8.60143432617187
8.09518127441406 9.57536010742187
9.06552429199219 10.8480926513672
10.0904388427734 11.7131164550781
11.0158996582031 12.6625122070312
12.0860198974609 14.0559478759766
13.0393676757813 14.6789276123047
14.0699066162109 16.3958953857422
15.0349884033203 16.2184478759766
16.0701873779297 17.0652526855469
17.0621643066406 19.6688018798828
18.0716796875 19.8661254882813
19.0795959472656 20.2305358886719
20.0910125732422 22.4033935546875
C:\>

最小自乗法で必要な∑x_n*x_n、∑y_n*x_n、∑x_n、∑y_n、∑1を求めるために、それぞれのxについてx_n*x_n、y_n*x_nを求めるんだな。

C:\>perl -alne "push @F,$F[0]**2,$F[0]*$F[1]; print \"@F\";" a.dat
1.01757202148438 2.11520385742188 1.03545281890781 2.15237226504834
2.05708923339844 3.05591735839844 4.23161611416378 6.28629469611683
3.03453979492188 4.1234130859375 9.20843176696452 12.512661100179
4.03664245605469 5.36311340332031 16.2944823180232 21.6489712604787
5.02697448730469 6.1601806640625 25.2704724960123 30.9670710354299
6.01442565917969 7.09348449707031 36.173316009799 42.663235172173
7.04530944824219 8.60143432617187 49.6363852214907 60.5997665266134
8.09518127441406 9.57536010742187 65.531959865624 77.5142758373729
9.06552429199219 10.8480926513672 82.1837306887005 98.3436474527513
10.0904388427734 11.7131164550781 101.81695603975 118.190485248248
11.0158996582031 12.6625122070312 121.350045279599 139.488963893428
12.0860198974609 14.0559478759766 146.071876961821 169.880465706726
13.0393676757813 14.6789276123047 170.02510938421 191.403934223019
14.0699066162109 16.3958953857422 197.962272188895 230.688716966556
15.0349884033203 16.2184478759766 226.050876287976 243.844175735163
16.0701873779297 17.0652526855469 258.250922361771 274.241808308457
17.0621643066406 19.6688018798828 291.117450826801 335.592329389322
18.0716796875 19.8661254882813 326.5856067276 359.014256455899
19.0795959472656 20.2305358886719 364.030981510914 385.990450552516
20.0910125732422 22.4033935546875 403.648786218176 450.10686159052
C:\>

確認できたらこいつをパイプ処理でb.datに出力してするんだな。

C:\>perl -alne "push @F,$F[0]**2,$F[0]*$F[1]; print \"@F\";" a.dat>b.dat

次に各行の和とデータの数を出力するんだな。

C:\>perl -alne "$x+=$F[0];$y+=$F[1];$xx+=$F[2];$xy+=$F[3];$n=$.;END{print \"$x $y $xx $xy $n\";}" b.dat>c.dat
C:\>perl -pe "" c.dat
211.10451965332 241.895156860352 2896.4767310872 3251.13074341602 20 
C:\>

最後に、下のように行列計算を行って結果を出力するんだな。

[∑y_n*x_n]=[∑x_n*x_n ∑x_n][a]
[∑y_n    ] [∑x_n     ∑1  ][b]
C:\WINDOWS\デスクトップ>perl -alne "$d=$F[2]*$F[4]-$F[0]**2;@a=(($F[4]*$F[3]-$F[0]*$F[1])/$d,($F[2]*$F[1]-$F[0]*$F[3])/$d);END{print \"@a\";}" c.dat
1.04437437081074 1.07115034860561
C:\WINDOWS\デスクトップ>

1つ目がa、2つ目がbなんだな。自信がないのでgnuplotで確かめてみるんだな。

gnuplot> fit a*x+b 'a.dat' via a,b
Final set of parameters            Asymptotic Standard Error
=======================            ==========================
a               = 1.04437          +/- 0.01597      (1.529%)
b               = 1.07115          +/- 0.1922       (17.94%)
correlation matrix of the fit parameters:
               a      b
a               1.000
b              -0.877  1.000
gnuplot> print a
1.04437437082676
gnuplot> print b
1.07115034838608
gnuplot> plot 'a.dat',a*x+b

確かに近い値となっていることがわかるんだな。完成したのでGIFで出力しておくんだな。

gnuplot> set terminal gif
Terminal type set to 'gif'
Options are 'small size 640,480 '
gnuplot> set output 'a.gif'
gnuplot> plot 'a.dat',a*x+b
gnuplot>

ソーシャルブックマーク

  1. はてなブックマーク
  2. Google Bookmarks
  3. del.icio.us

ChangeLog

  1. Posted: 2007-11-13T13:35:07+09:00
  2. Modified: 2007-11-13T04:43:05+09:00
  3. Generated: 2017-09-26T23:09:18+09:00