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

Question

我正在尝试去混合低分辨率光谱的发射线以获得高斯分量。该图代表我正在使用的数据类型:

经过一番搜索，我发现的唯一选择是应用 kmpfit 包中的 gauest 函数( http://www.astro.rug.nl/software/kapteyn/kmpfittutorial.html#gauest )。我复制了他们的示例，但无法使其工作。

我想知道是否有人可以为我提供任何替代方法来执行此操作或如何更正我的代码:

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

def 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, y

def 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_Maxima

def 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 Output

def Residuals_GaussianMixture(p, x, y, ZeroLevel):    
    return GaussianMixture_Model(p, x, ZeroLevel) - y

Wave, Flux  = CurveData()

Wave_Maxima, Flux_Maxima = FindMaxima(Wave, Flux)

EmLines_Number = len(Wave_Maxima)

ContinuumLevel = 1.64191e-16

# Define initial values
p_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 p1
Axis1.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()   

Answer 1

典型的简单拟合方式:

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

from scipy import optimize

p0 = [1e-15,3968,2,1e-15,3972,2]
p1,conv = optimize.leastsq(res,p0[:],args=(x,y))

plot(x,y,'+') # data
#fitted function
plot(arange(3962,3976,0.1),model(p1,arange(3962,3976,0.1)),'-')

其中 p0 是您最初的猜测。从表面上看，您可能想要使用洛伦兹函数......

如果您使用 full_output=True，您将获得有关拟合的所有信息。另请查看 scipy.optimize 中的 curve_fit 和 fmin* 函数。周围有很多包装器，但通常，像这里一样，直接使用它们会更容易。