積分関数を含む式を使って2つのパラメータでフィッティング

Question

前提・実現したいこと

pythonで積分関数（フォークト関数）を使ってフィッティングしようとしています。
単純な関数ではフィッティングができたのですが、積分を組み込むとうまくいきません。
フォークト関数では、zについて積分をします。
パラメータは2つでparam1とparam2です。

発生している問題・エラーメッセージ

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-10-baa10670dc7d> in <module>
     21 
     22 p0 = 0.5, 0.5
---> 23 popt, pcov = curve_fit(func, array_x, array_y, p0)

~\Anaconda3\lib\site-packages\scipy\optimize\minpack.py in curve_fit(f, xdata, ydata, p0, sigma, absolute_sigma, check_finite, bounds, method, jac, **kwargs)
    750         # Remove full_output from kwargs, otherwise we're passing it in twice.
    751         return_full = kwargs.pop('full_output', False)
--> 752         res = leastsq(func, p0, Dfun=jac, full_output=1, **kwargs)
    753         popt, pcov, infodict, errmsg, ier = res
    754         cost = np.sum(infodict['fvec'] ** 2)

~\Anaconda3\lib\site-packages\scipy\optimize\minpack.py in leastsq(func, x0, args, Dfun, full_output, col_deriv, ftol, xtol, gtol, maxfev, epsfcn, factor, diag)
    381     if not isinstance(args, tuple):
    382         args = (args,)
--> 383     shape, dtype = _check_func('leastsq', 'func', func, x0, args, n)
    384     m = shape[0]
    385 

~\Anaconda3\lib\site-packages\scipy\optimize\minpack.py in _check_func(checker, argname, thefunc, x0, args, numinputs, output_shape)
     24 def _check_func(checker, argname, thefunc, x0, args, numinputs,
     25                 output_shape=None):
---> 26     res = atleast_1d(thefunc(*((x0[:numinputs],) + args)))
     27     if (output_shape is not None) and (shape(res) != output_shape):
     28         if (output_shape[0] != 1):

~\Anaconda3\lib\site-packages\scipy\optimize\minpack.py in func_wrapped(params)
    456     if transform is None:
    457         def func_wrapped(params):
--> 458             return func(xdata, *params) - ydata
    459     elif transform.ndim == 1:
    460         def func_wrapped(params):

<ipython-input-10-baa10670dc7d> in func(x, param1, param2)
     17     def voigt(z):
     18         return param1 * 0.8256 / np.pi**(3/2) / param2 / widthL * np.exp((-1) * (2 * z / param2)**2 * np.log(2)) / (1 + (2 * (x - z) / widthL)**2)
---> 19     integration = integrate.quad(voigt, -np.inf, np.inf)[0]
     20     return integration
     21 

~\Anaconda3\lib\site-packages\scipy\integrate\quadpack.py in quad(func, a, b, args, full_output, epsabs, epsrel, limit, points, weight, wvar, wopts, maxp1, limlst)
    339     if weight is None:
    340         retval = _quad(func, a, b, args, full_output, epsabs, epsrel, limit,
--> 341                        points)
    342     else:
    343         retval = _quad_weight(func, a, b, args, full_output, epsabs, epsrel,

~\Anaconda3\lib\site-packages\scipy\integrate\quadpack.py in _quad(func, a, b, args, full_output, epsabs, epsrel, limit, points)
    448             return _quadpack._qagse(func,a,b,args,full_output,epsabs,epsrel,limit)
    449         else:
--> 450             return _quadpack._qagie(func,bound,infbounds,args,full_output,epsabs,epsrel,limit)
    451     else:
    452         if infbounds != 0:

TypeError: only size-1 arrays can be converted to Python scalars

該当のソースコード

import numpy as np
from scipy import integrate
from scipy.optimize import curve_fit
import math
import pandas as pd

Eins = 2.01 * 10**8
Lorentzbroad = Eins / 2 / math.pi
widthL = Lorentzbroad / 10**9

# フィッティング用データの宣言
df = pd.read_csv('book1.csv')
array_x=np.array(df['t'])
array_y=np.array(df['y'])

def func(x, param1, param2):
    def voigt(z):
        return param1 * 0.8256 / np.pi**(3/2) / param2 / widthL * np.exp((-1) * (2 * z / param2)**2 * np.log(2)) / (1 + (2 * (x - z) / widthL)**2)
    integration = integrate.quad(voigt, -np.inf, np.inf)[0]
    return integration

p0 = 0.5, 0.5
popt, pcov = curve_fit(func, array_x, array_y, p0)

試したこと

補足情報（FW/ツールのバージョンなど）

python3 Anaconda3

Comments

[meg_]

・コードは「コードの挿入」で記入しましょう。
・エラーメッセージは全文掲載しましょう。

[voigt]

2点修正いたしました。
宜しくお願いいたします。

[meg_]

エラーは「integration = integrate.quad(voigt, -np.inf, np.inf)[0]」で発生していますね。
未定義の変数(x)があるのとbook1.csvのデータが不明の為検証は出来ませんが、curve_fitを使う前の不具合ですね。

Answer 1

下記関数でnumpy配列を返しているのがエラーの原因ではないでしょうか？（要素数1のnumpy配列であればOKです）

def voigt(z):
        return param1 * 0.8256 / np.pi**(3/2) / param2 / widthL * np.exp((-1) * (2 * z / param2)**2 * np.log(2)) / (1 + (2 * (x - z) / widthL)**2)

---

【追記】

積分したい関数に複数引数を与えたいという意味であれば、下記サイト等が参考になるかと思います。
SciPyによる積分(scipy.integrate)①｜SciPy入門 #8

Comments

[voigt]

ご確認、ご回答いただきありがとうございます。
データには問題ないことは確認しております。
voigt関数はxの値毎にzで積分した値を得るために使っておりますが、
x毎の値を返すにはどうしたら良いのでしょうか？

[meg_]

param1とparam2が配列値ですよね？そのため計算結果も配列になっているかと思います。

[voigt]

理解力不足で申し訳ありませんが、どのように修正すべきなのでしょうか？
複数パラメータでのフィッティングのため、出力を確認しながらdef functionの引数は設定しました。
よろしくお願い致します。

Answer 2

積分関数を組み込んで計算したい人のために一応こちらもあげておきます。
numpy.vectrizeを使うことで計算することができました。
ただし、Faddeeva関数の方が格段に計算速度が速いので、そちらがおすすめです。

import math
import numpy as np
import scipy.optimize
import scipy.integrate
import matplotlib.pyplot as plt
import pandas as pd

df = pd.read_csv('data.csv')
array_x=np.array(df['t'])
array_y=np.array(df['y'])

def integrand(x, param):
    z, peak, gw, lw = param
    return peak / math.pi**(3/2) / gw / lw * math.exp((-1) * (2 * x / gw)**2 * np.log(2)) / (1 + (2 * (z - x) / lw)**2)

def curve(z, peak, gw, lw):
    model = scipy.integrate.quad(integrand, -np.inf, np.inf, [z, peak, gw, lw])
    return model[0]

vcurve = np.vectorize(curve, excluded=set([1]))

popt, pcov = scipy.optimize.curve_fit(vcurve, array_x, array_y, p0=[1, 1, 1])

Answer 3

Faddeeva関数を計算する方法に変更しました。

import numpy as np
import scipy.special
from scipy.optimize import curve_fit
import pandas as pd

df = pd.read_csv('book1.csv')
array_x=np.array(df['t'])
array_y=np.array(df['y'])

def voigt(xval,norm,center,lw,gw):
    z = (xval - 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

popt, pcov = curve_fit(voigt, array_x, array_y)