适配高斯函数

28

我有一个直方图(如下所示),我正在尝试找到平均值和标准差以及适合直方图的曲线拟合代码。 我认为SciPy或matplotlib中有一些可以帮助的东西,但我尝试过的每个示例都不起作用。

import matplotlib.pyplot as plt
import numpy as np

with open('gau_b_g_s.csv') as f:
    v = np.loadtxt(f, delimiter= ',', dtype="float", skiprows=1, usecols=None)

fig, ax = plt.subplots()

plt.hist(v, bins=500, color='#7F38EC', histtype='step')

plt.title("Gaussian")
plt.axis([-1, 2, 0, 20000])

plt.show()

7
"doesn't work" 的意思是什么?指它无法运行,还是输出结果不正确? - Jodaka
我无法运行从互联网上获取的代码,无法像它们应该做的那样制作曲线。 - user1496646
这很有可能是因为我刚开始学编程,一般都不知道我在做什么。 - user1496646
1
当您尝试运行它时,是否收到错误消息?还是程序在不产生任何结果的情况下完成了? - Jodaka
我只是不知道如何正确地将它与我的数据配合使用。 - user1496646
4个回答

46

请看这个答案,以将任意曲线拟合到数据上。基本上,您可以使用scipy.optimize.curve_fit来将任何函数拟合到您的数据上。下面的代码展示了如何将高斯函数拟合到一些随机数据上(感谢这篇文章在SciPy-User邮件列表中发布)。

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

# Define some test data which is close to Gaussian
data = numpy.random.normal(size=10000)

hist, bin_edges = numpy.histogram(data, density=True)
bin_centres = (bin_edges[:-1] + bin_edges[1:])/2

# Define model function to be used to fit to the data above:
def gauss(x, *p):
    A, mu, sigma = p
    return A*numpy.exp(-(x-mu)**2/(2.*sigma**2))

# p0 is the initial guess for the fitting coefficients (A, mu and sigma above)
p0 = [1., 0., 1.]

coeff, var_matrix = curve_fit(gauss, bin_centres, hist, p0=p0)

# Get the fitted curve
hist_fit = gauss(bin_centres, *coeff)

plt.plot(bin_centres, hist, label='Test data')
plt.plot(bin_centres, hist_fit, label='Fitted data')

# Finally, lets get the fitting parameters, i.e. the mean and standard deviation:
print 'Fitted mean = ', coeff[1]
print 'Fitted standard deviation = ', coeff[2]

plt.show()

谢谢,这样做可以得到均值和标准差,但是曲线拟合实际上并没有产生曲线,而是产生了直线。 - user1496646
1
你的意思是我的例子只能产生线条吗?还是当你将以上代码应用于你的数据时,会有线条出现?此外,直线和曲线之间有什么区别? - Chris
与钟形曲线不同,它看起来就像一个胡萝卜^ - user1496646
没有更多的信息,我无法真正地帮助你。你的数据看起来像胡萝卜吗?如果是这样,那么很可能是因为你的数据就是这个样子。当提问时,最好包括一个简短、自包含的示例 - Chris
3
我猜@user1496646的意思是,他的情况下没有那么多“<bin_centres>”,所以当您绘制(bin_centres, hist_fit)时,会出现采样不足的高斯分布图形(“胡萝卜”形)。他应该只需对bin_centers进行子采样,使用new_bin_centers = numpy.linspace(bin_centres[0], bin_centres[-1], 200),new_hist_fit = gauss(new_bin_centres, *coeff),然后绘制(new_bin_centres, new_hist_fit)。 - SuperElectric

16

您可以尝试使用sklearn高斯混合模型估计,如下所示:

import numpy as np
import sklearn.mixture

gmm = sklearn.mixture.GMM()

# sample data
a = np.random.randn(1000)

# result
r = gmm.fit(a[:, np.newaxis]) # GMM requires 2D data as of sklearn version 0.16
print("mean : %f, var : %f" % (r.means_[0, 0], r.covars_[0, 0]))

参考文献: http://scikit-learn.org/stable/modules/mixture.html#mixture

请注意,使用这种方法,您不需要使用直方图估计您的样本分布。


似乎在某个时候,sklearn.mixture.GMM 已被替换为 sklearn.mixture.GaussianMixture - jjc385

3

虽然这是一个老问题,但如果您只是想为一系列数据制作密度拟合图,可以尝试使用 matplotlib 的 .plot(kind='kde')。文档在这里

以下是 pandas 的示例:

mydf.x.plot(kind='kde')

哇,TIL matplotlib内置了核密度估计。+1 - Joseph Farah

1

我不确定您的输入是什么,但可能您的y轴刻度太大(20000),请尝试减小这个数字。以下代码适用于我:

import matplotlib.pyplot as plt
import numpy as np

#created my variable
v = np.random.normal(0,1,1000)


fig, ax = plt.subplots()


plt.hist(v, bins=500, normed=1, color='#7F38EC', histtype='step')

#plot
plt.title("Gaussian")
plt.axis([-1, 2, 0, 1]) #changed 20000 to 1

plt.show()

编辑:

如果您想要 y 轴上实际值的计数,可以设置 normed=0。并且可以删除 plt.axis([-1, 2, 0, 1])

import matplotlib.pyplot as plt
import numpy as np

#function
v = np.random.normal(0,1,500000)


fig, ax = plt.subplots()

# changed normed=1 to normed=0
plt.hist(v, bins=500, normed=0, color='#7F38EC', histtype='step')

#plot
plt.title("Gaussian")
#plt.axis([-1, 2, 0, 20000]) 

plt.show()

不,我正在处理超过50万个数据点,所以我希望比例尺足够大,因为我不想要只有5万个箱子。 - user1496646
我认为y轴上的值并不告诉你每个区间中观测值的数量,而是告诉你每个区间中的百分比。只需将整个“plt.axis([-1, 2, 0, 1])”行注释掉并运行它,你应该会得到一个分布图。 - Akavall
它肯定告诉我每个箱子中的数字,因为我可以看到直方图本身,y轴在20,000。 - user1496646
Downvoter,您能否解释一下为什么要点踩? - Akavall

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接