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

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

新規登録して質問してみよう
ただいま回答率
85.44%
Python

Pythonは、コードの読みやすさが特徴的なプログラミング言語の1つです。 強い型付け、動的型付けに対応しており、後方互換性がないバージョン2系とバージョン3系が使用されています。 商用製品の開発にも無料で使用でき、OSだけでなく仮想環境にも対応。Unicodeによる文字列操作をサポートしているため、日本語処理も標準で可能です。

Q&A

解決済

2回答

6725閲覧

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

saka_kei

総合スコア11

Python

Pythonは、コードの読みやすさが特徴的なプログラミング言語の1つです。 強い型付け、動的型付けに対応しており、後方互換性がないバージョン2系とバージョン3系が使用されています。 商用製品の開発にも無料で使用でき、OSだけでなく仮想環境にも対応。Unicodeによる文字列操作をサポートしているため、日本語処理も標準で可能です。

0グッド

1クリップ

投稿2017/08/01 15:35

編集2017/08/03 07:39

###前提・実現したいこと
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()

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

python

1import sys, csv 2import numpy as np 3from scipy import integrate 4sys.path.append("") 5import matplotlib.pyplot as plt 6 7cnt = 0 8x = np.array([]) 9y = np.array([]) 10da_int = np.array([]) 11 12f = open("beamsize_y_up_45mm.csv", "r", encoding="utf_8_sig") 13try: 14 reader = csv.reader(f) 15 for row in reader: 16 x = np.append(x, float(row[0])) 17 y = np.append(y, float(row[1])) 18 cnt = cnt + 1 19finally: 20 f.close() 21 22def 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 26param0 = [1.2, 1, 10.5, -24.3] 27param = optimize.leastsq(fanction_gauss, param0, args=(x, y))

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

python

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

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

python

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

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

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

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

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

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

バッドをするには、ログインかつ

こちらの条件を満たす必要があります。

guest

回答2

0

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

python

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

投稿2017/08/03 17:50

saka_kei

総合スコア11

バッドをするには、ログインかつ

こちらの条件を満たす必要があります。

magichan

2017/08/03 23:19 編集

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

2017/08/12 14:23

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

0

ベストアンサー

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

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

python

1def f(x): 2 ... 3 return (数値)

に対して

まちがい

python

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

せいかい

python

1integrate.quad(f, a, b)

追記を受けて:

python

1def function_gauss(param, xs, ys): 2 def f(x): 3 return x*param[0] 4 5 fxs = [integrate.quad(f, x-0.25, x+0.25)[0] for x in xs] 6 7 residual = ys - fxs 8 return residual

投稿2017/08/01 23:46

編集2017/08/03 13:45
ozwk

総合スコア13544

バッドをするには、ログインかつ

こちらの条件を満たす必要があります。

saka_kei

2017/08/02 02:53

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

2017/08/02 02:59

> 関数のつもりだったのですが、変数が被っていて値となってしまっていました。 変数が被っているから値になったのではありません。 ここで言う「関数」とは def f(x):...(省略)に対して f(x) <- 関数値 f <- 関数
saka_kei

2017/08/03 07: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 ``` 関数で書いたものでフィッティングはできないですか? すみません、よろしくお願いします。
ozwk

2017/08/03 13:47 編集

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

2017/08/03 17:49

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

あなたの回答

tips

太字

斜体

打ち消し線

見出し

引用テキストの挿入

コードの挿入

リンクの挿入

リストの挿入

番号リストの挿入

表の挿入

水平線の挿入

プレビュー

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

ただいまの回答率
85.44%

質問をまとめることで
思考を整理して素早く解決

テンプレート機能で
簡単に質問をまとめる

質問する

関連した質問