高斯拟合 Python

声明:本页面是StackOverFlow热门问题的中英对照翻译,遵循CC BY-SA 4.0协议,如果您需要使用它,必须同样遵循CC BY-SA许可,注明原文地址和作者信息,同时你必须将它归于原作者(不是我):StackOverFlow 原文地址: http://stackoverflow.com/questions/19206332/
Warning: these are provided under cc-by-sa 4.0 license. You are free to use/share it, But you must attribute it to the original authors (not me): StackOverFlow

提示:将鼠标放在中文语句上可以显示对应的英文。显示中英文
时间:2020-08-19 13:08:57  来源:igfitidea点击:

Gaussian fit for Python

pythongaussian

提问by Richard Hsia

I'm trying to fit a Gaussian for my data (which is already a rough gaussian). I've already taken the advice of those here and tried curve_fitand leastsqbut I think that I'm missing something more fundamental (in that I have no idea how to use the command). Here's a look at the script I have so far

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

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

What I get from this is a gaussian-ish shape which is my original data, and a straight horizontal line.

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

enter image description here

在此处输入图片说明

Also, I'd like to plot my graph using points, instead of having them connected. Any input is appreciated!

另外,我想使用点绘制我的图形,而不是将它们连接起来。任何输入表示赞赏!

采纳答案by Developer

Here is corrected code:

这是更正的代码:

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

result:
enter image description here

结果:
在此处输入图片说明

回答by user3406246

You get a horizontal straight line because it did not converge.

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

Better convergence is attained if the first parameter of the fitting (p0) is put as max(y), 5 in the example, instead of 1.

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

回答by James

After losing hours trying to find my error, the problem is your formula:

在花了几个小时试图找到我的错误之后,问题是你的公式:

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

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

This previous formula is wrong, the correct formula is the square root of this!;

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

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

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

Hope this helps

希望这可以帮助

回答by sravani vaddi

There is another way of performing the fit, which is by using the 'lmfit' package. It basically uses the cuve_fit but is much better in fitting and offers complex fitting as well. Detailed step by step instructions are given in the below link. http://cars9.uchicago.edu/software/python/lmfit/model.html#model.best_fit

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

回答by gastro

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

should be

应该

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

回答by strpeter

Explanation

解释

You need good starting values such that the curve_fitfunction converges at "good" values. I can not really say why your fit did not converge (even though the definition of your mean is strange - check below) but I will give you a strategy that works for non-normalized Gaussian-functions like your one.

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

Example

例子

The estimated parameters should be close to the final values (use the weighted arithmetic mean- divide by the sum of all values):

估计参数应该接近最终值(使用加权算术平均值- 除以所有值的总和):

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

I personally prefer using numpy.

我个人更喜欢使用 numpy。

Comment on the definition of the mean (including Developer's answer)

评论均值的定义(包括开发人员的回答)

Since the reviewers did not like my edit on #Developer's code, I will explain for what case I would suggest an improved code. The mean of developer does not correspond to one of the normal definitions of the mean.

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

Your definition returns:

您的定义返回:

>>> sum(x * y)
125

Developer's definition returns:

开发者定义返回:

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

The weighted arithmetic mean:

加权算术平均值:

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

Similarly you can compare the definitions of standard deviation (sigma). Compare with the figure of the resulting fit:

同样,您可以比较标准差 ( sigma)的定义。与得到的拟合图进行比较:

Resulting fit

结果拟合

Comment for Python 2.x users

给 Python 2.x 用户的评论

In Python 2.x you should additionally use the new division to not run into weird results or convert the the numbers before the division explicitly:

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

from __future__ import division

or e.g.

或例如

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