iminuit-由CERN所開發的科學擬合函式庫

[LeoTsai999]

歐洲核子物理組織CERN在C++環境下所開發的ROOT函式庫，是目前高能物理界做資料分析的首選。因應C++較高的使用門檻，CERN將ROOT的部分功能重新以Python編寫成iminuit，降低其使用門檻。相較於scipy僅能進行簡單的最小平方法擬合，iminuit能夠做到更複雜的功能，例如多變數擬合、二維error bar，亦可輕鬆看到各項統計學參數，如chi-squarevalue, Hesse-error等。
函式庫之document請見：https://scikit-hep.org/iminuit/about.html

我認為任何一位理學院學生，都可以學習此函式庫的使用。相較於Origin、Matlab等商用軟體，能使用程式語言自行設計數據分析程式，應該是一位科學家必備的能力。且自行設計之數據分析程式，能讓使用者更明瞭數據經過哪些數學上的處理。

以下為一個iminuit擬合的範例：

import numpy as np
import matplotlib.pyplot as plt
from iminuit import Minuit
from iminuit.cost import LeastSquares


# 計算各層沙重
box_weight = np.array([344.91, 342.84, 311.39, 307.90, 296.67, 285.04, 242.16, 361.03])
box_with_sant_weight = np.array([372.09, 373.98, 399.96, 451.94, 422.66, 322.83, 258.81, 389.57])
sand_weight = box_with_sant_weight - box_weight # 克

print(sand_weight)

# 篩網孔徑 單位mm
box_diameter = [2, 0.84, 0.42, 0.25, 0.149, 0.105, 0.074, 0.035]

# 計算各層沙重百分比
sum_of_sand_weight = np.sum(sand_weight)
for i in range(len(sand_weight)):
    sand_weight[i] = sand_weight[i] / sum_of_sand_weight * 100
    print('第', i+1, '層沙重百分比:', sand_weight[i], '%')

# 計算累積沙重百分比
accumulative_sand_weight = np.array([])
for i in range(len(sand_weight)):
    a = np.sum(sand_weight[:i+1])
    accumulative_sand_weight = np.append(accumulative_sand_weight, a)

# box_diameter = np.log(box_diameter)



ju = np.linspace(box_diameter[0], box_diameter[-1], 10000)

def fitting_function(x, a, b, c, d):
    x = np.array(x)
    return d* np.arctan(a * (x - b)) + c


least_squares = LeastSquares(box_diameter, accumulative_sand_weight, 0.001, fitting_function) #(x座標,y座標,y的誤差,想要拿來貼的線的形式[用def來做])
m = Minuit(least_squares, a=1, b=1, c=2, d=-1)  #"="用來設定起始值
m.migrad()  
m.hesse()  
plt.plot(ju, fitting_function(ju, *m.values), label="distribution", color = '#e76254', linewidth = 3.5, zorder=1)    #畫線,*m.values同時代表 阿法 和 碧塔
plt.scatter(box_diameter, accumulative_sand_weight, marker='x',label='data', color='black', zorder=2, s=30)


fit_info = [
    f"$\\chi^2$ / $n_\\mathrm{{dof}}$ = {m.fval:.1f} / {len(box_diameter) - m.nfit}",
]
for p, v, e in zip(m.parameters, m.values, m.errors):
    fit_info.append(f"{p} = ${v:.3f} \\pm {e:.3f}$")

plt.legend(title="\n".join(fit_info), fontsize=16)

plt.xlabel('Particle Diameter (mm)')
plt.ylabel('Percent Passing')

for i in range(len(m.values)):
    print('mvalues', m.values[i])

plt.legend(title="\n".join(fit_info), fontsize=16)
plt.show()

此程式在對數據進行預處理後，定義副程式fitting_function(x, a, b, c, d)為擬合方程式。在這個例子中，
我們使用arctan 函數進行擬合。接下來利用函式LeastSquares定義擬合之數據和誤差。接著使用Minuit設定擬合初始值，並使用migrad()和hesse()進行擬和數學運算。其中migrad是用來計算函數的最小值（在此案例中用來做最小平方法），而hesse()則和誤差有關。

[gigitilto]

我該怎麼定義擬合函數？是要根據理論公式還是根據數據找出最適合的方程式來擬合？