GSL Shell 学習:Fitting サブ

2017-01-09 :  PCクリニック
Python、C言語、Perl、グルコサミン、Firefox
前(2017-01-08)の記事「Perl から Lua : 関数Fitting ( By GSL )」では、
  ・・・・・
  ・・・・・
  はやり、関数フィッティングは GSL Shell
  と云うことですネ。
  ・・・・・
  ・・・・・

そうすると、もう gnuplot でのフィッティングは忘れる?


それで、
前回の様にたくさんの行をコーディングしたくない。


なので、
例えば、昨年2月(2016-02-24)の記事:
GSL Shell 学習:Nonlinear Least Squares Fit」のコード例なら、
X = matrix.vec{ 5, 7, 8, 9, 10, 13, 15, 17, 18 }
Y = matrix.vec{ 0.67, 0.84, 1.29, 1.25, 0.69, 0.35, 0.18, 0.09, 0.02 }
------- Y ~ a * X^3 + b * X^2 + c * X + d で近似
-------------------------------------------------
require'Fit_SUB'
s = Fit_GSL( X, Y, 'k[1]*x^3+k[2]*x^2+k[3]*x+k[4]', '{ 1, 0, 0, 0 }', 'x^3', 'x^2', 'x', '1')
-------------------------------------------------
print( s.x, s.chisq ) -- 求めた係数、カイ二乗
の様な、実質1行で済ませたい!

つまり、
 'Fit_GSL'関数のみを定義している、
 'Fit_SUB' サブプログラムを作る。

'Fit_GSL'関数の仕様
引数は、5つ以上で、
第1、2は、観測点データベクトル X、Y
第3は、近似関数(文字列)
第4は、初期推定値テーブル
第5以降は、ヤコビアン(文字列)

但し、ベクトル X、Y は、matrix.vec で確保。
近似関数とヤコビアンでの、係数は k[] で表記。
初期推定値テーブルの要素数とヤコビアンの数は、ちょうど係数の数で。


なお、ヤコビアンを求めたいときは、
Maxima」か「Euler」を使う。


結果、出来上がった'Fit_SUB' サブプログラム:
function Fit_GSL( x_, y_, pp1, pp2, ... )
local pj = {...}
local n, k = #x_, #pj
local X, Y = x_, y_
local v = {}
v[1]='function fdf( k, f, J )\n'
v[2]=' for i=1, %d do\n local x, y = X[i], Y[i]\n' -- n
v[3]=' if f then f[i] = ( %s - y ) end\n' -- pp1
v[4]=' if J then\n'
local vj=' J:set( i, %d, %s )\n' -- (pj)
v[5]=' end; end; end\n'
v[6]='s_ = num.nlinfit{ n=%d, p=%d }\n' -- n, k
v[7]='s_:set( fdf, matrix.vec%s )\n' -- pp2

local vv = v[1]..string.format( v[2], n )
vv = vv..string.format( v[3], pp1 )..v[4]
for q,pjq in ipairs(pj) do
vv = vv..string.format( vj, q, pjq )
end
vv = vv..v[5]..string.format( v[6], n, k )
vv = vv..string.format( v[7], pp2 )

assert( load( vv ) )() ----- function fdf 定義

for i=1, 100 do
s_:iterate()
if s_:test( 0, 1e-8 ) then break end
end
return s_
end
でしょうか?


なお、前の記事「Perl から Lua : 関数Fitting ( By GSL )」のコードなら、
X = matrix.vec{ ・・・・・ }
Y = matrix.vec{ ・・・・・ }
------- Y ~ exp( -b * log(X) + c ) で近似
-------------------------------------------------
exp=math.exp; log=math.log
require'Fit_SUB'
s = Fit_GSL( X, Y, 'exp(-k[1]*log(x)+k[2])', '{ 1, 0 }', '-log(x)*exp(-k[1]*log(x)+k[2])', 'exp(-k[1]*log(x)+k[2])')
-------------------------------------------------
b, c = s.x[1], s.x[2] -- 求めた係数
の様なものになる。


本日はここまで。


“Perl から Lua ( GSL Shell ) への移植”準備その2は終了?


見ていただいた序でとは厚かましい限りですが、
お帰りに投票して頂けるとなお嬉しいです。 ⇒ blogram投票ボタン


161210
関連記事
スポンサーサイト

コメントの投稿

管理者にだけ表示を許可する

おきてがみ/blogram
blogram投票ボタン



おきてがみ

最新記事
カレンダー
03 | 2017/04 | 05
- - - - - - 1
2 3 4 5 6 7 8
9 10 11 12 13 14 15
16 17 18 19 20 21 22
23 24 25 26 27 28 29
30 - - - - - -
月別アーカイブ
カテゴリ
最新コメント
検索フォーム
リンク
プロフィール

<紙>

Author:<紙>
ようこそ。
「パソコンヲタクの雑記帳」
もろもろなことを綴っています。
パソコン ヲタクってねくら?
画像は kami でなく kani です。

カウンター(fc2、i2i) /Google Analytics


i2i(from 2010-08-24)
Total =
Today  =  
Yesterday=
アンチエイジング

Google Analytics
ブックマーク