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

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

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

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

Q&A

解決済

1回答

1286閲覧

重み付きのガウシアンフィッティング

R_naniwa

総合スコア4

Python

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

0グッド

0クリップ

投稿2022/09/14 07:38

前提

先日teratailでヒストグラムのガウシアンフィッティングに関する質問をさせていただきました。scipy.optimize.curve_fitの機能であるwight(重み)を用いたフィッティングを用いた結果として困っていた複数ピークの同時フィッティングができるようになりました。
しかし、複数を同時にフィッティングすることで、高さが合わないという問題が発生しました。以下の図で示します。イメージ説明

原因として、重みをかけたことにより引っ張られてしまっているのではないかと考えています。

実現したいこと

高さが変わることでσなどのパラメータ、この関数の積分結果に影響を及ぼすため、解決したい。
イメージ説明
別のデータですが、これが理想です。

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

特になし

該当のソースコード

python

1from pylab import * 2from scipy.optimize import curve_fit 3import pandas as pd 4import matplotlib.pyplot as plt 5import numpy as np 6 7#csvファイルの読み込み 8data = pd.read_csv("C:/Users/rine/Desktop/CR-39_analysis/TD531/back/Minor.csv") 9#histの条件 alphaはグラフの透明性 10y,x,_ = plt.hist(data,bins=[i*0.48 for i in range(100)],alpha=.3,label='data',ec='black') 11plt.xlim(0,48) 12#plt.ylim(0,400) 13 14x = (x[1:]+x[:-1])/2 # for len(x)==len(y) 15 16def gauss(x,mu,sigma,A): 17 return A*np.exp(-(x-mu)**2/(2*sigma**2)) 18 19def bimodal(x,mu1,sigma1,A1,mu2,sigma2,A2,mu3,sigma3,A3,mu4,sigma4,A4): 20 return gauss(x,mu1,sigma1,A1)+gauss(x,mu2,sigma2,A2)+gauss(x,mu3,sigma3,A3)+gauss(x,mu4,sigma4,A4) 21 22expected = (13.7,2.87,87, 22.4,1.43,102, 28.5,1.41,367, 33.3,0.48,12000) 23params,cov = curve_fit(bimodal,x,y,expected) 24sigma = sqrt(diag(cov)) 25plt.plot(x, bimodal(x,*params), color='red', lw=3, label='fitting') 26plt.legend() 27plt.show() 28print(params,'\n',sigma) 29 30 31w = np.ones(len(x)) 32w[52:60] = 0.25 # 重みの設定(数値を小さく) 33params, cov = curve_fit(bimodal, x, y, expected, sigma = w) 34 35#print(x[52:60])#重みをかける位置の確認 36y,x,_ = plt.hist(data,bins=[i*0.48 for i in range(100)],alpha=.3,label='data',ec='black') 37plt.xlim(0,48) 38#plt.ylim(0,400) 39#plt.ylim(8000,12000) 40 41plt.plot(x, bimodal(x,*params), color='black', lw = 3, label = 'fitting') 42plt.show() 43print(params, '\n', sigma)

試したこと

重みをつける部分の確認及び変更(複数パターンで確認済み)
curve_fit以外のモジュール(GMM)などを試しましたが、いまだうまくいっていません。

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

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

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

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

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

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

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

jbpb0

2022/09/14 07:51

一番右の高い山と、その左のx=25〜32くらいの範囲の山の、それぞれに違う重みを設定しても、ダメでしょうか?
R_naniwa

2022/09/14 07:56

今回もコメントしていただきありがとうございます。 複数の設定は考慮していませんでした。ありがとうございます。
R_naniwa

2022/09/14 08:28

申し訳ありません。重みの複数セットの仕方を調べたのですが、いまいち見当たりませんでした。 w = np.ones(len(x)) z = np.ones(len(x)) w[52:60] = 0.2 # 重みの設定(数値を小さく) z[67:72] = 0.1 params, cov = curve_fit(bimodal, x, y, expected, sigma = w, sigma = z) このような書き方では keyword argument repeated: sigmaというエラーが出てしまうためできませんでした。
jbpb0

2022/09/14 09:19

w = np.ones(len(x)) w[((x>=25.0) & (x<32.0))] = 0.2 # 右から二番目の山 w[((x>=32.0) & (x<=34.0))] = 0.1 # 一番右の山 みたいな
jbpb0

2022/09/14 09:21

