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

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

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

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

Q&A

解決済

1回答

3596閲覧

Pythonでのガウシアンフィッティングが上手くフィットしない

mappys

総合スコア104

Python

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

0グッド

0クリップ

投稿2022/02/06 05:34

画像のように、Pythonでのガウシアンフィッティングが上手くフィットしません。(青い線がフィッティング関数)
イメージ説明
出力すると
OptimizeWarning: Covariance of the parameters could not be estimated
warnings.warn('Covariance of the parameters could not be estimated',
と出るので、初期値の設定が悪いのかな、と思うのですが、初期値はある程度設定できている気がしていて、もし他に何か原因があってわかる方いましたら教えていただきたいです。
コードは以下です。ネットの記事を参考にしました。

python

1#モジュールのインポート 2import pandas as pd 3import numpy as np 4from scipy.optimize import curve_fit 5import matplotlib.pyplot as plt 6import matplotlib.gridspec as gridspec 7 8#データの読み込み 9input_csv = pd.read_csv('元データ.csv') 10channels = input_csv[input_csv.keys()[0]] 11counts = input_csv[input_csv.keys()[1]] 12 13#フィッティング関数の定義 14def quadratic(x, a, b, c): 15 return a + b * x + c * x**2 16 17#ガウシアン 18def gauss(x, a, sigma, mean): 19 return a / np.sqrt(2.0*np.pi) / sigma * np.exp(-((x-mean)/sigma)**2/2) 20 21#フィッティング関数は二次関数とガウシアンをまとめて定義してもよい。 22#ガウシアンを三度呼び出すので、分けて記述した。 23#文頭のアスタリスクのついた*paramaters は引数を幾らでも渡せる。 24#本来ならば引数を個別に指定すべき。 25#今回はあまりにも冗長になるため、*paramaters と指定した。 26def fitting(x, *paramaters): 27 a0, b0, c0, a1, sigma1, mean1, a2, sigma2, mean2, a3, sigma3, mean3, = paramaters 28 return quadratic(x, a0, b0, c0) + \ 29 gauss(x, a1, sigma1, mean1) + \ 30 gauss(x, a2, sigma2, mean2) + \ 31 gauss(x, a3, sigma3, mean3) 32 33#初期値の設定。 34#ガウシアンの場合、meanを適切に設定しないと発散しやすい。 35#線形近似であっても適切に設定することが望ましい。 36#関数の引数が複数ある場合は、タプルで渡すこと。 37initial_quadratic = 10000, -0.02, 0.001 #a0, b0, c0 38initial_gauss1 = 889000, 236000, 200000 #a1, sigma1, mean1 39initial_gauss2 = 824000, 999000, 200000 #a2, sigma2, mean2 40initial_gauss3 = 741000, 1761000, 200000 #a3, sigma3, mean3 41 42#初期値が多いので、個別に指定して接続した。 43initial_parameter = initial_quadratic + initial_gauss1 + initial_gauss2 + initial_gauss3 44 45#フィッティングの実行 46#paramater_optimalが最適化されたパラメータで、covarianceは共分散。 47paramater_optimal, covariance = curve_fit(fitting, channels, counts, initial_parameter) 48 49#matplotlibの設定 50plt.rc('font', family='Arial', size=14) #フォントの設定 51plt.rc('xtick.major', width=1, size=6) #x軸の主目盛りの設定 52plt.rc('xtick', direction='in', top=True) #x軸目盛りの向き、上側に表示するか 53plt.rc('ytick.major', width=1, size=6) #y軸の主目盛りの設定 54plt.rc('ytick', direction='in', right=True) #y軸目盛りの向き、右側に表示するか 55plt.rc('axes', linewidth=1.5) #枠線の太さ 56plt.rc('lines', linewidth=1.0) #プロットの太さ 57 58 59fig = plt.figure(figsize=(8, 6), dpi=80) 60gs = gridspec.GridSpec(1, 1) 61ax1 = fig.add_subplot(gs[0]) 62 63#データをグラフに表示する。 64ax1.plot(channels, counts, ls='-', c='g', label='row data') 65ax1.set_xlabel('Channel', fontsize=16) 66ax1.set_ylabel('Count', fontsize=16) 67ax1.yaxis.set_data_interval(0,10) #set_ylim を使う方が一般的 68 69#フィッティングした曲線の表示 70#np.corrcoef で相関係数を求める。 71#結果が2×2の行列で返るので、1行目2列目もしくは2行目1列目の値を取得。 72ax1.plot(channels, fitting(channels, *paramater_optimal), 73 lw=1.5, ls='--', c='k', label='fitting curve\n' + \ 74 f'$R^2$ = {np.corrcoef(counts, fitting(channels, *paramater_optimal))[0][1]}') 75 76#二次関数のパラメーターを利用したいので、スライスで取り出す。 77a0, b0, c0 = paramater_optimal[0:3] 78 79ax1.plot(channels, quadratic(channels, a0, b0, c0), c='b', 80 label=r'$y = a_0 + b_0x + c_0x^2$' +'\n'+ f'$a_0$ = {a0}\nb$_0$ = {b0}\nc$_0$ = {c0}') 81 82#各meanをグラフ上に表示する。繰り返しのため for 文を利用。 83for i in range(3,12,3): 84#この関数内の、a0, b0, c0 はグローバル変数である。 85 def quadratic_gauss(x, *paramaters): 86 a1, sigma1, mean1, = paramaters 87 return quadratic(x, a0, b0, c0) + \ 88 gauss(x, a1, sigma1, mean1) 89 90#annotate で Axes 内にテキストや矢印を表示できる。 91 mean = paramater_optimal[i+2] 92 ax1.annotate(int(mean), fontsize = 14, 93 xy=(mean, quadratic_gauss(mean, *paramater_optimal[i:i+3])+3), 94 xytext=(mean, quadratic_gauss(mean, *paramater_optimal[i:i+3])+15), 95 ha="center", va='top', 96 arrowprops=dict(shrink=0.05, facecolor='k')) 97 98#凡例の表示。Figure の判例を出力することもできる。 99ax1.legend(loc='upper right', frameon=False) 100 101gs.tight_layout(fig) #グラフ描画前に行うべし。慣例では、fig.tight_layout() 102plt.show() #慣例では plt.show() と fig を指定しない。 103

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

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

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

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

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

jbpb0

2022/02/07 02:18 編集

> 初期値はある程度設定できている気がしていて 初期値が、かなりデータと合ってないと思います 確認方法を以下書きますので、参考にしてみてください まずは、二次関数の初期値を initial_quadratic = 0, 0, 0 #a0, b0, c0 として実質ガウシアンだけにして、 paramater_optimal, covariance = curve_fit(fitting, channels, counts, initial_parameter) をスキップして、その代わりに paramater_optimal = initial_parameter を追加して、コード全体を実行して、初期値のままでグラフを書かせた場合に、ガウシアンの初期値がデータにそこそこ合ってるかを確認してみてください そうすれば、a1〜3の初期値をもっと大きくしないと、ガウシアンの大きさがデータの極大値に全然合わないことが分かると思うので、グラフでガウシアンがデータにそこそこ合うように、ガウシアンの初期値を変えてみてください あと、初期値のsigma1〜3とmean1〜3の順番が逆です sigma1〜3もガウシアンの大きさに効くので、先に初期値のsigma1〜3とmean1〜3の順番を直してから、a1〜3を変えてみることをお勧めします
jbpb0

2022/02/07 02:00 編集

上記を行なって、ガウシアンの初期値がデータにそこそこ合ったら、二次関数の初期値を元に戻して、初期値のままでグラフを書かせた場合にどうなるのかを確認してみてください そうすれば、c0(二次の係数)が大きすぎることが分かると思うので、グラフで二次関数+ガウシアンがデータにそこそこ合うように、二次関数の初期値をを変えてみてください グラフで二次関数+ガウシアンがデータにそこそこ合うように初期値が設定できたら、 paramater_optimal = initial_parameter を無効にして、 paramater_optimal, covariance = curve_fit(fitting, channels, counts, initial_parameter) を有効にしてコード全体を実行したら、フィッティングができると思います 【追記】 当方で、適当に似た感じのデータを作って確認したら、ガウシアンの初期値だけ適切に直したら、二次関数の初期値は質問のコードの数値と同じままでも、フィッティングができました ただし、ガウシアンの初期値も質問のコードの数値と同じにしたら、ダメでした
mappys

2022/02/07 10:32

jbpb0様 回答誠にありがとうございます!! 教えていただいた手順通りに進めたら上手くフィットできました!本当にありがとうございます!! 自分一人ではa1~a3の値をいじることが頭になかったので絶対解決できませんでした。本当に感謝です。 初期値はこんな感じでうまくいきました。 initial_gauss1 = 110000000000, 52993, 236000 #a1, sigma1, mean1 initial_gauss2 = 110000000000, 59083, 999000 #a2, sigma2, mean2 initial_gauss3 = 110000000000, 67539,1761000 あと、一点疑問なのですが、a1~a3っていうのはガウシアンの高さですよね?グラフを見た感じ800000くらいだと思ってそれくらいを初期値で入れてたのですが、実際もっと大きくしないとフィットしなかったですよね。これは一体どういうことなのでしょうか。。 a1~a3=ガウシアンの高さという認識がそもそもおかしいのでしょうか?
mappys

2022/02/07 10:38

たびたび申し訳ないのですが、もう一点質問させてください。 フィットした後のmenanとsigmaの値も一緒に出力したいのですが、meanは #annotate で Axes 内にテキストや矢印を表示できる。 mean = paramater_optimal[i+2] ax1.annotate(int(mean), fontsize = 14, xy=(mean, quadratic_gauss(mean, *paramater_optimal[i:i+3])+3), xytext=(mean, quadratic_gauss(mean, *paramater_optimal[i:i+3])+15), ha="center", va='top', arrowprops=dict(shrink=0.05, facecolor='k')) print(mean)#追加 の部分の最後にprint(mean)を追加で解決したのですが、sigmaを出力する場合どうしたらいいですか?
mappys

2022/02/07 10:41

>>あと、一点疑問なのですが、a1~a3っていうのはガウシアンの高さですよね?グラフを見た感じ800000くらいだと思ってそれくらいを初期値で入れてたのですが、実際もっと大きくしないとフィットしなかったですよね。これは一体どういうことなのでしょうか。。 a1~a3=ガウシアンの高さという認識がそもそもおかしいのでしょうか? すみませんガウシアンの式をちゃんと見ていませんでした。ガウシアンの式 a / np.sqrt(2.0*np.pi) / sigma * np.exp(-((x-mean)/sigma)**2/2) のa / np.sqrt(2.0*np.pi) / sigmaで大きさ1に規格化してるんですね。 これを a * np.exp(-((x-mean)/sigma)**2/2) としたらa=ガウシアンの高さですよね
mappys

2022/02/07 13:50

何度もすみません! >>sigmaを出力する場合どうしたらいいですか? こちらも解決しました! ありがとうございました!! あとは誤差と残差も表示できるよう頑張ってみます!本当にありがとうございました。
guest

回答1

0

自己解決

jbpb0様の回答で、初期値お見直し、解決しました。

投稿2022/02/07 13:51

mappys

総合スコア104

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

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

あなたの回答

tips

太字

斜体

打ち消し線

見出し

引用テキストの挿入

コードの挿入

リンクの挿入

リストの挿入

番号リストの挿入

表の挿入

水平線の挿入

プレビュー

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

ただいまの回答率
85.50%

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

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

質問する

関連した質問