使用python分離曲線的高斯分量

Question

我想解耦低分辨率光譜的發射線以得到高斯分量。此圖代表我使用的數據類型：

搜索了一下後，我發現的唯一的選擇是從kmpfit包（http://www.astro.rug.nl/software/kapteyn/kmpfittutorial.html#gauest）的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 *函數。這些周圍有很多包裝，但是經常像這裏一樣，直接使用它們更容易。