gpt4 book ai didi

python - 如何在 Python 中拟合双高斯分布?

转载 作者:太空狗 更新时间:2023-10-30 02:18:58 25 4
gpt4 key购买 nike

我正在尝试使用 Python 获得数据 ( link ) 的双高斯分布。原始数据的形式是:

enter image description here

对于给定的数据,我想获得图中所示峰值的两个高斯分布。我尝试使用以下代码 ( source ):

from sklearn import mixture
import matplotlib.pyplot
import matplotlib.mlab
import numpy as np
from pylab import *
data = np.genfromtxt('gaussian_fit.dat', skiprows = 1)
x = data[:, 0]
y = data[:, 1]
clf = mixture.GMM(n_components=2, covariance_type='full')
clf.fit((y, x))
m1, m2 = clf.means_
w1, w2 = clf.weights_
c1, c2 = clf.covars_
fig = plt.figure(figsize = (5, 5))
plt.subplot(111)
plotgauss1 = lambda x: plot(x,w1*matplotlib.mlab.normpdf(x,m1,np.sqrt(c1))[0], linewidth=3)
plotgauss2 = lambda x: plot(x,w2*matplotlib.mlab.normpdf(x,m2,np.sqrt(c2))[0], linewidth=3)
fig.savefig('gaussian_fit.pdf')

但我无法获得所需的输出。那么,如何在Python中得到双高斯分布呢?

更新

我能够使用以下代码拟合单个高斯分布:

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

data = np.genfromtxt('gaussian_fit.dat', skiprows = 1)
x = data[:, 0]
y = data[:, 1]
n = len(x)
mean = sum(x*y)/n
sigma = sum(y*(x-mean)**2)/n


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


fig = plt.figure(figsize = (5, 5))
plt.subplot(111)
plt.plot(x, y, label='Raw')
plt.plot(x, gaus(x, *popt), 'o', markersize = 4, label='Gaussian fit')
plt.xlabel('X')
plt.ylabel('Y')
plt.legend()
fig.savefig('gaussian_fit.pdf')

enter image description here

最佳答案

您不能为此使用 scikit-learn,因为您处理的不是要估计其分布的一组样本。您当然可以将曲线转换为 PDF,对其进行采样,然后尝试使用高斯混合模型对其进行拟合,但这对我来说似乎有点矫枉过正。

这是一个使用简单的最小二乘曲线拟合的解决方案。为了让它工作,我必须删除背景,即忽略所有带有 y < 5 的数据点。 ,并为 leastsq 提供了一个很好的起始向量,可以通过数据图进行估算。

寻找起始向量

用最小二乘法求出的参数向量就是向量

params = [c1, mu1, sigma1, c2, mu2, sigma2]

在这里,c1c2是两个高斯的比例因子,即它们的高度,mu1mu2是平均值,即峰的水平位置和 sigma1sigma2确定高斯宽度的标准偏差。为了找到起始向量,我只查看了数据图并估计了两个峰的高度(分别为 = c1c2)及其水平位置(分别为 = mu1mu1) . sigma1sigma2简单地设置为 1.0 .

代码

from sklearn import mixture
import matplotlib.pyplot
import matplotlib.mlab
import numpy as np
from pylab import *
from scipy.optimize import leastsq

data = np.genfromtxt('gaussian_fit.dat', skiprows = 1)
x = data[:, 0]
y = data[:, 1]

def double_gaussian( x, params ):
(c1, mu1, sigma1, c2, mu2, sigma2) = params
res = c1 * np.exp( - (x - mu1)**2.0 / (2.0 * sigma1**2.0) ) \
+ c2 * np.exp( - (x - mu2)**2.0 / (2.0 * sigma2**2.0) )
return res

def double_gaussian_fit( params ):
fit = double_gaussian( x, params )
return (fit - y_proc)

# Remove background.
y_proc = np.copy(y)
y_proc[y_proc < 5] = 0.0

# Least squares fit. Starting values found by inspection.
fit = leastsq( double_gaussian_fit, [13.0,-13.0,1.0,60.0,3.0,1.0] )
plot( x, y, c='b' )
plot( x, double_gaussian( x, fit[0] ), c='r' )

关于python - 如何在 Python 中拟合双高斯分布?,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/33178790/

25 4 0
Copyright 2021 - 2024 cfsdn All Rights Reserved 蜀ICP备2022000587号
广告合作:1813099741@qq.com 6ren.com