Как подогнать кривую к двухстороннему / наклонному гауссовскому с использованием LMfit - PullRequest
Новая
       56

Как подогнать кривую к двухстороннему / наклонному гауссовскому с использованием LMfit

0 голосов
/ 25 марта 2020

У меня есть набор кода, который изолирует пик масс-спектрометрии от спектра и поместил значения пика в два списка. Gaussian_x и gaussian_x. Теперь мне нужно подогнать эту кривую, в идеале используя Levenberg-Marquardt или аналогичный алгоритм, чтобы я мог вычислить площадь пика после его подгонки. Пики, которые имеют разные значения сигмы для первой и второй половины пика), и поэтому я изо всех сил пытаюсь найти способ заставить его работать. 

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, а затем график двух списков. Я понятия не имею, почему значения y не регистрируются в 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)

В отчете Fit указано: 

[[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'

1 Ответ

0 голосов
/ 27 марта 2020

Похоже, начальные параметры вашей посадки недостаточно хороши. Смотрите следующее обсуждение для лучшего предположения ваших параметров; { ссылка }

...