Python 模拟高斯噪声数据的高斯拟合
Python gaussian fit on simulated gaussian noisy data
我需要使用高斯拟合对来自仪器的数据进行插值。为此,我考虑使用 scipy
中的 curve_fit
函数。
由于我想在仪器上尝试之前在假数据上测试此功能,因此我编写了以下代码来生成嘈杂的高斯数据并对其进行拟合:
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
函数?提前致谢!
看起来一些数值不稳定性正悄悄进入优化器。尝试缩放数据。具有以下数据:
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
sigma的负号当然可以忽略
这给出了以下情节
就规模问题而言,我同意 Olaf 的观点。最佳参数相差许多数量级。但是,缩放用于生成玩具数据的参数似乎并不能解决实际应用程序的问题。 curve_fit
uses lestsq
,它在数值上近似于雅可比行列式,其中由于尺度的差异而出现数值问题(尝试使用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()
正如我在评论中所说,如果您提供合理的初始猜测,则拟合成功,即调用 curve_fit
那样:
popt, pcov = curve_fit(gaussian, x, y, [50000,0.00012,0.00002])
我需要使用高斯拟合对来自仪器的数据进行插值。为此,我考虑使用 scipy
中的 curve_fit
函数。
由于我想在仪器上尝试之前在假数据上测试此功能,因此我编写了以下代码来生成嘈杂的高斯数据并对其进行拟合:
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
函数?提前致谢!
看起来一些数值不稳定性正悄悄进入优化器。尝试缩放数据。具有以下数据:
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
sigma的负号当然可以忽略
这给出了以下情节
就规模问题而言,我同意 Olaf 的观点。最佳参数相差许多数量级。但是,缩放用于生成玩具数据的参数似乎并不能解决实际应用程序的问题。 curve_fit
uses lestsq
,它在数值上近似于雅可比行列式,其中由于尺度的差异而出现数值问题(尝试使用curve_fit
中的full_output
关键字)。
根据我的经验,通常最好使用 fmin
,它不依赖于近似导数,而仅使用函数值。您现在必须编写自己的要优化的最小二乘函数。
初始值仍然很重要。在您的情况下,您可以通过获取 a
的最大幅度和 b
和 c
.
在代码中,它看起来像这样:
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()
正如我在评论中所说,如果您提供合理的初始猜测,则拟合成功,即调用 curve_fit
那样:
popt, pcov = curve_fit(gaussian, x, y, [50000,0.00012,0.00002])