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

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

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

Python 3はPythonプログラミング言語の最新バージョンであり、2008年12月3日にリリースされました。

Q&A

解決済

2回答

1282閲覧

ガウシアンフィッティング

R_naniwa

総合スコア4

Python 3.x

Python 3はPythonプログラミング言語の最新バージョンであり、2008年12月3日にリリースされました。

0グッド

0クリップ

投稿2022/09/08 05:03

前提

私はpythonを使ってデータの解析を行っています。
楕円の短径分布をヒストグラムで作成し、できたピークに対しガウシアンフィッティングを行っています。
今回は4つのピークがあり、それらをすべて同時にフィッティングしたいのですがうまくできません。
一枚目はグラフの全体像。
二枚目は拡大し、問題の箇所を見やすくした図。

イメージ説明
イメージ説明

実現したいこと

・以下のコードを改良または新しいコードを用いてすべてのピークをフィッティングしたい。
・特に重要視しているピークは右の二つなので、4つすべてが不可能な場合それらをフィッティングしたい。

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

エラーは特にありません。

該当のソースコード

python

from pylab import * from scipy.optimize import curve_fit import pandas as pd import matplotlib.pyplot as plt import numpy as np #csvファイルの読み込み data = pd.read_csv("Minor.csv") y,x,_ = plt.hist(data,bins=[i*0.48 for i in range(100)],alpha=.3,label='data',ec='black') plt.xlim(0,48) plt.ylim(0,400) x = (x[1:]+x[:-1])/2 # for len(x)==len(y) #gaussian function def gauss(x,mu,sigma,A): return A*np.exp(-(x-mu)**2/(2*sigma**2)) #summation gaussian function def bimodal(x,mu1,sigma1,A1,mu2,sigma2,A2,mu3,sigma3,A3,mu4,sigma4,A4): return gauss(x,mu1,sigma1,A1)+gauss(x,mu2,sigma2,A2)+gauss(x,mu3,sigma3,A3)+gauss(x,mu4,sigma4,A4) #Guess parameters expected = (13.6,2.50,74.0, 21.0,1.96,80.9, 26.8,1.96,340, 32.76,0.48,9085) params,cov = curve_fit(bimodal,x,y,expected) sigma = np.sqrt(np.diag(cov)) plt.plot(x,bimodal(x,*params),color='red',lw=3,label='fitting') plt.legend() plt.show() print(params,'\n',sigma)

試したこと

フィッティングするピークの数を変更して、その数分(1~3ピーク)でかかるかなどを試した。

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

ここにより詳細な情報を記載してください。

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

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

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

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

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

jbpb0

2022/09/08 07:01 編集

curve_fit()を実行する前に、expected(bimodal()の初期値)でヒストグラムとbimodal()のグラフを書いてみて、bimodal()がヒストグラムの形にだいたい沿うようにexpectedを調整してから、curve_fit()を実行してみたら、いかがでしょうか? 質問の二つ目のグラフ画像を見ると、右から二番目のピークでヒストグラムの形とbimodal()があまり合ってないように思います 右から二番目のピークに相当するガウス分布の標準偏差をもっと小さくする等して、ヒストグラムの形ともっと合うように初期値を設定してからcurve_fit()を実行したら、うまくいきそうな気がします
R_naniwa

2022/09/08 07:25

ご回答ありがとうございます。 申し訳ありません。説明が不足しておりました。 推定値のうちAとmuについてはご指摘のように先に描いたヒストグラム内から読み取ったものを使用しています。 標準偏差や高さなどのパラメータを変更しても同じようなグラフが出力されてしまっています。 そのため、推定(初期値)以外の理由を何か考えていたのですが、何が考えられるのでしょうか? ご教授いただければ幸いです。
jbpb0

2022/09/08 10:54 編集

