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

Question

前提

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

実現したいこと

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

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

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

該当のソースコード

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

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

Comments

[jbpb0]

curve_fit()を実行する前に、expected(bimodal()の初期値)でヒストグラムとbimodal()のグラフを書いてみて、bimodal()がヒストグラムの形にだいたい沿うようにexpectedを調整してから、curve_fit()を実行してみたら、いかがでしょうか？

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

[R_naniwa]

ご回答ありがとうございます。
申し訳ありません。説明が不足しておりました。
推定値のうちAとmuについてはご指摘のように先に描いたヒストグラム内から読み取ったものを使用しています。

標準偏差や高さなどのパラメータを変更しても同じようなグラフが出力されてしまっています。
そのため、推定（初期値）以外の理由を何か考えていたのですが、何が考えられるのでしょうか？
ご教授いただければ幸いです。

[jbpb0]

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]

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

[jbpb0]

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

重視したいデータの「不確かさ」を小さく設定してフィットで重視させることで、「重み」として使うわけです

[R_naniwa]

丁寧な解説ありがとうございます。
早速重みづけを行ってみたのですが、フィッテイング結果に変化がありませんでした。
何か間違っているのでしょうか…？
質問ばかりになってしまい申し訳ありません。

以下にコードを添付いたします。
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]

> w[25:30] = 0.1

だと、25番目〜29番目のデータの重みを増やしてます
それのxの範囲を、下記を実行して確認してください

print(x[25:30])

だいたいx=25〜30になってますでしょうか？

[R_naniwa]

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

[jbpb0]

w = np.ones(len(x))
w[((x>=25.0) & (x<=29.0))] = 0.1

みたいな

25.0, 29.0, 0.1は、フィット結果を確認しながら、適宜調整してみてください

[jbpb0]

https://teratail.com/questions/9cxf2so0cwd5km#reply-6tqo96a60pi3fr
に

> scipy.optimize.curve_fitの機能であるwight(重み)を用いたフィッティングを用いた結果として困っていた複数ピークの同時フィッティングができるようになりました。

> curve_fit以外のモジュール（GMM）などを試しましたが、いまだうまくいっていません。

と書かれてますが、この質問のcan110さんの回答の方法ではうまくいかなかった、ということでしょうか？

[R_naniwa]

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

言い訳にはなりますが、十分に時間のある時に勉強し、使えるようになろうと考えています。

[jbpb0]

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

今のままだと、解決するかどうか分からない方法をまず参考にしてしまいますよね
それで解決しなければ、無駄な時間を使ってしまいます
また、「ベストアンサーではない」という理由で、私の回答を参考にしないかもしれませんし

[R_naniwa]

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

Answer 1

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

↓ 変更

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

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

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

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

Answer 2

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

Comments

[R_naniwa]

ご回答ありがとうございます。
そのモジュールについて勉強してみます。

ピークを二個に割り切るというのは、あるｘ（２，３個目のピークの間）でrowデータを区切ってしまって、フィッテイングをおこなうということでしょうか？試してみます。
ありがとうございます。

[8524ba23]

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

[R_naniwa]

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