対象読者
- Pythonの基本文法を知っている人
- ヒストグラムを複数のガウス関数の和でフィッティングしたい人
序論
とある実験の結果として２つ以上のピークがありそうなヒストグラムを得たため、ガウス関数の和でフィッティングしようと思いました。方法をググってみたのですが、検索結果のスクリプトは複雑そうなものでした。そこで、よりシンプルなものを作成したので共有します。基本文法以上の部分には解説を加えました。
スクリプト
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.stats import norm
from scipy.optimize import curve_fit
# ガウス関数の定義
def gauss(x, *p):
v = np.zeros_like(x)
for igauss in range(0, int(len(p) / 3)):
v += p[3*igauss] * np.exp(-(x - p[3*igauss+1]) ** 2 / (2 * p[3*igauss+2] ** 2))
return v
# データの作成
hist, bin_edges = np.histogram(data, bins=20)
bin_center = [(e1 + e2) * 0.5 for e1, e2 in zip(bin_edges[:-1], bin_edges[1:])]
# フィッティング
popt, _ = curve_fit(gauss, bin_center, hist, p0=[10.0, 0.0, 0.1, 5, 0.5, 0.1])
# グラフの描画
## 棒グラフ
plt.bar(bin_center, hist, width=0.8 * (bin_center[1] - bin_center[0]))
## フィッティングカーブのデータ取得
x = np.arange(0, 2, 0.01)
y1 = gauss(x, *popt[:])
y2 = gauss(x, *popt[:3])
y3 = gauss(x, *popt[3:])
## フィッティングカーブの描画
plt.plot(x, y1, c='red')
plt.plot(x, y2, c='red')
plt.plot(x, y3, c='red')
## グラフの出力
plt.show()
出力例
解説
ガウス関数の定義
def gauss(x, *p):
v = np.zeros_like(x)
for igauss in range(0, int(len(p) / 3)):
v += p[3*igauss] * np.exp(-(x - p[3*igauss+1]) ** 2 / (2 * p[3*igauss+2] ** 2))
return v
ポイントはこの関数の第２引数です。引数名の前にアスタリスクをつけることで、任意の数の位置引数を受け取れるようにしています。３個のパラメータを受け取ったときには1つのガウス関数で、６個のパラメータを受け取ったときには2つのガウス関数の和を、というように受け取るパラメータの個数に応じてガウス関数の個数を自動変更する仕様にしました。シンプルさ重視で、３の倍数以外の個数のパラメータが指定された場合の対応は実装していません。
データの作成
hist, bin_edges = np.histogram(data, bins=20)
bin_center = [(e1 + e2) * 0.5 for e1, e2 in zip(bin_edges[:-1], bin_edges[1:])]
numpyのhistogram関数の出力は、ヒストグラムと、ビン端の値です。２行目で、ビン端の値からビン中央の値を計算しています。２行目は中級者向けのイディオムで、リストの要素をひとつずらして足し合わせることで、ある要素と次の要素の平均を計算しています。
フィッティング
popt, _ = curve_fit(gauss, bin_center, hist, p0=[10.0, 0.0, 0.1, 5, 0.5, 0.1])
gauss関数は任意の数のパラメータを受け取ります。ここではp0に6つのパラメータを渡しているので、2つのガウス関数の和でフィッティングが行われます。3n個のパラメータを渡すと、n個のガウス関数の和でフィッティングが行われます。
フィッティングカーブのデータ取得
x = np.arange(0, 2, 0.01)
y1 = gauss(x, *popt[:])
y2 = gauss(x, *popt[:3])
y3 = gauss(x, *popt[3:])
２行目以降のgauss関数の呼び出しで、第２引数にアスタリスクがついていることに注目してください。こうすることで、リストのスライスを展開して複数の引数としてgauss関数に渡します。gauss関数の第２引数はアスタリスクつきで定義されているので、複数の引数をタプルとして受け取ります。