использование lmfit ExponentialGaussianModel ()

Question

Попытка подбора ExponentialGaussianModel() из lmfit, но при этом появляется следующее сообщение об ошибке: The input contains nan values

Я использую Jupyternotebook на Windows, и я новичок в Python и lmfit. Я нахожу документацию lmfit немного неясной для новичка и надеюсь найти здесь помощь. Ниже приведен мой код: я хотел бы сгенерировать экспоненциально-гауссовскую гистограмму, извлечь точки данных и попрактиковаться в подборе библиотеки lmfit. (Я хотел бы попрактиковаться в подборе и поиске наименьшего количества точек, которые бы воспроизводили параметры, используемые при создании гистограммы)

from scipy.stats import exponnorm
from lmfit.models import ExponentialGaussianModel

K2 = 1.5
r2 = exponnorm.rvs(K2, size=500, loc = 205, scale = 40)

Q           =  np.histogram(r2,500)
exp_gaus_x  =  Q[1]
exp_gaus_y  =  Q[0]

tof_x       =  exp_gaus_x[1:]
tof_y       =  exp_gaus_y

mod =  ExponentialGaussianModel()
pars = mod.guess(tof_y, x=tof_x)
out  = mod.fit(tof_y, pars, x=tof_x)
(out.fit_report(min_correl=0.25))

Я получаю ошибку, что есть входные значения nan. Я ожидал отчет, подобный показанному в руководстве.

Answer 1

Определение экспоненциального гауссиана, используемое в lmfit, взято из https://en.wikipedia.org/wiki/Exponentially_modified_Gaussian_distribution. При экспоненциальном члене

exp[center*gamma + (gamma*sigma)**2/2 - gamma*x]

это имеет тенденцию давать Infдля больших значений sigma и gamma и / или center.Я полагаю, что вы получаете такие Inf значения и что это является причиной сообщения NaN, которое вы видите.Подходящие процедуры (в Фортране) не обрабатывают NaN или Inf изящно (фактически, «вообще»).Это реальное ограничение этой конкретной модели.Вы увидите, что все примеры на странице википедии имеют значения x, гораздо более близкие к 1, чем 200, а gamma и sigma также порядка 1, а не около 50, чтоавтоматизированный guess дает.

Я думаю, что более простое определение для экспоненциально модифицированного гауссова было бы лучше для вас.С

def expgaussian(x, amplitude=1, center=0, sigma=1.0, gamma=1.0):
    """ an alternative exponentially modified Gaussian."""
    dx = center-x
    return amplitude* np.exp(gamma*dx) * erfc( dx/(np.sqrt(2)*sigma))

Вы получите приличное соответствие, хотя значение параметров немного изменилось, и вам нужно будет явно указать начальные значения, а не полагаться на процедуру guess().Они не должны быть очень близкими, просто не очень далеко.

Полный сценарий может быть:

import numpy as np
from scipy.stats import exponnorm
from scipy.special import erfc
from lmfit import Model
import matplotlib.pyplot as plt

def expgaussian(x, amplitude=1, center=0, sigma=1.0, gamma=1.0):
    """ an alternative exponentially modified Gaussian."""
    dx = center-x
    return amplitude* np.exp(gamma*dx) * erfc( dx/(np.sqrt(2)*sigma))

K2 = 1.5
r2 = exponnorm.rvs(K2, size=500, loc = 205, scale = 40)
Q           =  np.histogram(r2,500)
exp_gaus_x  =  Q[1]
exp_gaus_y  =  Q[0]
tof_x       =  exp_gaus_x[1:]
tof_y       =  exp_gaus_y

mod =  Model(expgaussian)
pars = mod.make_params(sigma=20, gamma=0.1, amplitude=2, center=250)

out  = mod.fit(tof_y, pars, x=tof_x)

print(out.fit_report())

plt.plot(tof_x, tof_y, label='data')
plt.plot(tof_x, out.best_fit, label='fit')
plt.legend()
plt.show()

, который распечатает

[[Model]]
    Model(expgaussian)
[[Fit Statistics]]
    # fitting method   = leastsq
    # function evals   = 65
    # data points      = 500
    # variables        = 4
    chi-square         = 487.546692
    reduced chi-square = 0.98295704
    Akaike info crit   = -4.61101662
    Bayesian info crit = 12.2474158
[[Variables]]
    gamma:      0.01664876 +/- 0.00239048 (14.36%) (init = 0.1)
    sigma:      39.6914678 +/- 3.90960254 (9.85%) (init = 20)
    center:     235.600396 +/- 11.8976560 (5.05%) (init = 250)
    amplitude:  3.43975318 +/- 0.15675053 (4.56%) (init = 2)
[[Correlations]] (unreported correlations are < 0.100)
    C(gamma, center)     =  0.930
    C(sigma, center)     =  0.870
    C(gamma, amplitude)  =  0.712
    C(gamma, sigma)      =  0.693
    C(center, amplitude) =  0.572
    C(sigma, amplitude)  =  0.285

и покажите сюжет вроде

Надеюсь, что поможет.