[SciPy] 2. ガウス分布によるcurve_fit

python

はじめに

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)に適用すれば良い。

決定係数R2

決定係数R2を求めるには、

  1. 各点で近似曲線とデータの残差を求める
  2. 残差の2乗を求める(rss:残差平方和)
  3. 各データと全体の平均の差の和の2乗をもとめる(tss:全平方和)
  4. 1-(rss/tss)→決定係数

と計算すれば良い。

近似曲線の図示

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’)でカーブフィッティングしたガウス分布を表示している。

ax.annotate(“$R^2$=”+str(r_squared)[0:5], xy=(0.6, 0.6), xycoords=’axes fraction’)
で求めた決定係数(R2)を図中に挿入している。xycoordsがaxes fractionなので、左下を(0,0)右上を(1,1)とした座標の(0.6,0.6)の位置にR2の値を挿入している。

pythonSciPy
この記事をシェアする
sabopy.comをフォローする
サボテンパイソン

コメント

  1. […] [SciPy] 2. ガウス分布によるcurve_fitSciPyのcurve_fitによるヒストグラムのガウシアンフィッティングsabopy.com2019.01.02 […]