使用python分离曲线的高斯分量

本文介绍了使用python分离曲线的高斯分量的处理方法，对大家解决问题具有一定的参考价值，需要的朋友们下面随着小编来一起学习吧！ 问题描述 我正在尝试低分辨率光谱的发射线，以便获得高斯分量。此图表示我正在使用的数据类型：I am trying to deblend the emission lines of low resolution spectrum in order to get the gaussian components. This plot represents the kind of data I am using: 经过一番搜索，我发现的唯一选择是从kmpfit包中应用gauest函数（ http://www.astro.rug.nl/software/kapteyn/kmpfittutorial .html＃gauest ）。我已经复制了他们的示例，但我无法使它正常工作。After searching a bit, the only option I found was the application of the gauest function from the kmpfit package (http://www.astro.rug.nl/software/kapteyn/kmpfittutorial.html#gauest). I have copied their example but I cannot make it work.我想知道是否有人可以提供任何替代方法，或者如何纠正我的代码：I wonder if anyone could please offer me any alternative to do this or how to correct my code:import numpy as npimport matplotlib.pyplot as pltfrom scipy import optimizedef CurveData(): x = np.array([3963.67285156, 3964.49560547, 3965.31835938, 3966.14111328, 3966.96362305, 3967.78637695, 3968.60913086, 3969.43188477, 3970.25463867, 3971.07714844, 3971.89990234, 3972.72265625, 3973.54541016, 3974.36791992, 3975.19067383]) y = np.array([1.75001533e-16, 2.15520995e-16, 2.85030769e-16, 4.10072843e-16, 7.17558032e-16, 1.27759917e-15, 1.57074192e-15, 1.40802933e-15, 1.45038722e-15, 1.55195653e-15, 1.09280316e-15, 4.96611341e-16, 2.68777266e-16, 1.87075114e-16, 1.64335999e-16]) return x, ydef FindMaxima(xval, yval): xval = np.asarray(xval) yval = np.asarray(yval) sort_idx = np.argsort(xval) yval = yval[sort_idx] gradient = np.diff(yval) maxima = np.diff((gradient > 0).view(np.int8)) ListIndeces = np.concatenate((([0],) if gradient[0] < 0 else ()) + (np.where(maxima == -1)[0] + 1,) + (([len(yval)-1],) if gradient[-1] > 0 else ())) X_Maxima, Y_Maxima = [], [] for index in ListIndeces: X_Maxima.append(xval[index]) Y_Maxima.append(yval[index]) return X_Maxima, Y_Maximadef GaussianMixture_Model(p, x, ZeroLevel): y = 0.0 N_Comps = int(len(p) / 3) for i in range(N_Comps): A, mu, sigma = p[i*3:(i+1)*3] y += A * np.exp(-(x-mu)*(x-mu)/(2.0*sigma*sigma)) Output = y + ZeroLevel return Outputdef Residuals_GaussianMixture(p, x, y, ZeroLevel): return GaussianMixture_Model(p, x, ZeroLevel) - yWave, Flux = CurveData()Wave_Maxima, Flux_Maxima = FindMaxima(Wave, Flux)EmLines_Number = len(Wave_Maxima)ContinuumLevel = 1.64191e-16# Define initial valuesp_0 = []for i in range(EmLines_Number): p_0.append(Flux_Maxima[i]) p_0.append(Wave_Maxima[i]) p_0.append(2.0)p1, conv = optimize.leastsq(Residuals_GaussianMixture, p_0[:],args=(Wave, Flux, ContinuumLevel))Fig = plt.figure(figsize = (16, 10)) Axis1 = Fig.add_subplot(111) Axis1.plot(Wave, Flux, label='Emission line')Axis1.plot(Wave, GaussianMixture_Model(p1, Wave, ContinuumLevel), 'r', label='Fit with optimize.leastsq')print p1Axis1.plot(Wave, GaussianMixture_Model([p1[0],p1[1],p1[2]], Wave, ContinuumLevel), 'g:', label='Gaussian components')Axis1.plot(Wave, GaussianMixture_Model([p1[3],p1[4],p1[5]], Wave, ContinuumLevel), 'g:')Axis1.set_xlabel( r'Wavelength $(\AA)$',)Axis1.set_ylabel('Flux' + r'$(erg\,cm^{-2} s^{-1} \AA^{-1})$')plt.legend()plt.show() 推荐答案一种典型的简单拟合方法：A typical simplistic way to fit:def model(p,x): A,x1,sig1,B,x2,sig2 = p return A*np.exp(-(x-x1)**2/sig1**2) + B*np.exp(-(x-x2)**2/sig2**2)def res(p,x,y): return model(p,x) - yfrom scipy import optimizep0 = [1e-15,3968,2,1e-15,3972,2]p1,conv = optimize.leastsq(res,p0[:],args=(x,y))plot(x,y,'+') # data#fitted functionplot(arange(3962,3976,0.1),model(p1,arange(3962,3976,0.1)),'-')其中p0是您最初的猜测。从外观上看，您可能想使用Lorentzian函数... Where p0 is your initial guess. By the looks of things, you might want to use Lorentzian functions... 如果使用full_output = True，则会获得有关拟合的所有信息。还要在scipy.optimize中检出curve_fit和fmin *函数。周围有很多包装纸，但是通常，像这里一样，直接使用它们比较容易。If you use full_output=True, you get all kind of info about the fitting. Also check out curve_fit and the fmin* functions in scipy.optimize. There are plenty of wrappers around these around, but often, like here, it's easier to use them directly. 这篇关于使用python分离曲线的高斯分量的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持！