【SciPy】curve_fitを用いてピークフィッティングする方法[Python]

2024年6月7日2024年10月20日

URLをコピーしました！

SciPy

前回、PythonのSciPyでcurve_fitを用いてカーブフィッティングをする方法を紹介しました。

curve_fitでピークフィッティングする方法

curve_fitでピークフィッティングするには、前回のカーブフィッティングと同様「curve_fit(フィッティングする関数, フィッティングしたいX値のリスト, フィッティングしたいY値のリスト)」とします。

そして実際のフィッティングカーブを取得するには、先ほどのフィッティングする関数のX値のリストはそのままに変数のところに「*popt」を追加します。

from scipy.optimize import curve_fit

popt, pcov = curve_fit(gauss, x_list, y_list)

fig = plt.figure()
plt.clf()

plt.scatter(x_list, y_list)
plt.plot(x_list, gauss(x_list, *popt))

plt.show()

実行結果

グラフを作成する部分とピークフィッティングの部分をあわせて書くとこんな感じになります。

import numpy as np
import matplotlib.pyplot as plt
import random
from scipy.optimize import curve_fit

x_min = 0
x_max = 100

def gauss(x_list, a=1, m=0, s=1):
    y_list = [a * np.exp(-(x - m)**2 / (2*s**2)) for x in x_list]
    return y_list

x_list = np.arange(x_min, x_max, 1)
y_list = gauss(x_list, random.random(), random.randint(x_min, x_max+1), random.uniform(0, 10))

popt, pcov = curve_fit(gauss, x_list, y_list)

fig = plt.figure()
plt.clf()

plt.scatter(x_list, y_list)
plt.plot(x_list, gauss(x_list, *popt))

plt.show()

find_peaksと合わせてフィッティング精度を上げる方法

ただ実は上記の方法ではうまくフィッティングできずピークが検出されないことも多いです。

import numpy as np
import matplotlib.pyplot as plt
import random
from scipy.optimize import curve_fit

x_min = 0
x_max = 100

def gauss(x_list, a=1, m=0, s=1):
    y_list = [a * np.exp(-(x - m)**2 / (2*s**2)) for x in x_list]
    return y_list

x_list = np.arange(x_min, x_max, 1)
y_list = gauss(x_list, random.random(), random.randint(x_min, x_max+1), random.uniform(0, 10))

popt, pcov = curve_fit(gauss, x_list, y_list)

fig = plt.figure()
plt.clf()

plt.scatter(x_list, y_list)
plt.plot(x_list, gauss(x_list, *popt))

plt.show()

実行結果

もしフィッティングするピークの変数の値が大体でも分かっているのなら、それを引数に与えると精度が増します。

その際は「curve_fit(フィッティングする関数, Xの値のリスト, Yの値のリスト, p0=変数の初期値)」とします。

import numpy as np
import matplotlib.pyplot as plt
import random
from scipy.optimize import curve_fit

x_min = 0
x_max = 100

def gauss(x_list, a=1, m=0, s=1):
    y_list = [a * np.exp(-(x - m)**2 / (2*s**2)) for x in x_list]
    return y_list

x_list = np.arange(x_min, x_max, 1)

parameter = [random.random(), random.randint(x_min, x_max+1), random.uniform(0, 10)]
y_list = gauss(x_list, parameter[0], parameter[1], parameter[2])

popt, pcov = curve_fit(gauss, x_list, y_list, p0=parameter)

fig = plt.figure()
plt.clf()

plt.scatter(x_list, y_list)
plt.plot(x_list, gauss(x_list, *popt))

plt.show()

実行結果

全てが分からなくても、ピーク位置のような重要なパラメータの初期値を与えるだけでも精度はかなり向上します。

import numpy as np
import matplotlib.pyplot as plt
import random
from scipy.optimize import curve_fit

x_min = 0
x_max = 100

def gauss(x_list, a=1, m=0, s=1):
    y_list = [a * np.exp(-(x - m)**2 / (2*s**2)) for x in x_list]
    return y_list

x_list = np.arange(x_min, x_max, 1)

parameter = [random.random(), random.randint(x_min, x_max+1), random.uniform(0, 10)]
y_list = gauss(x_list, parameter[0], parameter[1], parameter[2])

popt, pcov = curve_fit(gauss, x_list, y_list, p0=([1, parameter[1], 1]))

fig = plt.figure()
plt.clf()

plt.scatter(x_list, y_list)
plt.plot(x_list, gauss(x_list, *popt))

plt.show()

実行結果

ピーク位置に関しては「SciPy.signal」の「find_peaks()」を使えば簡単に取得できます。

複数のピークをフィッティングする方法

ここまで来たら複数のピークフィッティングも試してみましょう。

まず複数のピークをもつグラフを作成するには、それぞれのピークを作成し、足し合わせることで作成できます。

y_list = np.zeros_like(x_list)
for _ in range(num_peaks):
    parameter = [random.random(), random.randint(x_min, x_max+1), random.uniform(0, 10)]
    y_list_temp = gauss(x_list, parameter[0], parameter[1], parameter[2])
    y_list = y_list + np.array(y_list_temp)

また「find_peaks()」で取得したピーク情報を使って、複数のピークにフィッティングさせるには、for関数を使って一つずつピークフィッティングさせていきます。

fit_results = []
for peak in peaks:
    popt, pcov = curve_fit(gauss, x_list, y_list, p0=(1, x_list[peak], 1))
    fit_results.append(gauss(x_list, *popt))

ということで全体としてはこんな感じになりました。

import numpy as np
import matplotlib.pyplot as plt
import random
from scipy.optimize import curve_fit
from scipy.signal import find_peaks

x_min = 0
x_max = 100

num_peaks =3

def gauss(x_list, a=1, m=0, s=1):
    y_list = [a * np.exp(-(x - m)**2 / (2*s**2)) for x in x_list]
    return y_list

x_list = np.arange(x_min, x_max, 1)

y_list = np.zeros_like(x_list)
for _ in range(num_peaks):
    parameter = [random.random(), random.randint(x_min, x_max+1), random.uniform(0, 10)]
    y_list_temp = gauss(x_list, parameter[0], parameter[1], parameter[2])
    y_list = y_list + np.array(y_list_temp)

peaks, _ = find_peaks(y_list)

fit_results = []
for peak in peaks:
    popt, pcov = curve_fit(gauss, x_list, y_list, p0=(1, x_list[peak], 1))
    fit_results.append(gauss(x_list, *popt))

fig = plt.figure()
plt.clf()

plt.scatter(x_list, y_list)
for fit in fit_results:
    plt.plot(x_list, fit)

plt.show()

実行結果