読者です 読者をやめる 読者になる 読者になる

pythonで非線形最小二乗法

scipy.optimize.leastsq()の使い方で苦戦したのでメモ

LM(Levenberg-Marquardt)法を実装するのが面倒だったので既存のleastsq(LM法で非線形最小二乗してくれるやつ)を使うことにした.
というわけで、良くわからないまま適当に以下のようなコードを書いた(研究用のコードを載せるのはアレなので適当に書き換えてる)

"""  
minimize y - a x^2 + b x + c 
params : a,b,c
array_of_Point2d : (x,y) data
"""
def cost_func(param,array_of_Point2d):
    a,b,c = params
    sum_of_err = 0.0
    for point in array_of_Point2d:
        x,y = point.x,point.y
        sum_of_err += (y - (a*x**2 + b*x + c))**2
    return cost_func

params = [0,0,0]
array_of_Point2d = [Point2d(i,i) for i in range(10)]
leastsq(cost_func,x0=params,args=array_of_Point2d)

で、実行してみると案の定エラーを吐かれた

Improper input: N=6 must not exceed M=2

ググる
→ http://stackoverflow.com/questions/21592348/scipy-leastsq-fitting-a-square-grid-to-experimental-points-in-2d

The function passed to leastsq (e.g. err_func) should return an array of values of the same shape as xeand ye

ようするにleastsqに渡す関数の戻り値は、関数に渡す引数と同じ形じゃないとアカンとのこと.
x,yデータの数だけ誤差を計算し、その誤差を全て返すのが正しい形っぽい

修正版

"""
minimize y - a x^2 + b x + c 
params : a,b,c
array_of_Point2d : (x,y) data
"""
def cost_func(param,array_of_Point2d):
    a,b,c = params
    errs = []
    for point in array_of_Point2d:
        x,y = point.x,point.y
        errs.append(y - (a*x**2 + b*x + c))
    return errs

params = [0,0,0]
array_of_Point2d = [Point2d(i,i) for i in range(10)]
leastsq(cost_func,x0=params,args=array_of_Point2d

というか、僕がミスったようなことは全部リファレンスに書いてあったscipy.optimize.leastsq - SciPy v0.18.1 Reference
英語だからって面倒くさがらずにちゃんと読もうな!

追記

どうやら使っているScipyが古かった模様 scipy 0.17.1以降はleast_squares()という関数があって、こいつを使うよう推奨されている
least squares - SciPy: leastsq vs least_squares - Stack Overflow

least_squares()ではキーワード引数methodを指定することで以下の三つのアルゴリズムを使い分けることができる
- ‘trf’ : Trust Region Reflective algorithm
- ‘dogbox’ : dogleg algorithm
- ‘lm’ : Levenberg-Marquardt

ついでに言えばleast_squaresの方は解に上下限を付けることができ、解(x,y,z)について
bounds = ([x_l_lim,y_l_lim,z_l_lim],[x_u_lim,y_u_lim,z_u_lim])
のように解の下限リストと上限リストをタプルにして渡せばよい
(注:LM法を使う場合は設定できない)