2010年3月24日水曜日

Pythonで非線形最適化計算する

私自身がPythonを便利に思って重点的に利用させてもらっている技術分野は,科学技術計算とUSBやシリアルなどとの通信に関するもの.

特にお世話になっている最適化計算についてちょっとご紹介.

最適化で使っているパッケージは, scipy.
他に Numericが必要.

excel等で近似多項式を求める事が出来るけど,奇数次数のみの多項式に
フィットさせるなどは不可能.

そういうときにPythonの scipyは便利.
次のスクリプトは9次項までの奇数次数のみの多項式へのフィッティングを行うもの.
私は魚眼カメラの射影式を扱うとき重宝している.

import scipy
import scipy.optimize
import Numeric as n

def func(p, x):
k1,k2,k3,k4,k5 = p
p0=[k5,0.,k4,0.,k3,0.,k2,0.,k1,0.]
return scipy.polyval(p0, x)

def residue(p,y,x,sigma):
#print "residual_func ", y, x
err = y - func(p, x)
#print "error ", err
return err

def fitting(x, y):
sigma = 0
r = scipy.optimize.leastsq(
residue,
[1.0,0,0,0,0],args=(y, x, sigma), full_output=1, ftol=1e-32, xtol=1e-32, maxfev=10000)
return r[0].tolist()

if __name__== '__main__':
sample = n.array([
0.16581,0.128181366,
0.33161,0.258930716,
0.49742,0.393306338,
0.66323,0.530482793,
0.82903,0.666571512,
0.99484,0.793397256,
1.16064,0.899001855,
1.32645,0.970166103,
1.49226,0.99910374,
1.65806,0.991575672])

xy = n.reshape(sample, (10,2))
x = xy[:,0]
y = xy[:,1]
optimized_parameter = fitting(x, y)

print optimized_parameter

scipy.optimize.leastsqは minpackを移植したもののようである.
Levenberg-Marquardt法を用いた最適化をお手軽に行うにはもってこいだと思う.