Help us understand the problem. What is going on with this article?

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

パラメータの制限付きfitting

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

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

voigt 関数のフィットの例を使って、パラメータの制限付きフィットの例を示す。
エラーの計算方法がscipy.optimize.leastsqscipy.optimize.least_squires で若干異なる。

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

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],
         [100, 100, 100, 100]]

# do fit 

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

# 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

一般的な注意

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

Why do not you register as a user and use Qiita more conveniently?
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away
Comments
Sign up for free and join this conversation.
If you already have a Qiita account
Why do not you register as a user and use Qiita more conveniently?
You need to log in to use this function. Qiita can be used more conveniently after logging in.
You seem to be reading articles frequently this month. Qiita can be used more conveniently after logging in.
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away