測定データを積分関数でフィッティングしたいです。

Question

###前提・実現したいこと
scipy.optimizeのlastsqを用いたフィッティングについて質問です。
フィッティングする関数に、積分関数を使いたいと思ってます。
残渣を求める行にintegrate.quad関数を用いて、その引数にパラメーターを指定してフィッティングをしようと思いましたが、下記のエラーが出てつまづいています。

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

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-32-02de193ee868> in <module>()
     25 
     26 param0 = [1.2, 1, 10.5, -24.3]
---> 27 param = optimize.leastsq(fanction_gauss, param0, args=(x, y))

~/.pyenv/versions/anaconda3-2.4.1/lib/python3.5/site-packages/scipy/optimize/minpack.py in leastsq(func, x0, args, Dfun, full_output, col_deriv, ftol, xtol, gtol, maxfev, epsfcn, factor, diag)
    375     if not isinstance(args, tuple):
    376         args = (args,)
--> 377     shape, dtype = _check_func('leastsq', 'func', func, x0, args, n)
    378     m = shape[0]
    379     if n > m:

~/.pyenv/versions/anaconda3-2.4.1/lib/python3.5/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):

<ipython-input-32-02de193ee868> in fanction_gauss(param, x, y)
     21 
     22 def fanction_gauss(param, x, y):
---> 23     residual = y - integrate.quad((param[0]*np.exp(-(x-param[2])**2/(2*param[1]))+param[3]), x-0.25, x+0.25)
     24     return residual
     25 

~/.pyenv/versions/anaconda3-2.4.1/lib/python3.5/site-packages/scipy/integrate/quadpack.py in quad(func, a, b, args, full_output, epsabs, epsrel, limit, points, weight, wvar, wopts, maxp1, limlst)
    321     if (weight is None):
    322         retval = _quad(func, a, b, args, full_output, epsabs, epsrel, limit,
--> 323                        points)
    324     else:
    325         retval = _quad_weight(func, a, b, args, full_output, epsabs, epsrel,

~/.pyenv/versions/anaconda3-2.4.1/lib/python3.5/site-packages/scipy/integrate/quadpack.py in _quad(func, a, b, args, full_output, epsabs, epsrel, limit, points)
    370 def _quad(func,a,b,args,full_output,epsabs,epsrel,limit,points):
    371     infbounds = 0
--> 372     if (b != Inf and a != -Inf):
    373         pass   # standard integration
    374     elif (b == Inf and a != -Inf):

ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()

###該当のソースコード

import sys, csv
import numpy as np
from scipy import integrate
sys.path.append("")
import matplotlib.pyplot as plt

cnt = 0
x = np.array([])
y = np.array([])
da_int = np.array([])

f = open("beamsize_y_up_45mm.csv", "r", encoding="utf_8_sig")    
try:
    reader = csv.reader(f)
    for row in reader:
        x = np.append(x, float(row[0]))
        y = np.append(y, float(row[1]))
        cnt = cnt + 1
finally:
    f.close()
    
def fanction_gauss(param, x, y):
    residual = y - integrate.quad((param[0]*np.exp(-(x-param[2])**2/(2*param[1]))+param[3]), x-0.25, x+0.25)
    return residual

param0 = [1.2, 1, 10.5, -24.3]
param = optimize.leastsq(fanction_gauss, param0, args=(x, y))

###試したこと
フィッティング関数を単なるgaussianにするとうまくいきました。

def gaussian(param, x, y):
    residual = y - (((param[0]/(np.sqrt(2*np.pi)*param[1]))*np.exp(-(x-param[2])**2/(2*param[1]**2)))+param[3])
    return residual

回答をいただいた後、次のようにfanction_gaussを変更しました。

def fanction_gauss(param, x, y):
    def f(x):
        return param[0]*np.exp(-(x-param[2])**2/(2*param[1]))+param[3]
    residual = y - integrate.quad(f, x-0.25, x+0.25)[0]
    return residual

array型のyからfloatを引いてしまっている？と、書きながら思いました。

###補足情報(言語/FW/ツール等のバージョンなど)
Python 3.5.3 :: Anaconda 2.4.1 (x86_64)を使っています。
python初心者です。よろしくお願いいたします。

Answer 1

自己解決と言っていいのか微妙ですが、、
回答者様の教育的な回答により、正しいプログラムを少し自分で考えながら組むことができました。

def fanction_gauss(param, x, y):
    ll = np.array([])
    def f(x):
        return param[0]*np.exp(-(x-param[2])**2/(2*param[1]))+param[3]
    for i in range(cnt):
        l = integrate.quad(f, x[i]-0.25, x[i]+0.25)
        ll = np.append(ll, l[0])
    residual = y - ll
    return residual

Comments

[magichan]

引数としてベクトルデータに対応する方法として、
回等のように内部でループを回してもよいのですが、numpyには vectorize() というベクトル化関数を作成する関数が用意されておりますので、それを使用するのが簡単です。
今回の場合は、下記のように新しくベクトル対応関数を作成して、これを leastsq() に渡します。

def vectorized_fanction_gauss(param, x, y):
__f = np.vectorize(fanction_gauss, excluded=['param'])
__return f(param=param,x=x,y=y)

[saka_kei]

返信遅れてしまいすみません。
同様の結果を得ることができました！ありがとうございます！

Answer 2

integrate.quadの使い方が間違っています
第一引数は値ではなく関数ですし、戻り値はタプルです。

「値ではなく関数」の意味:

def f(x):
    ...
    return (数値)

に対して

まちがい

integrate.quad(f(x), a, b)

せいかい

integrate.quad(f, a, b)

追記を受けて:

def function_gauss(param, xs, ys):
    def f(x):
        return x*param[0]
    
    fxs = [integrate.quad(f, x-0.25, x+0.25)[0] for x in xs]
    
    residual = ys - fxs
    return residual

Comments

[saka_kei]

早速の回答ありがとうございます。
関数のつもりだったのですが、変数が被っていて値となってしまっていました。
積分範囲を変えながら、上の関数の積分値とyの残渣を取りたいのですが、どのようにしたら良いでしょうか？

[ozwk]

> 関数のつもりだったのですが、変数が被っていて値となってしまっていました。
変数が被っているから値になったのではありません。

ここで言う「関数」とは
def f(x):...(省略)に対して

f(x) <- 関数値
f <- 関数

[saka_kei]

詳しい回答ありがとうございます。

optimize.leastsqに用いる関数の中で、integrate.quadに用いる関数を定義し、指摘していただいたようにプログラムをし直してみましたが、同様のエラーが出てしまいました。

```python
def fanction_gauss(param, x, y):
def f(x):
return param[0]*np.exp(-(x-param[2])**2/(2*param[1]))+param[3]
residual = y - integrate.quad(f, x-0.25, x+0.25)[0]
return residual
```

関数で書いたものでフィッティングはできないですか？
すみません、よろしくお願いします。

[ozwk]

1. fanction_gaussに何が渡っているか確認するため、内部にprint(x)を入れる
2. np配列が渡っていることがわかるので、np配列同士で引き算できるように調整
3. やったー

[saka_kei]

教育的な回答ありがとうございます。
今回の質問により、使い勝手のわからない関数へのアプローチの方法も学ぶことができました。重ねて御礼申し上げますm(_ _)m