如何在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是非常好的。