Scipy與ROOT等的擬合（高斯）

Question

我現在已經多次偶然發現python中的擬合，scipy.curve_fit在某種程度上比其他工具（例如， ROOT（https://root.cern.ch/）

例如，擬合高斯時，與SciPy的我大多得到一條直線：

相應的代碼：

def fit_gauss(y, x = None): 
    n = len(y) # the number of data 
    if x is None: 
     x = np.arange(0,n,1) 
    mean = y.mean() 
    sigma = y.std() 

    def gauss(x, a, x0, sigma): 
     return a * np.exp(-(x - x0) ** 2/(2 * sigma ** 2)) 

    popt, pcov = curve_fit(gauss, x, y, p0=[max(y), mean, sigma]) 

    plt.plot(x, y, 'b+:', label='data') 
    plt.plot(x, gauss(x, *popt), 'ro:', label='fit') 
    plt.legend() 
    plt.title('Gauss fit for spot') 
    plt.xlabel('Pixel (px)') 
    plt.ylabel('Intensity (a.u.)') 
    plt.show() 

使用ROOT，我得到了完美的結合，甚至沒有給出起始參數：

再次，相應的代碼：

import ROOT 
import numpy as np 

y = np.array([2., 2., 11., 0., 5., 7., 18., 12., 19., 20., 36., 11., 21., 8., 13., 14., 8., 3., 21., 0., 24., 0., 12., 0., 8., 11., 18., 0., 9., 21., 17., 21., 28., 36., 51., 36., 47., 69., 78., 73., 52., 81., 96., 71., 92., 70., 84.,72., 88., 82., 106., 101., 88., 74., 94., 80., 83., 70., 78., 85., 85., 56., 59., 56., 73., 33., 49., 50., 40., 22., 37., 26., 6., 11., 7., 26., 0., 3., 0., 0., 0., 0., 0., 3., 9., 0., 31., 0., 11., 0., 8., 0., 9., 18.,9., 14., 0., 0., 6., 0.]) 
x = np.arange(0,len(y),1) 
#yerr= np.array([0.1,0.2,0.1,0.2,0.2]) 
graph = ROOT.TGraphErrors() 
for i in range(len(y)): 
    graph.SetPoint(i, x[i], y[i]) 
    #graph.SetPointError(i, yerr[i], yerr[i]) 
func = ROOT.TF1("Name", "gaus") 
graph.Fit(func) 

canvas = ROOT.TCanvas("name", "title", 1024, 768) 
graph.GetXaxis().SetTitle("x") # set x-axis title 
graph.GetYaxis().SetTitle("y") # set y-axis title 
graph.Draw("AP") 

有人可以向我解釋，爲什麼結果差異很大？在scipy中的實現是壞的/依賴於良好的啓動參數？ 有沒有辦法解決它？我需要自動處理大量的擬合，但無法訪問目標計算機上的ROOT，因此它只能使用python。

當考慮從根本上契合的結果，並給他們SciPy的作爲啓動參數，配合正常工作與SciPy的，以及...

Answer 1

沒有實際的數據是不容易複製的結果，但創制噪聲數據，它看起來沒什麼問題：

這是我使用的代碼：

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

# your gauss function 
def gauss(x, a, x0, sigma): 
    return a * np.exp(-(x - x0) ** 2/(2 * sigma ** 2)) 

# create some noisy data 
xdata = np.linspace(0, 4, 50) 
y = gauss(xdata, 2.5, 1.3, 0.5) 
y_noise = 0.4 * np.random.normal(size=xdata.size) 
ydata = y + y_noise 
# plot the noisy data 
plt.plot(xdata, ydata, 'bo', label='data') 

# do the curve fit using your idea for the initial guess 
popt, pcov = curve_fit(gauss, xdata, ydata, p0=[ydata.max(), ydata.mean(), ydata.std()]) 