curve_fit()は引数「sigma」で「重み」的なものを設定できるので、うまくいってないx=25〜30の「重み」が大きくなるように設定してみたらいかがでしょうか? x, yと同じサイズの配列を用意して、それに値を入れます 値の逆数が「重み」になるので、x=25〜30に相当するところは値を小さくします 参考 http://www.yamamo10.jp/yamamoto/comp/Python/library/SciPy/fit/index.php 以下、「重み」を設定する例です 下記のように、対称な山型のデータ(上がって下がる)に直線をフィットさせようとしても、うまくいきません def lf(x, a0, a1): # 直線 return x * a1 + a0 # この行はインデント有り x = np.arange(22) y = np.array([range(11)] + [range(10, -1, -1)]).ravel() # 対称な山型のデータ(上がって下がる) p0 = [0.0, 0.0] params, cov = curve_fit(lf, x, y, p0) plt.plot(x, y) plt.plot(x, lf(x,*params), color='red') plt.show() データ前半の「重み」を大きくすれば、前半のデータの傾きに合うように直線がフィットされます w = np.ones(len(x)) w[0:11] = 0.1 # 前半の重みを大きく(数値を小さく) params, cov = curve_fit(lf, x, y, p0, sigma=w) plt.plot(x, y) plt.plot(x, lf(x,*params), color='red') plt.show()
R_naniwa

2022/09/09 00:57

ありがとうございます。 重みを設定できることを知りませんでした。 試してみます。
jbpb0

2022/09/09 01:53 編集

curve_fit()の引数「sigma」の本来の意味は、「y」データの「不確かさ」(誤差)です ・「不確かさ」が大きいデータは、比較的に信用できないので、フィットで軽視する(合わなくてもあまり気にしない) ・「不確かさ」が小さいデータは、比較的に信用できるので、フィットで重視する 重視したいデータの「不確かさ」を小さく設定してフィットで重視させることで、「重み」として使うわけです
R_naniwa

2022/09/09 09:31

丁寧な解説ありがとうございます。 早速重みづけを行ってみたのですが、フィッテイング結果に変化がありませんでした。 何か間違っているのでしょうか…? 質問ばかりになってしまい申し訳ありません。 以下にコードを添付いたします。 from pylab import * from scipy.optimize import curve_fit import pandas as pd import matplotlib.pyplot as plt import numpy as np #csvファイルの読み込み data = pd.read_csv("Minor.csv") y,x,_ = plt.hist(data,bins=[i*0.48 for i in range(100)],alpha=.3,label='data',ec='black') plt.xlim(0,48) plt.ylim(0,400) x = (x[1:]+x[:-1])/2 # for len(x)==len(y) def gauss(x,mu,sigma,A): return A*np.exp(-(x-mu)**2/(2*sigma**2)) def bimodal(x,mu1,sigma1,A1,mu2,sigma2,A2,mu3,sigma3,A3,mu4,sigma4,A4): return gauss(x,mu1,sigma1,A1)+gauss(x,mu2,sigma2,A2)+gauss(x,mu3,sigma3,A3)+gauss(x,mu4,sigma4,A4) expected = (13.6,2.67,74.0, 21.0,1.39,80.9, 26.8,0.48,340, 32.76,0.48,9085) params,cov = curve_fit(bimodal, x, y, expected) sigma = np.sqrt(np.diag(cov)) plt.plot(x, bimodal(x,*params) ,color = 'red', lw = 3, label = 'fitting') plt.legend() plt.show() print(params,'\n',sigma) w = np.ones(len(x)) w[25:30] = 0.1 # 前半の重みを大きく(数値を小さく) params, cov = curve_fit(bimodal, x, y, expected, sigma = w) y,x,_ = plt.hist(data,bins=[i*0.48 for i in range(100)],alpha=.3,label='data',ec='black') plt.xlim(0,48) plt.ylim(0,400) plt.plot(x, bimodal(x,*params), color='black', lw = 3, label = 'fitting') plt.show()
jbpb0

2022/09/09 11:16 編集

> w[25:30] = 0.1 だと、25番目〜29番目のデータの重みを増やしてます それのxの範囲を、下記を実行して確認してください print(x[25:30]) だいたいx=25〜30になってますでしょうか?
R_naniwa

2022/09/10 14:40

ありがとうございます 完全に解釈を間違えておりました。
jbpb0

