模拟高斯噪声数据的Python高斯拟合

Question

我需要使用高斯拟合来插值来自仪器的数据。为此，我考虑使用来自scipy的scipy函数。由于我想在虚拟数据上测试这一功能，然后再在仪器上试用，所以我编写了以下代码来生成含噪声的高斯数据并对其进行拟合：

from scipy.optimize import curve_fit
import numpy
import pylab

# Create a gaussian function

def gaussian(x, a, b, c):
    val = a * numpy.exp(-(x - b)**2 / (2*c**2))
    return val

# Generate fake data.
zMinEntry = 80.0*1E-06
zMaxEntry = 180.0*1E-06
zStepEntry = 0.2*1E-06

x = numpy.arange(zMinEntry, 
                 zMaxEntry, 
                 zStepEntry, 
                 dtype = numpy.float64)    

n = len(x)                 

meanY = zMinEntry + (zMaxEntry - zMinEntry)/2
sigmaY = 10.0E-06

a = 1.0/(sigmaY*numpy.sqrt(2*numpy.pi))
y = gaussian(x, a, meanY, sigmaY) + a*0.1*numpy.random.normal(0, 1, size=len(x))


# Fit
popt, pcov = curve_fit(gaussian, x, y)

# Print results
print("Scale =  %.3f +/- %.3f" % (popt[0], numpy.sqrt(pcov[0, 0])))
print("Offset = %.3f +/- %.3f" % (popt[1], numpy.sqrt(pcov[1, 1])))
print("Sigma =  %.3f +/- %.3f" % (popt[2], numpy.sqrt(pcov[2, 2])))


pylab.plot(x, y, 'ro')
pylab.plot(x, gaussian(x, popt[0], popt[1], popt[2]))
pylab.grid(True)
pylab.show()

不幸的是，这不能正常工作，代码的输出如下：

Scale =  6174.816 +/- 7114424813.672
Offset = 429.319 +/- 3919751917.830
Sigma =  1602.869 +/- 17923909301.176

绘制的结果是(蓝色是拟合函数，红色点是噪声输入数据)：

​

​

我也试着看this的答案，但不知道我的问题出在哪里。我是不是漏掉了什么？还是我以错误的方式使用了curve_fit函数？提前感谢！

Answer 1

我同意奥拉夫的观点，因为这是一个规模问题。最优参数相差许多数量级。然而，缩放生成玩具数据的参数似乎并不能解决实际应用程序的问题。curve_fit使用lestsq，它在数值上近似于Jacobian，其中数值问题是由于尺度的不同而产生的(尝试在curve_fit中使用full_output关键字)。

根据我的经验，通常最好使用fmin，它不依赖近似导数，而只使用函数值。您现在必须编写您自己的最小二乘函数，这是要优化的。

起始值仍然很重要。在您的情况下，您可以通过获取a的最大振幅以及对应的b和c的x值来进行足够好的猜测。

在代码中，如下所示：

from scipy.optimize import curve_fit,fmin
import numpy
import pylab

# Create a gaussian function

def gaussian(x, a, b, c):
    val = a * numpy.exp(-(x - b)**2 / (2*c**2))
    return val

# Generate fake data.
zMinEntry = 80.0*1E-06
zMaxEntry = 180.0*1E-06
zStepEntry = 0.2*1E-06

x = numpy.arange(zMinEntry, 
                 zMaxEntry, 
                 zStepEntry, 
                 dtype = numpy.float64)    

n = len(x)                 

meanY = zMinEntry + (zMaxEntry - zMinEntry)/2
sigmaY = 10.0E-06

a = 1.0/(sigmaY*numpy.sqrt(2*numpy.pi))
y = gaussian(x, a, meanY, sigmaY) + a*0.1*numpy.random.normal(0, 1, size=len(x))

print a, meanY, sigmaY
# estimate starting values from the data
a = y.max()
b = x[numpy.argmax(a)]
c = b

# define a least squares function to optimize
def minfunc(params):
    return sum((y-gaussian(x,params[0],params[1],params[2]))**2)

# fit
popt = fmin(minfunc,[a,b,c])

# Print results
print("Scale =  %.3f" % (popt[0]))
print("Offset = %.3f" % (popt[1]))
print("Sigma =  %.3f" % (popt[2]))


pylab.plot(x, y, 'ro')
pylab.plot(x, gaussian(x, popt[0], popt[1], popt[2]),lw = 2)
pylab.xlim(x.min(),x.max())
pylab.grid(True)
pylab.show()

​

​

Answer 2

看起来一些数值不稳定正在悄悄地进入优化器。试着缩放数据。有以下数据：

zMinEntry = 80.0*1E-06 * 1000
zMaxEntry = 180.0*1E-06 * 1000
zStepEntry = 0.2*1E-06 * 1000
sigmaY = 10.0E-06 * 1000

我得到的估计

Scale =  39.697 +/- 0.526
Offset = 0.130 +/- 0.000
Sigma =  -0.010 +/- 0.000

将其与真值进行比较：

Scale = 39.894228
Offset = 0.13
Sigma = 0.01

西格玛的负号当然可以忽略。

这给出了下面的情节

​

​

Answer 3

正如我在一条评论中所说的，如果您提供了合理的初始猜测，则fit成功，也就是这样调用curve_fit：

popt, pcov = curve_fit(gaussian, x, y, [50000,0.00012,0.00002])

​

​