相補誤差関数でフィットしたい（scipy.special.erfc）
1import pandas as pd
2import numpy as np
3from scipy.optimize import curve_fit
4from scipy.special import erfc
5
6def erfc_function(x, amp, mu, sigma, offset):
7 """相補誤差補正関数
8 Args:
9 x(array-like): 測定データ
10 amp(float): 振幅
11 mu(float): 平均
12 sigma(float): 標準偏差
13 offset(float): 切片
14 Returns:
15 f(array-like): 誤差補正関数
16 """
17 f = amp * erfc((x - mu) / sigma) + offset
18 return f
19
20def fit(data: pd.DataFrame, x: str, y: str):
21 """
22 Args:
23 data(pd.DataFrame): 測定データ
24 x(str): X軸のカラム名
25 y(str): Y軸のカラム名
26 Returns:
27 popt(array_like): フィット値
28 perr(array_like): フィット値の誤差
29 x_fit(array_like): フィットした結果のX軸の値
30 y_fit(array_like): フィットした結果のY軸の値
31 """
32
33 copied = data.copy()
34 func = erfc_function
35 x_data = copied[x]
36 y_data = copied[y]
37 p_init = [初期パラメータ]
38
39 popt, pcov = curve_fit(func, x_data, y_data, p_init)
40 perr = np.sqrt(np.diag(pcov))
41
42 xmin = x_data.min()
43 xmax = x_data.max()
44
45 amp, mu, sigma, offset = popt
46 amp_e, mu_e, sigma_e, offset_e = perr
47
48 print(f"Amplitute: {amp} +/- {amp_e}")
49 print(f"Mean: {mu} +/- {mu_e}")
50 print(f"Sigma: {sigma} +/- {sigma_e}")
51 print(f"Offset: {offset} +/- {offset_e}")
52
53 x_fit = np.arange(xmin, xmax, 0.1)
54 y_fit = func(x_fit, amp, mu, sigma, offset)
55
56 return popt, perr, x_fit, y_fit
宇宙線測定時のスレッショルドを決定したときのサンプルです。 スレッショルド値をスライドさせ、得られた検出率の曲線に対して、 相補誤差関数でフィッティングしました。
誤差関数
\[\text{erf}(x) = \frac{2}{\sqrt(\pi)} \int_{t=0}^{x} \exp(-t^{2}) \text{d}t
\]
誤差関数は、ガウス分布（正規分布）の累積分布関数です。 宇宙線測定時の信号に乗るノイズはガウス分布にしたがうと仮定できます。 あるスレッショルド値 \(T\) に対して、信号の大きさ \(x\) が \(x < T\) （\(T\)より小さい）となる確率を算出できます。
相補誤差関数
\[\text{erfc}(x) = 1 - \text{erf}(x)
\]
あるスレッショルド値 \(T\) に対して、信号の大きさ \(x\) が \(T < x\) （\(T\)より大きい）となる確率を算出できます。