SciPyリファレンス scipy.optimize 日本語訳にいろいろな最適化の関数が書いてあったので、いくつか試してみた。
y = c + a*(x - b)**2の2次関数にガウスノイズを乗せて、これを2次関数で最適化してパラメータ求めてみた。

import numpy as np
from matplotlib import pyplot as plt
import scipy.optimize as optimize
%matplotlib inline

この後で決定係数のR2や、平均2乗誤差、平均絶対誤差を計算するためにScikit-learnもインポートする。

from sklearn.metrics import r2_score
from sklearn.metrics import mean_squared_error
from sklearn.metrics import mean_absolute_error

モデルになる2次関数の式

def base_func(x,a,b,c):
    y = c + a*(x - b)**2
    return y

パラメータを入れて計算して、グラフ表示する。

x = np.arange(-30, 30, 1)
para = [2.0,5.0,10.0]
np.random.seed(seed=10)
y = base_func(x,para[0],para[1],para[2])+np.random.normal(0, 60, len(x))

plt.scatter(x , y)
plt.show()

最適化したグラフを書くための関数を用意。

def plot_fig(x,y,para):
    fit = base_func(x,para[0],para[1],para[2])
    r2 = r2_score(y, fit)
    print()
    print('R2:',r2)
    print('Results:', para)
    plt.scatter(x,y,label='data')
    plt.plot(x,fit,'r',label='fit')
    plt.legend()
    plt.show()

最小二乗法（RMSEを最小化する）

Python SciPy : 非線形最小二乗問題の最適化アルゴリズム

# カーブフィッティング
cf = optimize.curve_fit(f=base_func, xdata=x, ydata=y, p0=para)
print(cf)
print(cf[0])
plot_fig(x,y,cf[0])


(array([ 2.01296897,  4.92997267, 16.69809461]), array([[ 6.97990322e-04, -1.88282482e-03, -1.88758970e-01],
[-1.88282482e-03,  1.54028062e-02,  2.83488209e-01],
[-1.88758970e-01,  2.83488209e-01,  1.06165700e+02]]))
[ 2.01296897  4.92997267 16.69809461]

R2: 0.9934619969757706
Results: [ 2.01296897 4.92997267 16.69809461]

非線形最小二乗法

残差関数を定義する必要がある。

# residual function
def ls_res(param,x,y):
    residual = y - base_func(x,param[0],param[1],param[2])
    return residual

# レガシー関数 leastsq(func, x0[, args, Dfun, full_output, …]) 方程式セットの平方和を最小化する。
le_lsq=optimize.leastsq(ls_res, para, args=(x, y))
print(le_lsq)
print(le_lsq[0])
plot_fig(x,y,le_lsq[0])


(array([ 2.01296897,  4.92997267, 16.69809461]), 1)
[ 2.01296897  4.92997267 16.69809461]

R2: 0.9934619969757706
Results: [ 2.01296897 4.92997267 16.69809461]

# 推奨　least_squares(fun, x0[, jac, bounds, …]) 変数境界付きの非線形最小二乗法問題を解く。
lsq=optimize.leastsq(ls_res, para, args=(x, y))
print(lsq)
print(lsq[0])
plot_fig(x,y,lsq[0])

(array([ 2.01296897,  4.92997267, 16.69809461]), 1)
[ 2.01296897  4.92997267 16.69809461]

R2: 0.9934619969757706
Results: [ 2.01296897 4.92997267 16.69809461]

最適化

最小化する関数を定義する
Root Mean Squared Error (RMSE)
$$RMSE = \sqrt{ \frac{ \sum_{i} (y_{\mathrm{obs}, i} - y_{\mathrm{pred}, i})^2 }{ n }}$$

# RMSEの関数を定義
def rmse_res(param,x,y):
    residual = y - base_func(x,param[0],param[1],param[2])
    resid_=np.sqrt(((residual*residual).sum())/len(x))
    return resid_

# 　参考：Scikitlearnの関数を使った場合。rmse_resと同じ結果になる
# def rmse_res2(param,x,y):
#    fit = base_func(x,param[0],param[1],param[2])
#    rmse = np.sqrt(mean_squared_error(y, fit))
#    return rmse

# 局所（多変数）最適化 minimize(fun, x0[, args, method, jac, hess, …]) １つ以上の変数のスカラー関数の最小化
mi_rmse=optimize.minimize(rmse_res, para, args=(x, y),method='Nelder-Mead')
print(mi_rmse)

# 　結果の出力に注意
print(mi_rmse.x)
plot_fig(x,y,mi_rmse.x)

     final_simplex: (array([[ 2.01296895,  4.92997307, 16.69812532],
           [ 2.01296886,  4.92997279, 16.69816863],
           [ 2.01296906,  4.92997222, 16.6981314 ],
           [ 2.01296907,  4.92997182, 16.69803887]]), array([53.48433856, 53.48433856, 53.48433856, 53.48433856]))
               fun: 53.48433855800152
           message: 'Optimization terminated successfully.'
              nfev: 173
               nit: 93
            status: 0
           success: True
                 x: array([ 2.01296895,  4.92997307, 16.69812532])
    [ 2.01296895  4.92997307 16.69812532]

R2: 0.9934619969757669
Results: [ 2.01296895 4.92997307 16.69812532]

#  fmin(func, x0[, args, xtol, ftol, maxiter, …])　レガシー関数 非推奨
fmin_rmse=optimize.fmin(rmse_res, para, args=(x, y))
print(fmin_rmse)
plot_fig(x,y,fmin_rmse)

    Optimization terminated successfully.
             Current function value: 53.484339
             Iterations: 93
             Function evaluations: 173
    [ 2.01296895  4.92997307 16.69812532]

