質問をすることでしか得られない、回答やアドバイスがある。

15分調べてもわからないことは、質問しよう!

ただいまの
回答率

87.35%

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

解決済

回答 2

投稿 編集

  • 評価
  • クリップ 1
  • VIEW 4,853

score 11

前提・実現したいこと

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初心者です。よろしくお願いいたします。

  • 気になる質問をクリップする

    クリップした質問は、後からいつでもマイページで確認できます。

    またクリップした質問に回答があった際、通知やメールを受け取ることができます。

    クリップを取り消します

  • 良い質問の評価を上げる

    以下のような質問は評価を上げましょう

    • 質問内容が明確
    • 自分も答えを知りたい
    • 質問者以外のユーザにも役立つ

    評価が高い質問は、TOPページの「注目」タブのフィードに表示されやすくなります。

    質問の評価を上げたことを取り消します

  • 評価を下げられる数の上限に達しました

    評価を下げることができません

    • 1日5回まで評価を下げられます
    • 1日に1ユーザに対して2回まで評価を下げられます

    質問の評価を下げる

    teratailでは下記のような質問を「具体的に困っていることがない質問」、「サイトポリシーに違反する質問」と定義し、推奨していません。

    • プログラミングに関係のない質問
    • やってほしいことだけを記載した丸投げの質問
    • 問題・課題が含まれていない質問
    • 意図的に内容が抹消された質問
    • 過去に投稿した質問と同じ内容の質問
    • 広告と受け取られるような投稿

    評価が下がると、TOPページの「アクティブ」「注目」タブのフィードに表示されにくくなります。

    質問の評価を下げたことを取り消します

    この機能は開放されていません

    評価を下げる条件を満たしてません

    評価を下げる理由を選択してください

    詳細な説明はこちら

    上記に当てはまらず、質問内容が明確になっていない質問には「情報の追加・修正依頼」機能からコメントをしてください。

    質問の評価を下げる機能の利用条件

    この機能を利用するためには、以下の事項を行う必要があります。

回答 2

checkベストアンサー

0

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

投稿

編集

  • 回答の評価を上げる

    以下のような回答は評価を上げましょう

    • 正しい回答
    • わかりやすい回答
    • ためになる回答

    評価が高い回答ほどページの上位に表示されます。

  • 回答の評価を下げる

    下記のような回答は推奨されていません。

    • 間違っている回答
    • 質問の回答になっていない投稿
    • スパムや攻撃的な表現を用いた投稿

    評価を下げる際はその理由を明確に伝え、適切な回答に修正してもらいましょう。

  • 2017/08/03 16:35

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

    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
    ```

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

    キャンセル

  • 2017/08/03 22:24 編集

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

    キャンセル

  • 2017/08/04 02:49

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

    キャンセル

0

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

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

投稿

  • 回答の評価を上げる

    以下のような回答は評価を上げましょう

    • 正しい回答
    • わかりやすい回答
    • ためになる回答

    評価が高い回答ほどページの上位に表示されます。

  • 回答の評価を下げる

    下記のような回答は推奨されていません。

    • 間違っている回答
    • 質問の回答になっていない投稿
    • スパムや攻撃的な表現を用いた投稿

    評価を下げる際はその理由を明確に伝え、適切な回答に修正してもらいましょう。

  • 2017/08/04 08:13 編集

    引数としてベクトルデータに対応する方法として、
    回等のように内部でループを回してもよいのですが、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)

    キャンセル

  • 2017/08/12 23:23

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

    キャンセル

15分調べてもわからないことは、teratailで質問しよう!

  • ただいまの回答率 87.35%
  • 質問をまとめることで、思考を整理して素早く解決
  • テンプレート機能で、簡単に質問をまとめられる

関連した質問

同じタグがついた質問を見る