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() 

给：