pythonのフィッテング、least_squaresを試すべし

データに対して関数のフィッテングをしていたんですが、データも関数も同じなのにcurve_fitではうまくいかなかったのがleast_squaresを使ったらうまくいきました。使っているアルゴリズムが違うのだと思いますが、理由はわかんないです。

フィッテングさせたいデータと関数

・データ
指数関数的に減少していくデータの漸近部分で次の関数でフィッテングしたい。プログラム内では$mm$と置いていて、後半の漸近部分のみをフィッテングしたいので$mm[150:]$に対してフィッテングしてます。

・関数

y=(1-c)\exp(-\frac{c}{a}t)

変数は$t$、パラメタは$a$と$c$の二つです。
値の範囲は$a>0$、$1>c>0$です。およそ$[a, c]=[1000, 0.9]$が正解です。これらを範囲と初期値に設定しました。

うまくいかなかったscipy.optimize.curve_fitの場合

import numpy as np
from scipy.optimize import curve_fit

def func_c1(x, a, c):
    return (1-c)*np.exp(-x*c/a)
     
popt, pcov = curve_fit(func_c1, mm_fit[150:1000], times[150:1000], p0=[1000,0.9], bounds=([0,0], [np.inf,1]))
    
xd = np.linspace(min(times), max(times), 1000)
fitting = func_c1(xd, popt[0],popt[1])

print("Optimal parameters:", popt)

結果を確認しましょう。

Optimal parameters: [4000.000901077    0.000000000]

プロットしてみます。オレンジの点に対してフィッテングを行なっていて、緑の点線がフィッティングされた曲線です。

うまくできてないですね。

うまくいったscipy.optimize.least_squaresの場合

import numpy as np
from scipy.optimize import least_squares

def func_c1(t, params):
    a, c = params
    return (1 - c) * np.exp(-t * c / a)
    
def residuals(params, x, y):
    return y - func_c1(x, params)
    
initial_params = [800, 0.8]
    
result = least_squares(residuals, initial_params, args=(times[150:], mm[150:]), bounds=([0, 0], [np.inf, 1]))

print("Optimal parameters:", result.x)

出力も良さげです。

Optimal parameters: [1351.5720327910    0.8665103412]

同様にプロットしてみます。

うまくいってます。

まとめ

pythonのscipy.optimizeのcurve_fitを用いたフィッティングがうまくできなかったらleast_squaresを試してみようというメモでした。

実装の違いとしては、curve_fitは当てはめ関数をそのまま使えるのに対して、least_squaresは差分を使うところが違うんで書き換えるときはそこだけ注意です。