とりあえずデータファイルを次のようにして生成したんだな。
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>