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

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

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

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

解決済

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

R_naniwa
R_naniwa

総合スコア2

Python 3.x

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

2回答

0リアクション

0クリップ

327閲覧

投稿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/ツールのバージョンなど)

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

以下のような質問にはリアクションをつけましょう

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

リアクションが多い質問は、TOPページの「注目」タブのフィードに表示されやすくなります。

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

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

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

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

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

適切な質問に修正を依頼しましょう。

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

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

まだ回答がついていません

会員登録して回答してみよう

アカウントをお持ちの方は

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

ただいまの回答率
86.12%

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

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

質問する

関連した質問

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

Python 3.x

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