python - 如何使用 LMfit 将曲线拟合到双高斯/偏斜高斯
问题描述
我有一堆代码可以将质谱峰从光谱中分离出来,并将峰的值放入两个列表中。Gaussian_x 和 gaussian_x。我现在需要拟合这条曲线,最好使用 Levenberg-Marquardt 或类似的算法,这样我就可以在拟合后计算峰面积。
对于双高斯峰(在峰的前半部分和后半部分具有不同 sigma 值的峰)并没有太大帮助,因此我正在努力寻找使其工作的方法。
yvals = np.asarray(gaussian_y)
model = SkewedGaussianModel()
# set initial parameter values
params = model.make_params(amplitude=a, center=b, sigma=c, gamma=d)
# adjust parameters to best fit data.
result = model.fit(yvals, params, x=xvals)
print(result.fit_report())
plt.plot(xvals, yvals)
plt.plot(xvals, result.init_fit)
plt.plot(xvals, result.best_fit)
plt.show()
当我绘制这个图时,在 y=0 处只有一条直线,然后是两个列表的图。我不知道为什么 yvalue 没有向 best_fit 和 init_fit 注册。
这是完整的代码:
import pymzml
import numpy as np
import matplotlib.pyplot as plt
from lmfit.models import SkewedGaussianModel
import filepath
TARGET_MASS = 152
MASS_WIDTH = 1
PEAK_CENTER = 0.9
PEAK_WIDTH = 0.1
gaussian_x = []
gaussian_y = []
a = 1
b = 400
c = 100
d = 280
def main(in_path):
x_array = []
y_array = []
run = pymzml.run.Reader(in_path)
for spectrum in run:
if spectrum.ms_level == 1:
mass_array = []
intensity_array = []
target_array = []
# change tuple into array.
my_tuple = spectrum.peaks("raw")
my_array = np.asarray(my_tuple)
# append the first element of each index to mass array and second element to intensity array.
for i in my_array:
mass_array.append(i[0])
intensity_array.append(i[1])
# for index and value in mass_array, absolute value of mass - target mass taken.
# append the abs_mass values to a separate array.
for i, mass in enumerate(mass_array):
abs_mass = abs(mass - TARGET_MASS)
target_array.append(abs_mass)
# if the absolute mass values are less than selected mass_width, append the value in the intensity_array
# at that same index same index
if abs_mass < MASS_WIDTH:
y_array.append(intensity_array[i])
x_array.append(spectrum.scan_time_in_minutes())
new_x = []
new_y = []
# loop through x values
for i, x1 in enumerate(x_array):
# temporary store for where there are multiple y values for the same x value
y_temp = []
# ensure we only run through a search for each UNIQUE x value once
if x1 in new_x:
# moves onto the next for loop interation without executing the below
continue
# add unique x value to new_x
new_x.append(x1)
# loop through the x values again, looking for where there are duplicates
for j, x2 in enumerate(x_array):
if x1 == x2:
# where we find a duplicate entry, add the value of y to a temporary array
y_temp.append(y_array[j])
# after accumulating all values of y for the same x value,
# add it to the new_y array
new_y.append(max(y_temp))
lower_bound = PEAK_CENTER - (PEAK_WIDTH / 2)
upper_bound = PEAK_CENTER + (PEAK_WIDTH / 2)
# for each index and value in retention time array
for i, x in enumerate(new_x):
# if x greater than the lower bound and smaller than the upper bound then append x and y value
# to new lists
if lower_bound < x < upper_bound:
gaussian_x.append(x)
gaussian_y.append(new_y[i])
# if x greater than upper bound, stop function
if x > upper_bound:
break
xvals = np.asarray(gaussian_x)
yvals = np.asarray(gaussian_y)
model = SkewedGaussianModel()
# set initial parameter values
params = model.make_params(amplitude=a, center=b, sigma=c, gamma=d)
# adjust parameters to best fit data.
result = model.fit(yvals, params, x=xvals)
print(result.fit_report())
plt.plot(xvals, yvals)
plt.plot(xvals, result.init_fit)
plt.plot(xvals, result.best_fit)
plt.show()
if __name__ == "__main__":
main(in_path)
适合报告说:
[[Model]]
Model(skewed_gaussian)
[[Fit Statistics]]
# fitting method = leastsq
# function evals = 5
# data points = 22
# variables = 4
chi-square = 1.6213e+10
reduced chi-square = 9.0074e+08
Akaike info crit = 457.197124
Bayesian info crit = 461.561294
## Warning: uncertainties could not be estimated:
amplitude: at initial value
center: at initial value
sigma: at initial value
gamma: at initial value
[[Variables]]
amplitude: 1.00000000 (init = 1)
center: 400.000000 (init = 400)
sigma: 100.000000 (init = 100)
gamma: 280.000000 (init = 280)
height: 0.00398942 == '0.3989423*amplitude/max(2.220446049250313e-16, sigma)'
fwhm: 235.482000 == '2.3548200*sigma'
解决方案
似乎您适合的初始参数不够好。请参阅以下讨论以更好地猜测您的参数; https://stackoverflow.com/a/19207683/13131172
推荐阅读
- gdb - 带有 Eclipse Oxygen 3a 的远程 gdb
- php - 在 PHP 中转义一个充满“和”的字符串
- php - 如果存在则更新行或插入
- angular - 尝试使用带有自定义 tslint 规则的 typescript 库,但在 ng lint 时出现“无法在模块外使用 import 语句”错误
- java - 基于值数组的绘画
- javascript - Java 特殊字符验证
- php - 我如何在 laravel 中修复此错误“1005 无法创建表`englishcollage`.`role_user`(错误号:150“外键约束格式不正确”
- angular - 谁/如何/何时使用 SSR 在 Angular 中创建 app.browser.module.ts 文件?
- python - 尝试使用 Python 列出 IBM Cloud Object Storage 中的对象时出现“指定的存储桶不存在”
- sql - PL/SQL 执行中间更新查询,给出 SQLCODE -932