特にお世話になっている最適化計算についてちょっとご紹介．
最適化で使っているパッケージは, scipy.
他に Numericが必要.
excel等で近似多項式を求める事が出来るけど，奇数次数のみの多項式に
フィットさせるなどは不可能．
そういうときにPythonの scipyは便利．
次のスクリプトは９次項までの奇数次数のみの多項式へのフィッティングを行うもの．
私は魚眼カメラの射影式を扱うとき重宝している．
01 | import scipy |
02 | import scipy.optimize |
03 | import Numeric as n |
04 |
05 | def func(p, x): |
06 | k1,k2,k3,k4,k5 = p |
07 | p0=[k5,0.,k4,0.,k3,0.,k2,0.,k1,0.] |
08 | return scipy.polyval(p0, x) |
09 |
10 | def residue(p,y,x,sigma): |
11 | #print "residual_func ", y, x |
12 | err = y - func(p, x) |
13 | #print "error ", err |
14 | return err |
15 |
16 | def fitting(x, y): |
17 | sigma = 0 |
18 | r = scipy.optimize.leastsq( |
19 | residue, |
20 | [1.0,0,0,0,0],args=(y, x, sigma), full_output=1, ftol=1e-32, xtol=1e-32, maxfev=10000) |
21 | return r[0].tolist() |
22 |
23 | if __name__== '__main__': |
24 | sample = n.array([ |
25 | 0.16581,0.128181366, |
26 | 0.33161,0.258930716, |
27 | 0.49742,0.393306338, |
28 | 0.66323,0.530482793, |
29 | 0.82903,0.666571512, |
30 | 0.99484,0.793397256, |
31 | 1.16064,0.899001855, |
32 | 1.32645,0.970166103, |
33 | 1.49226,0.99910374, |
34 | 1.65806,0.991575672]) |
35 |
36 | xy = n.reshape(sample, (10,2)) |
37 | x = xy[:,0] |
38 | y = xy[:,1] |
39 | optimized_parameter = fitting(x, y) |
40 |
41 | print optimized_parameter |
scipy.optimize.leastsqは minpackを移植したもののようである．
Levenberg-Marquardt法を用いた最適化をお手軽に行うにはもってこいだと思う．
