Python 的高斯拟合

Question

我正在尝试为我的数据拟合一个高斯（这已经是一个粗略的高斯）。我已经听取了这里的建议并尝试了curve_fit 和leastsq，但我认为我缺少一些更基本的东西（因为我不知道如何使用该命令）。 这是我到目前为止的脚本

import pylab as plb
import matplotlib.pyplot as plt

# Read in data -- first 2 rows are header in this example. 
data = plb.loadtxt('part 2.csv', skiprows=2, delimiter=',')

x = data[:,2]
y = data[:,3]
mean = sum(x*y)
sigma = sum(y*(x - mean)**2)

def gauss_function(x, a, x0, sigma):
    return a*np.exp(-(x-x0)**2/(2*sigma**2))
popt, pcov = curve_fit(gauss_function, x, y, p0 = [1, mean, sigma])
plt.plot(x, gauss_function(x, *popt), label='fit')

# plot data

plt.plot(x, y,'b')

# Add some axis labels

plt.legend()
plt.title('Fig. 3 - Fit for Time Constant')
plt.xlabel('Time (s)')
plt.ylabel('Voltage (V)')
plt.show()

我从中得到的是一个高斯形状，这是我的原始数据，以及一条水平直线。

另外，我想用点来绘制我的图表，而不是让它们连接起来。 任何输入表示赞赏！

Answer 1

这里是更正的代码：

import pylab as plb
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
from scipy import asarray as ar,exp

x = ar(range(10))
y = ar([0,1,2,3,4,5,4,3,2,1])

n = len(x)                          #the number of data
mean = sum(x*y)/n                   #note this correction
sigma = sum(y*(x-mean)**2)/n        #note this correction

def gaus(x,a,x0,sigma):
    return a*exp(-(x-x0)**2/(2*sigma**2))

popt,pcov = curve_fit(gaus,x,y,p0=[1,mean,sigma])

plt.plot(x,y,'b+:',label='data')
plt.plot(x,gaus(x,*popt),'ro:',label='fit')
plt.legend()
plt.title('Fig. 3 - Fit for Time Constant')
plt.xlabel('Time (s)')
plt.ylabel('Voltage (V)')
plt.show()

结果：

Answer 2

说明

您需要良好的起始值，以便curve_fit 函数收敛于“良好”值。我真的不能说为什么你的拟合没有收敛（即使你的平均值的定义很奇怪 - 请检查下面），但我会给你一个适用于像你这样的非归一化高斯函数的策略。

示例

估计的参数应该接近最终值（使用weighted arithmetic mean - 除以所有值的总和）：

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

x = np.arange(10)
y = np.array([0, 1, 2, 3, 4, 5, 4, 3, 2, 1])

# weighted arithmetic mean (corrected - check the section below)
mean = sum(x * y) / sum(y)
sigma = np.sqrt(sum(y * (x - mean)**2) / sum(y))

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), 'r-', label='fit')
plt.legend()
plt.title('Fig. 3 - Fit for Time Constant')
plt.xlabel('Time (s)')
plt.ylabel('Voltage (V)')
plt.show()

我个人更喜欢使用 numpy。

评论均值的定义（包括开发者的回答）

由于审阅者不喜欢我的edit on #Developer's code，我将解释在什么情况下我会建议改进代码。开发人员的平均值不对应于平均值的正常定义之一。

你的定义返回：

>>> sum(x * y)
125

开发者定义返回：

>>> sum(x * y) / len(x)
12.5 #for Python 3.x

加权算术平均值：

>>> sum(x * y) / sum(y)
5.0

同样，您可以比较标准差的定义 (sigma)。与结果拟合的数字进行比较：

Python 2.x 用户的评论

在 Python 2.x 中，您应该另外使用新的除法来避免出现奇怪的结果或显式转换除法之前的数字：

from __future__ import division

或例如

sum(x * y) * 1. / sum(y)

Answer 3

你得到一条水平直线，因为它没有收敛。

如果将拟合的第一个参数 (p0) 设置为 max(y)，在示例中为 5，而不是 1，则可以实现更好的收敛。

Answer 4

在花费数小时试图找出我的错误之后，问题在于你的公式：

sigma = sum(y*(x-mean)**2)/n

前面这个公式是错误的，正确的公式是这个的平方根！;

sqrt(sum(y*(x-mean)**2)/n)

希望对你有帮助

Answer 5

还有另一种执行拟合的方法，即使用“lmfit”包。它基本上使用cuve_fit，但在拟合方面要好得多，并且还提供复杂的拟合。 下面的链接中给出了详细的分步说明。 http://cars9.uchicago.edu/software/python/lmfit/model.html#model.best_fit

Answer 6

sigma = sum(y*(x - mean)**2)

应该是

sigma = np.sqrt(sum(y*(x - mean)**2))

Answer 7

实际上，您不需要进行第一次猜测。简单的做

import matplotlib.pyplot as plt  
from scipy.optimize import curve_fit
from scipy import asarray as ar,exp

x = ar(range(10))
y = ar([0,1,2,3,4,5,4,3,2,1])

n = len(x)                          #the number of data
mean = sum(x*y)/n                   #note this correction
sigma = sum(y*(x-mean)**2)/n        #note this correction

def gaus(x,a,x0,sigma):
    return a*exp(-(x-x0)**2/(2*sigma**2))

popt,pcov = curve_fit(gaus,x,y)
#popt,pcov = curve_fit(gaus,x,y,p0=[1,mean,sigma])

plt.plot(x,y,'b+:',label='data')
plt.plot(x,gaus(x,*popt),'ro:',label='fit')
plt.legend()
plt.title('Fig. 3 - Fit for Time Constant')
plt.xlabel('Time (s)')
plt.ylabel('Voltage (V)')
plt.show()

工作正常。这更简单，因为进行猜测并非易事。我有更复杂的数据，没有做出正确的第一次猜测，但简单地删除第一次猜测就可以了：)

P.S.：使用 numpy.exp() 更好，表示 scipy 的警告

Answer 8

这个问题也在这里得到解答：How can I fit a gaussian curve in python?

在 MSeifert 的回答中，还提到了使用 astropy 包的简单高斯建模。