2022/09/12 07:48

w = np.ones(len(x)) w[((x>=25.0) & (x<=29.0))] = 0.1 みたいな 25.0, 29.0, 0.1は、フィット結果を確認しながら、適宜調整してみてください
jbpb0

2022/09/21 03:22

https://teratail.com/questions/9cxf2so0cwd5km#reply-6tqo96a60pi3fr に > scipy.optimize.curve_fitの機能であるwight(重み)を用いたフィッティングを用いた結果として困っていた複数ピークの同時フィッティングができるようになりました。 > curve_fit以外のモジュール(GMM)などを試しましたが、いまだうまくいっていません。 と書かれてますが、この質問のcan110さんの回答の方法ではうまくいかなかった、ということでしょうか?
R_naniwa

2022/09/21 03:35

コメントありがとうございます。 質問の答えとしては”うまくいかなかったというよりかは、重みを用いたフィッティングで解決したため、途中でこの方法をやめた”というのが適切です。 現状ほかのデータに対しても、weightでできていますので、GMMによるフィッティングは行っていない感じです。 言い訳にはなりますが、十分に時間のある時に勉強し、使えるようになろうと考えています。
jbpb0

2022/09/26 03:31

https://teratail.com/tour に 「質問・回答によって生まれたコンテンツを、同じ問題を持った人に最適な形で届けます。」 と書かれてるように、teratailには同じことに悩む第三者の手助けになる、という役割もありますので、実際に解決した回答をベストアンサーにした方がいいと思います 同じことに悩んで検索してここに来た人は、まずはベストアンサーを参考にする場合が多いと思うからです 今のままだと、解決するかどうか分からない方法をまず参考にしてしまいますよね それで解決しなければ、無駄な時間を使ってしまいます また、「ベストアンサーではない」という理由で、私の回答を参考にしないかもしれませんし
R_naniwa

2022/09/26 05:46

申し訳ありません。 解決後について配慮がかけておりました。 ベストアンサーの方変更させて頂きました。
guest

回答2

0

ベストアンサー

python

1params,cov = curve_fit(bimodal,x,y,expected)

↓ 変更

python

1w = np.ones(len(x)) 2w[((x>=25.0) & (x<=29.0))] = 0.1 3params,cov = curve_fit(bimodal,x,y,expected,sigma=w)

とすれば、x=25.0〜29.0の範囲がフィットで重視されます

25.0, 29.0, 0.1は、フィット結果を確認しながら、適宜調整してみてください
「0.1」は、数値が小さいほど重視されます

参考
SciPyフィッティング (fitting)

投稿2022/09/12 08:01

jbpb0

総合スコア7651

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

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

0

sklearnにそのものズバリのsklearn.mixture.GaussianMixtureというものがあるので、それを利用してみてはいかがでしょうか。
具体的な使い方はガウス混合分布のパラメータをscikit-learnで推定するなどが参考になりそうです。
なお、右端の値のでかいピークに引っ張られて、分散も大きく推定されすぎているように見えますが、ピーク数を2個に割り切って限定してみると結果は変わりそうな気がします。

投稿2022/09/08 07:11

can110

総合スコア38233

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

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

R_naniwa

2022/09/08 07:28

ご回答ありがとうございます。 そのモジュールについて勉強してみます。 ピークを二個に割り切るというのは、あるx(2,3個目のピークの間)でrowデータを区切ってしまって、フィッテイングをおこなうということでしょうか?試してみます。 ありがとうございます。
can110

2022/09/08 07:31

おそらく与えるデータはそのままで、n_components=2にするだけでよいかと思います。 右側の大きな値だけでうまく推定してくれないものかとの希望的観測ですが。
R_naniwa

2022/09/08 07:55

なるほど!ありがとうございます。 試してみます。
guest

あなたの回答

tips

太字

斜体

打ち消し線

見出し

引用テキストの挿入

コードの挿入

リンクの挿入

リストの挿入

番号リストの挿入

表の挿入

水平線の挿入

プレビュー

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

ただいまの回答率
85.50%

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

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

質問する

関連した質問