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

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

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

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

解決済

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

R_naniwa
R_naniwa

総合スコア2

Python

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

1回答

0リアクション

0クリップ

249閲覧

投稿2022/09/14 07:38

前提

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

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

実現したいこと

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

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

特になし

該当のソースコード

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("C:/Users/rine/Desktop/CR-39_analysis/TD531/back/Minor.csv") #histの条件 alphaはグラフの透明性 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.7,2.87,87, 22.4,1.43,102, 28.5,1.41,367, 33.3,0.48,12000) params,cov = curve_fit(bimodal,x,y,expected) sigma = sqrt(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[52:60] = 0.25 # 重みの設定(数値を小さく) params, cov = curve_fit(bimodal, x, y, expected, sigma = w) #print(x[52:60])#重みをかける位置の確認 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.ylim(8000,12000) plt.plot(x, bimodal(x,*params), color='black', lw = 3, label = 'fitting') plt.show() print(params, '\n', sigma)

試したこと

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

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

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

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

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

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

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

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

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

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

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

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

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が一致してます

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

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

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

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

ただいまの回答率
86.12%

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

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

質問する

関連した質問

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

Python

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