R2: 0.9934619969757669
Results: [ 2.01296895 4.92997307 16.69812532]

大域最適化

# basinhopping(func, x0[, niter, T, stepsize, …]) basin-hoppingアルゴリズムを使って関数の大域最小値を求める。
bh_rmse=optimize.basinhopping(rmse_res, para, minimizer_kwargs={"method": "L-BFGS-B","args":(x, y)})
print(bh_rmse)
print(bh_rmse.x)
plot_fig(x,y,bh_rmse.x)

     fun: 53.48433855798703
     lowest_optimization_result:       fun: 53.48433855798703
     hess_inv: <3x3 LbfgsInvHessProduct with dtype=float64>
          jac: array([9.66338121e-05, 5.68434189e-06, 0.00000000e+00])
      message: b'CONVERGENCE: REL_REDUCTION_OF_F_<=_FACTR*EPSMCH'
         nfev: 120
          nit: 24
       status: 0
      success: True
            x: array([ 2.012969  ,  4.92997261, 16.6980913 ])
                        message: ['requested number of basinhopping iterations completed successfully']
          minimization_failures: 0
                           nfev: 6564
                            nit: 100
                              x: array([ 2.012969  ,  4.92997261, 16.6980913 ])
    [ 2.012969    4.92997261 16.6980913 ]

R2: 0.9934619969757704
Results: [ 2.012969 4.92997261 16.6980913 ]

# differential_evolution(func, bounds[, args, …]) 多変数関数の大域最小値を求める。
bounds = [(-0,10), (-10,10), (-50,50)]
de_rmse=optimize.differential_evolution(rmse_res, bounds, args=(x, y))
print(de_rmse)
print(de_rmse.x)
plot_fig(x,y,de_rmse.x)


         fun: 53.48433860476771
         jac: array([-1.77564630e-03, -5.75539616e-04, -2.98427949e-05])
     message: 'Optimization terminated successfully.'
        nfev: 935
         nit: 18
     success: True
           x: array([ 2.01297481,  4.92995775, 16.6948975 ])
    [ 2.01297481  4.92995775 16.6948975 ]

R2: 0.9934619969643333
Results: [ 2.01297481 4.92995775 16.6948975 ]

# brute(func, ranges[, args, Ns, full_output, …]) 与えられた範囲でブルートフォース法によって関数を最小化する。
br_rmse=optimize.brute(rmse_res, ranges=bounds, args=(x, y))
print(br_rmse)
plot_fig(x,y,br_rmse)

    [ 2.01296928  4.92997168 16.69807946]

R2: 0.9934619969757458
Results: [ 2.01296928 4.92997168 16.69807946]

最小化する関数を絶対誤差に変更

Mean Absolute Error (MAE)
$$ MAE = \frac{ \sum_{i} | y_{\mathrm{obs}, i} - y_{\mathrm{pred}, i} | }{ n }$$

# MAEの関数を定義
def mae_res(param,x,y):
    residual = np.abs(y - base_func(x,param[0],param[1],param[2]))
    resid_=((residual).sum())/len(x)
    return resid_

# 参考：Scikitlearnの関数を使った場合。mae_resと同じ結果になる
# def mae_res2(param,x,y):
#    fit = base_func(x,param[0],param[1],param[2])
#    mae = mean_absolute_error(y,fit)
#    return  mae

mi_mae=optimize.minimize(mae_res, para, args=(x, y),method='Nelder-Mead')
print(mi_mae)
print(mi_mae.x)
plot_fig(x,y,mi_mae.x)


     final_simplex: (array([[ 2.00441621,  5.08263527, 12.59609788],
           [ 2.00441602,  5.0826368 , 12.59609761],
           [ 2.00441614,  5.08263639, 12.59604168],
           [ 2.00441602,  5.0826364 , 12.5961183 ]]), array([41.23285905, 41.23285928, 41.23285957, 41.2328597 ]))
               fun: 41.23285904630758
           message: 'Optimization terminated successfully.'
              nfev: 197
               nit: 108
            status: 0
           success: True
                 x: array([ 2.00441621,  5.08263527, 12.59609788])
    [ 2.00441621  5.08263527 12.59609788]

R2: 0.9932364583092926
Results: [ 2.00441621 5.08263527 12.59609788]

fmin_mae=optimize.fmin(mae_res, para, args=(x, y))
print(fmin_mae)
plot_fig(x,y,fmin_mae)

    Optimization terminated successfully.
             Current function value: 41.232859
             Iterations: 108
             Function evaluations: 197
    [ 2.00441621  5.08263527 12.59609788]

R2: 0.9932364583092926
Results: [ 2.00441621 5.08263527 12.59609788]

# differential_evolution(func, bounds[, args, …]) 多変数関数の大域最小値を求める。
de_mae=optimize.differential_evolution(mae_res, bounds, args=(x, y))
print(de_mae)
print(de_mae.x)
plot_fig(x,y,de_mae.x)

         fun: 41.22302370179955
         jac: array([ 8.8972655 , -2.14848157, -0.03333298])
     message: 'Optimization terminated successfully.'
        nfev: 730
         nit: 13
     success: True
           x: array([ 2.0129207 ,  5.01018698, 13.87485399])
    [ 2.0129207   5.01018698 13.87485399]

R2: 0.993388114525508
Results: [ 2.0129207 5.01018698 13.87485399]

最小二乗法 （RMSEを最小化する）

非線形最小二乗法

最適化

大域最適化

最小化する関数を絶対誤差に変更

最小二乗法（RMSEを最小化する）