[SciPy] 2. ガウス分布によるカーブフィッティング

はじめに
コード
解説
関連記事
参考

はじめに

SciPyのcurve_fitによりガウシアンフィッティングをデータに適用する方法について解説する。

コード

解説

モジュールのインポート

バージョン

データの生成

norm.rvs(loc=100,scale=5,size=1000)
normは正規分布(ガウス分布)のことで、rvs(Random variates)でランダムな確率変数を得ることができる。ここでは、期待値100、標準偏差5の確率変数を1000個取得している。

hist_1, bins = np.histogram(sample_1, 100, range=(50,150))np.histogramでヒストグラムデータを得ることができる。ヒストグラムの度数とx軸のデータが得られるので、hist_1, binsで受け止めている。
sample_1がヒストグラムにしたいデータで、100はbinsの要素数, rangeでヒストグラムを取る範囲を指定している。

bins = bins[:-1]binsの要素数は境界値が得られるので、この場合、101個ある。そのため、普通にplt.plot(bins,hist_1)でプロットすると、
ValueError: operands could not be broadcast together with shapes (101,) (100,)というエラーが帰ってくる。
そのため、binsの最後尾の要素を1個削除している。

フィッティングする関数の定義

$$\displaystyle f(x)=a \exp\left(-\frac{(x-\mu)^2}{2\sigma^2} \right)$$

ガウス分布関数を定義する。

カーブフィッティング

param_iniはフィッティングパラメータの初期値。ガウス分布の関数の引数を3個にしているので、(強度, 期待値, 標準偏差)の順になっている。

popt, pcov = curve_fit(func, bins, hist_1, p0=param_ini)
funcにフィッティングしたい関数, x軸のデータ(bins)とy軸のデータ(hist_1)を入れて、p0にフィッティングパラメータの初期値を入れる。
poptにフィッティングした結果得られたパラメータ、pcovに共分散が返ってくる。

フィッティングしたデータを生成するには、poptのパラメータをガウス分布(func)に適用すれば良い。

近似曲線の図示

ax.bar(bins,hist_1,width=100/100,alpha=0.5,color='m',align='edge')でヒストグラムを表示している。widthは(ヒストグラムの範囲/binsの要素数)とすれば隙間なく表示できる。
align='edge'とすることで棒グラフが左による。つまり、デフォルトでは、x=50のときは棒のセンターがx=50となるが、'edge'とすることで棒の右側がx=50となる。

ax.plot(bins,fitting,'k')でカーブフィッティングしたガウス分布を表示している。

コードをダウンロード(.pyファイル)

コードをダウンロード(.ipynbファイル)

curve_fitの拡張版であるlmfitでガウシアンフィッティングした例は下記記事を参考願います。

[lmfit] 1. lmfitによる非線形最小２乗ガウシアンフィッティング

lmfitは非線形最小二乗法を用いてカーブフィットするためのライブラリであり、Scipy.optimize.curve_fitの拡張版に位置する。ここでは、lmfitでガウシアンフィッティングする方法について説明する。

参考

Bitly

scipy.optimize.curve_fit — SciPy v1.12.0 Manual