あるいは、別の手ですが、「def bimodal(...」を、四つの「gauss()」の和ではなく、一つの「gauss()」だけにして、xの値によって「gauss(x,mu,sigma,A)」のmu, sigma, Aを切り替えるようにするとか たとえば、xが32以上の場合は、 mu=mu4 sigma=sigma4 A=A4 として、「gauss(x,mu,sigma,A)」を計算 x=25〜32の場合は、 mu=mu3 sigma=sigma3 A=A3 として、「gauss(x,mu,sigma,A)」を計算 みたいに
R_naniwa

2022/09/15 02:01

ありがとうございます。 コードのほかに、出力を見直してみたところ [1.35777807e+01 2.71907485e+00 8.88311099e+01 2.17379257e+01 1.43040496e+00 7.42942013e+01 3.07321557e+01 3.62666499e+00 2.81729152e+02 3.33337823e+01 3.92227071e-01 1.08561073e+04] [7.30938486e-01 7.81206557e-01 2.03320914e+01 6.81463962e-01 7.14845753e-01 2.78762079e+01 3.15555329e-01 3.06862505e-01 1.96258149e+01 2.24755284e-03 2.21541510e-03 5.33108722e+01] [32.4 32.88 33.36 33.84 34.32] [1.35311340e+01 2.65789526e+00 8.94413237e+01 2.23831855e+01 1.86849556e+00 9.24460480e+01 2.82521290e+01 1.16988907e+00 3.42557601e+02 3.33320450e+01 4.02449420e-01 1.10254714e+04] [7.30938486e-01 7.81206557e-01 2.03320914e+01 6.81463962e-01 7.14845753e-01 2.78762079e+01 3.15555329e-01 3.06862505e-01 1.96258149e+01 2.24755284e-03 2.21541510e-03 5.33108722e+01] 見づらくて申し訳ないのですが、このような結果が得られていることがわかりました。 一つ目のカッコは、weightsを用いないフィッテイング結果のパラメータ2つ目はその誤差(?)、3つ目はweightsをかけているxの座標、4,5はweightsを用いたフィッテイング結果のパラメータと誤差(?)を示しています。 weightsを用いないフィッテイングの方はピークトップまできれいにできているのですがA=10856をしめしており、逆に二つ目のグラフ(前提にのっているグラフ)の方はA=11025となっています。 描画がおかしいだけなのでしょうか。 ちなみに、一番右のピークのAの推定値は11110程度です。
jbpb0

2022/09/15 03:20 編集

> weightsを用いないフィッテイングの方はピークトップまできれいにできているのですがA=10856をしめしており、 params = [1.35777807e+01, 2.71907485e+00, 8.88311099e+01, 2.17379257e+01, 1.43040496e+00, 7.42942013e+01, 3.07321557e+01, 3.62666499e+00, 2.81729152e+02, 3.33337823e+01, 3.92227071e-01, 1.08561073e+04] xx = np.arange(0, 50, 0.01) yy = bimodal(xx,*params) print(yy.max()) の結果は「11073.57979265338」 > 二つ目のグラフ(前提にのっているグラフ)の方はA=11025となっています。 params2 = [1.35311340e+01, 2.65789526e+00, 8.94413237e+01, 2.23831855e+01, 1.86849556e+00, 9.24460480e+01, 2.82521290e+01, 1.16988907e+00, 3.42557601e+02, 3.33320450e+01, 4.02449420e-01, 1.10254714e+04] xx = np.arange(0, 50, 0.01) yy2 = bimodal(xx,*params2) print(yy2.max()) の結果は「11025.35684360307」 > 描画がおかしいだけなのでしょうか。 xを0.01間隔で計算した場合の最大値は、大差無いです 上記のxx, yyやxx, yy2を、ヒストグラムに重ね書きしてみてください なお、xの間隔をもっと細かくしても、最大値はほとんど変わりませんでした 【追記】 plt.plot(xx, yy, color='red') plt.plot(xx, yy2, color='black') plt.show() を実行したグラフを見たら分かりますが、「重み」無し(グラフ赤線)は一番高い山のところに別の低い山が重なってるので、A=10856とyyの最大値11074が乖離してます 「重み」有り(グラフ黒線)は山が分離してるので、A=11025とyy2の最大値11025が一致してます
guest

回答1

0

ベストアンサー

グラフ表示での関数プロットの横軸の間隔が荒すぎるのかもしれませんので、間隔をもっと細かくしてみてください

 

python

1plt.plot(x, bimodal(x,*params), color='red', lw=3, label='fitting')

↓ 変更

python

1xx = np.arange(x[0], x[-1]+0.01, 0.01) 2plt.plot(xx, bimodal(xx,*params), color='red', lw=3, label='fitting')

 

python

1plt.plot(x, bimodal(x,*params), color='black', lw = 3, label = 'fitting')

↓ 変更

python

1xx = np.arange(x[0], x[-1]+0.01, 0.01) 2plt.plot(xx, bimodal(xx,*params), color='black', lw = 3, label = 'fitting')

 
「0.01」は、適宜調整してください

投稿2022/09/15 22:29

編集2022/09/21 03:28
jbpb0

総合スコア7651

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

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

R_naniwa

2022/09/21 03:31

コメントありがとうございます。 解決後のフィードバックをかけていなかったこと申し訳ありません。 jbpb0さんの方法で視覚的に綺麗なフィッティングを行うことができました。
guest

あなたの回答

tips

太字

斜体

打ち消し線

見出し

引用テキストの挿入

コードの挿入

リンクの挿入

リストの挿入

番号リストの挿入

表の挿入

水平線の挿入

プレビュー

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

ただいまの回答率
85.49%

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

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

質問する

関連した質問