9
16

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

More than 3 years have passed since last update.

pythonでパラメータの制限付きでfittingを行う方法

Last updated at Posted at 2019-09-10

パラメータの制限付きfitting

python で直線フィットとガウシアンフィットをする簡単な方法 では、[scipy.optimize.leastsq] (https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.leastsq.html) を用いた簡単なフィット方法を紹介したが、パラメータの制限付きでフィットしたい場合は、[scipy.optimize.least_squires] (https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.least_squares.html) を用いる必要がある。基本的な使い方は同じで、bounds に上限値と下限値を入れるだけである。

制限方法

scipy.optimize の least_squares には bounds というオプションがあるので、ここにパラメータの上限と下限の配列を入れるだけである。least_squares は2乗値の最小化問題なので、一般的な最小化問題よりは速く収束する。

voigt 関数のパラメータ制限付き fitting の例

voigt 関数のフィットの例を使って、パラメータの制限付きフィットの例を示す。
エラーの計算方法が[scipy.optimize.leastsq] (https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.leastsq.html) と[scipy.optimize.least_squires] (https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.least_squares.html) で若干異なる。

python で voigt 関数をプロット方法は、pythonでフォークト関数(Voigt function)をプロットする方法を参照。

この例では、bounds = [[0, 0, 0, 0],[300, 100, 100, 100]] として、4つのパラメータ(norm,center,lw,gw)を適当に正の範囲で制限する例を示した。

fit_voigt_lq.py

#!/usr/bin/env python

__version__= '1.0'
# This is created from the code made by M.Sawada during SXS tests. 

import matplotlib.pyplot as plt
plt.rcParams['font.family'] = 'serif'
import numpy as np
import scipy as sp
import scipy.optimize as so
import scipy.special

def calcchi(params,consts,model_func,xvalues,yvalues,yerrors):
    model = model_func(xvalues,params,consts)
    chi = (yvalues - model) / yerrors
    return(chi)

# optimizer
def solve_least_squares(xvalues, yvalues, yerrors, param_init, consts, model_func, bounds=(-np.inf, np.inf)):
    param_output = so.least_squares(
        calcchi,
        param_init,
        bounds=bounds,
        args=(consts, model_func, xvalues, yvalues, yerrors))
    param_result = param_output.x
    hessian = np.dot(param_output.jac.T, param_output.jac)
    covar_output = np.array(np.linalg.inv(hessian))
    error_result = np.sqrt(covar_output.diagonal())
    dof = len(xvalues) - 1 - len(param_init)
    chi2 = np.sum(np.power(calcchi(param_result,consts,model_func,xvalues,yvalues,yerrors),2.0))
    return([param_result, error_result, chi2, dof])

# model
def mymodel(xvalues,params,consts):
    norm,center,lw,gw = params
    # norm : normalization 
    # center : center of Lorentzian line
    # lw : HWFM of Lorentzian 
    # gw : sigma of the gaussian 
    c0,c1 = consts
    z = (xvalues - center + 1j*lw)/(gw * np.sqrt(2.0))
    w = scipy.special.wofz(z)
    model_y = norm * (w.real)/(gw * np.sqrt(2.0*np.pi))
    return model_y

# set data
x = np.array([ 1, 2, 3,  4,  5, 6,   7,  8,  9, 10, 11, 12, 13, 14, 15])
y = np.array([ 1, 2, 5, 10, 28, 32, 37, 51, 41, 25,  8,  3,  1,  2,  1])

# initialize parameter 
norm_init = 60.
lorent_init = 2.
gauss_sigma_init = 2.
init_params = np.array([norm_init, np.median(x), lorent_init, gauss_sigma_init])
# initialize constants 
consts = np.array([0,0]) 
bounds = [[0, 0, 0, 0],
         [300, 100, 100, 100]]

# do fit 
result, error, chi2, dof = solve_least_squares(x, y, np.sqrt(y), init_params, consts, mymodel, bounds)

# get results 
ene = result[1]
lw = np.abs(result[2])
gw = np.abs(result[3])

enee = error[1]
lwe = np.abs(error[2])
gwe = np.abs(error[3])

fwhm = 2.35 * gw
fwhme = 2.35 * gwe

label = "E = " + str("%4.2f(+/-%4.2f)" % (ene,enee)) + " FWHM = " + str("%4.2f(+/-%4.2f)" % (fwhm,fwhme) + " (FWHM)")
print(label)

# plot results
plt.title("Voigt fit : " + label)
model_y = mymodel(x,result,consts)
plt.errorbar(x, y, yerr = np.sqrt(y), fmt="bo", label = "data")
plt.plot(x, model_y, 'r-', label = "model")
plt.legend(numpoints=1, frameon=False, loc="best")
plt.grid(linestyle='dotted',alpha=0.5)

plt.savefig("voigtfit_lq.png")
plt.show()

生成される図はこちら。

voigtfit_lq.png

一般的な注意

パラメータに制限をつけると、計算が遅くなったり、変な値に収束することも多い。本当に制限すべきパラメータなのか、データからそもそも制限できないパラメータをフリーにしてないか、あるいは何らかの方法でパラメータを固定できないのかをよく考えた上で、必要な場合のみ制限をつけよう。

9
16
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
9
16

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?