【发布时间】:2019-11-15 10:21:16
【问题描述】:
我正在尝试概括一些代码,以便能够在单个数据集中拟合多个(n 从 1 到 >10)高斯曲线/峰值。
使用 Scipy Optimize Curve_fit 当我对 1-3 个高斯函数进行硬编码时,我可以得到很好的拟合,并且我已经成功地生成了对于一般的、任意数量的高斯函数运行没有错误的函数。但是,输出拟合很差。尽管提供的输入参数与用于生成“原始”数据的输入参数相同 - 即最佳情况。
此外,在某些时候可能需要从简单的高斯函数修改特定函数的可能性非零,但现在应该没问题。
下面是我的代码示例,输出图如下所示。
import numpy as np
import pandas as pd
import scipy
import scipy.optimize
import matplotlib.pyplot as plt
from matplotlib import gridspec
amp1 = 1
cen1 = 1
sigma1 = 0.05
df=pd.DataFrame(index=np.linspace(0,10,num=1000),columns=['int'])
def _ngaussian(x, amps,cens,sigmas):
fn = 0
if len(amps)== len(cens)== len(sigmas):
for i in range(len(amps)):
fn = fn+amps[i]*(1/(sigmas[i]*(np.sqrt(2*np.pi))))*\
(np.exp((-1.0/2.0)*(((x-cens[i])/sigmas[i])**2)))
else:
print('Your inputs have unequal lengths')
return fn
amps = [1,1.1,0.9]
cens = [1,2,1.7]
sigmas=[0.05]*3
popt_peaks = [amps,cens,sigmas]
df['peaks'] = _ngaussian(df.index, *popt_peaks)
# Optionally adding noise to the raw data
#noise = np.random.normal(0,0.1,len(df['peaks']))
#df['peaks'] = df['peaks']+noise
def wrapper_fit_func(x, *args):
N = len(args)
a, b, c = list(args[0][:N]),list(args[0][N:N*2]),list(args[0][2*N:3*N])
return _ngaussian(x, a, b, c)
def unwrapper_fit_func(x, *args):
N = int(len(args)/3)
a, b, c = list(args[:N]),list(args[N:N*2]),list(args[2*N:3*N])
return _ngaussian(x, a, b, c)
popt_fitpeaks, pcov_fitpeaks = scipy.optimize.curve_fit(lambda x, *popt_peaks: wrapper_fit_func(x, popt_peaks),
df.index, df['peaks'], p0=popt_peaks,
method='lm')
df['peaks_fit'] = unwrapper_fit_func(df.index, *popt_fitpeaks)
fig = plt.figure(figsize=(8,8))
gs = gridspec.GridSpec(1,1)
ax1 = fig.add_subplot(gs[0])
ax1.set_xlim(0,3)
ax1.plot(df.index, df['peaks'], "b",label='ideal data')
ax1.plot(df.index, df['peaks_fit'], "g",label='fit data')
ax1.legend(loc='upper right')
如果您有兴趣,可以参考分析化学、核磁共振 (NMR) 和傅里叶变换离子回旋共振质谱 (FTICR MS) 信号处理。
【问题讨论】:
-
我建议首先拟合单个峰,例如使用“index=np.linspace(0,1.25,num=1000)”。在那之后工作,然后尝试两个峰值,然后最后全部三个。使用您发布的代码,我个人不太适合单峰。
标签: python curve-fitting scipy-optimize