如何在Python中擬合多個高斯曲線到質譜數據？

Question

我想在Python中將多條高斯曲線擬合爲質譜數據。現在，我一次將數據擬合成一個高斯 - 一次一個範圍。

有沒有更簡化的方法來做到這一點？有沒有辦法通過循環運行數據來在每個峯值處繪製高斯？我猜想有一個更好的方法，但我已經梳理了互聯網。

我的兩個高斯圖如下所示。

我的示例數據，可以發現：http://txt.do/dooxv

下面是我當前的代碼：

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

from scipy.interpolate import interp1d 

RGAdata = np.loadtxt("/Users/ilenemitchell/Desktop/RGAscan.txt", skiprows=14) 
RGAdata=RGAdata.transpose() 

x=RGAdata[0] 
y=RGAdata[1] 

# graph labels 
plt.ylabel('ion current') 
plt.xlabel('mass/charge ratio') 
plt.xticks(np.arange(min(RGAdata[0]), max(RGAdata[0])+2, 2.0)) 
plt.ylim([10**-12.5, 10**-9]) 
plt.title('RGA Data Jul 25, 2017') 

plt.semilogy(x, y,'b') 

#fitting a guassian to a peak 

def gauss(x, a, mu, sig): 
return a*np.exp(-(x-mu)**2/(2*sig**2)) 


fitx=x[(x>40)*(x<43)] 
fity=y[(x>40)*(x<43)] 
mu=np.sum(fitx*fity)/np.sum(fity) 
sig=np.sqrt(np.sum(fity*(fitx-mu)**2)/np.sum(fity)) 

print (mu, sig, max(fity)) 

popt, pcov = opt.curve_fit(gauss, fitx, fity, p0=[max(fity),mu, sig]) 
plt.semilogy(x, gauss(x, popt[0],popt[1],popt[2]), 'r-', label='fit') 

#second guassian 

fitx2=x[(x>26)*(x<31)] 
fity2=y[(x>26)*(x<31)] 
mu=np.sum(fitx2*fity2)/np.sum(fity2) 
sig=np.sqrt(np.sum(fity2*(fitx2-mu)**2)/np.sum(fity2)) 

print (mu, sig, max(fity2)) 

popt2, pcov2 = opt.curve_fit(gauss, fitx2, fity2, p0=[max(fity2),mu, sig]) 
plt.semilogy(x, gauss(x, popt2[0],popt2[1],popt2[2]), 'm', label='fit2') 

plt.show() 

Answer 1

這裏有一個數據集，讓你開始識別峯的一些示例代碼。你可以找到所有例子的鏈接here。

import numpy as np 
import peakutils 
cb = np.array([-0.010223, ... ]) 
indexes = peakutils.indexes(cb, thres=0.02/max(cb), min_dist=100) 
# [ 333 693 1234 1600] 

interpolatedIndexes = peakutils.interpolate(range(0, len(cb)), cb, ind=indexes) 
# [ 332.61234263 694.94831376 1231.92840845 1600.52446335] 

Answer 2

除了亞歷克斯·F公司的答案，你需要確定高峯和分析周圍的環境來識別xmin和xmax值。

如果你這樣做，你可以使用這個範圍內稍微重構代碼和循環繪製的所有相關數據

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

from scipy.interpolate import interp1d 

def _gauss(x, a, mu, sig): 
    return a*np.exp(-(x-mu)**2/(2*sig**2)) 

def gauss(x, y, xmin, xmax): 
    fitx = x[(x>xmin)*(x<xmax)] 
    fity = y[(x>xmin)*(x<xmax)] 
    mu = np.sum(fitx*fity)/np.sum(fity) 
    sig = np.sqrt(np.sum(fity*(fitx-mu)**2)/np.sum(fity)) 

    print (mu, sig, max(fity)) 

    popt, pcov = opt.curve_fit(_gauss, fitx, fity, p0=[max(fity), mu, sig]) 
    return _gauss(x, popt[0], popt[1], popt[2]) 

# Load data and define x - y 
RGAdata = np.loadtxt("/Users/ilenemitchell/Desktop/RGAscan.txt", skiprows=14) 
x, y = RGAdata.T 

# Create the plot 
fig, ax = plt.subplots() 
ax.semilogy(x, y, 'b') 

# Plot the Gaussian's between xmin and xmax 
for xmin, xmax in [(40, 43), (26, 31)]: 
    yG = gauss(x, y, xmin, xmax) 
    ax.semilogy(x, yG) 

# Prettify the graph 
ax.set_xlabel("mass/charge ratio") 
ax.set_ylabel("ion current") 
ax.set_xticks(np.arange(min(x), max(x)+2, 2.0)) 
ax.set_ylim([10**-12.5, 10**-9]) 
ax.set_title("RGA Data Jul 25, 2017") 

plt.show() 

Answer 3

您可能會發現lmfit模塊（https://lmfit.github.io/lmfit-py/）有幫助。這提供了一個預先構建的GaussianModel類，用於將峯值擬合爲單個高斯，並支持向複合模型中添加多個模型（不一定是高斯，還包括其他峯值模型和其他可能對背景有用的函數）立即適合。

Lmfit支持固定或給予一定範圍的一些參數，這樣就可以建立一個模型，高斯的固定位置的總和，限制值的質心與一定範圍內變化（這樣峯不能混淆） 。另外，您可以對參數值施加簡單的數學約束，以便您可能要求所有峯寬都是相同的大小（或以某種簡單形式相關）。

特別是，你可以看看https://lmfit.github.io/lmfit-py/builtin_models.html#example-3-fitting-multiple-peaks-and-using-prefixes的一個例子，使用2個高斯和一個背景函數擬合。我發現scipy.signal.find_peaks_cwt是非常好的。