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

Question

前提

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

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

実現したいこと

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

別のデータですが、これが理想です。

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

特になし

該当のソースコード

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

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

Comments

[jbpb0]

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

[R_naniwa]

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

[R_naniwa]

申し訳ありません。重みの複数セットの仕方を調べたのですが、いまいち見当たりませんでした。
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]

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

みたいな

[jbpb0]

あるいは、別の手ですが、「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]

ありがとうございます。
コードのほかに、出力を見直してみたところ
[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を用いないフィッテイング結果のパラメータ２つ目はその誤差(?)、３つ目はweightsをかけているｘの座標、４，５はweightsを用いたフィッテイング結果のパラメータと誤差(?)を示しています。
weightsを用いないフィッテイングの方はピークトップまできれいにできているのですがA=10856をしめしており、逆に二つ目のグラフ（前提にのっているグラフ）の方はA=11025となっています。
描画がおかしいだけなのでしょうか。
ちなみに、一番右のピークのAの推定値は11110程度です。

[jbpb0]

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

Answer 1

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

　

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

↓ 変更

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

　

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

↓ 変更

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

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

Comments

[R_naniwa]

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