# plot the fit as well 
plt.plot(xdata, gauss(xdata, *popt), 'r-', label='fit') 

plt.show() 

和你一樣，我也使用p0=[ydata.max(), ydata.mean(), ydata.std()]作爲初始猜測，並且對於不同的噪聲強度，它似乎可以很好地工作。

編輯

我剛剛意識到你實際上提供了數據;那麼結果如下所示：

代碼：

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


def gauss(x, a, x0, sigma): 
    return a * np.exp(-(x - x0) ** 2/(2 * sigma ** 2)) 

ydata = np.array([2., 2., 11., 0., 5., 7., 18., 12., 19., 20., 36., 11., 21., 8., 13., 14., 8., 3., 21., 0., 24., 0., 12., 
0., 8., 11., 18., 0., 9., 21., 17., 21., 28., 36., 51., 36., 47., 69., 78., 73., 52., 81., 96., 71., 92., 70., 84.,72., 
88., 82., 106., 101., 88., 74., 94., 80., 83., 70., 78., 85., 85., 56., 59., 56., 73., 33., 49., 50., 40., 22., 37., 26., 
6., 11., 7., 26., 0., 3., 0., 0., 0., 0., 0., 3., 9., 0., 31., 0., 11., 0., 8., 0., 9., 18.,9., 14., 0., 0., 6., 0.]) 

xdata = np.arange(0, len(ydata), 1) 

plt.plot(xdata, ydata, 'bo', label='data') 

popt, pcov = curve_fit(gauss, xdata, ydata, p0=[ydata.max(), ydata.mean(), ydata.std()]) 
plt.plot(xdata, gauss(xdata, *popt), 'r-', label='fit') 

plt.show() 

Answer 2

您可能沒有真正想用ydata.mean()爲高斯質心爲初始值的初始值，或ydata.std()差異 - 這些可能更好地從xdata中猜出。我不知道這是什麼原因造成了最初的麻煩。

您可能會發現lmfit庫有用。這提供了一種方法，可以將模型函數gauss轉換爲Model類，並使用fit()方法使用從模型函數確定的命名參數。使用它，你的配合可能看起來像：

import numpy as np 
import matplotlib.pyplot as plt 

from lmfit import Model 

def gauss(x, a, x0, sigma): 
    return a * np.exp(-(x - x0) ** 2/(2 * sigma ** 2)) 

ydata = np.array([2., 2., 11., 0., 5., 7., 18., 12., 19., 20., 36., 11., 21., 8., 13., 14., 8., 3., 21., 0., 24., 0., 12., 0., 8., 11., 18., 0., 9., 21., 17., 21., 28., 36., 51., 36., 47., 69., 78., 73., 52., 81., 96., 71., 92., 70., 84.,72., 88., 82., 106., 101., 88., 74., 94., 80., 83., 70., 78., 85., 85., 56., 59., 56., 73., 33., 49., 50., 40., 22., 37., 26., 6., 11., 7., 26., 0., 3., 0., 0., 0., 0., 0., 3., 9., 0., 31., 0., 11., 0., 8., 0., 9., 18.,9., 14., 0., 0., 6., 0.]) 

xdata = np.arange(0, len(ydata), 1) 

# wrap your gauss function into a Model 
gmodel = Model(gauss) 
result = gmodel.fit(ydata, x=xdata, 
        a=ydata.max(), x0=xdata.mean(), sigma=xdata.std()) 

print(result.fit_report()) 

plt.plot(xdata, ydata, 'bo', label='data') 
plt.plot(xdata, result.best_fit, 'r-', label='fit') 
plt.show() 

還有幾個附加功能。例如，您可能希望看到的信心在最適合的，這將是（主機版，即將要發佈的）：

# add estimated band of uncertainty: 
dely = result.eval_uncertainty(sigma=3) 
plt.fill_between(xdata, result.best_fit-dely, result.best_fit+dely, color="#ABABAB") 
plt.show() 

給：