Ajuster un pic gaussien à mesuré

Question

J'ai un petit pic de spectre et j'essaie de lui associer une fonction gaussienne. J'ai cherché un exemple en ligne et mélangé le code avec celui que j'ai fait.

wveleng=[ 639.188 639.454 639.719 639.985 640.25 640.516 640.781 641.046 
     641.312 641.577] 
    counts=[ 778. 1613.8 12977.4 32990. 33165.2 13171. 2067.2 900.8 
     788.8 747.8] 

mon premier code est le suivant

def gaus(x,a,mu,sigma): 
    return a*exp(-(x-mu)**2/(2*sigma**2)) 

    a=ydata.max() 
x0=ydata.mean() 
sigm=ydata.std() 

mean = sum(ydata*xdata)/len(ydata) 
sigma = np.sqrt(sum(ydata*(xdata-mean)**2)/len(ydata)) 

#print(ydata.max()) 
popt, pcov = curve_fit(Gauss, xdata,ydata,maxfev=991,p0=[a,x0,sigm])  
#gmodel = Model(Gauss) 
#result = gmodel.fit(ydata, x=xdata, a=ydata.max(),x0=ydata.mean(),sigm=ydata.std()) 
print(popt) 
#plt.scatter(xdata,ydata,label='data points') 
#plt.plot(xdata, result.best_fit, 'r-') 
#popt, pcov = curve_fit(gauss, xdata, ydata,p0=[ydata.max(), ydata.mean(), ydata.std()]) 
xx = np.linspace(639,642, 10) 
plt.plot(xx, gauss(xdata, *popt), 'r-', label='fit') 

avec terrain i obtenir ce qui suit.

je pense qu'il a à faire avec des paramètres estimation initiale

un second code que je trouve plus compact et convient mieux pour moi.

def gauss(x, a, x0, sigma): 
    return a * np.exp(-(x - x0) ** 2/(2 * sigma ** 2)) 

ydata = np.array([778.,1613.8,12977.4,32990.,33165.2,13171.,2067.2,900.8,788.8,747.8]) 

xx = np.arange(639,642, 100) 
xdata=np.array([639.188,639.454,639.719,639.985,640.250,640.516,640.781,641.046,641.312,641.577]) 


#plt.plot(xdata, ydata, 'bo', label='data') 
def Gauss(x, a, x0, sigm): 
    return a * np.exp(-(x - x0)**2/(2 * sigm**2)) 

gmodel = Model(Gauss) 
result = gmodel.fit(ydata, x=xdata, a=ydata.max(),x0=ydata.mean(),sigm=ydata.std()) 
plt.scatter(xdata,ydata,label='data points') 
plt.plot(xdata, result.best_fit, 'r-') 

Je reçois exactement la même forme que la première méthode. Y a-t-il un moyen d'ajuster plus de points que les données elles-mêmes?

Answer 1

Je pense que vous êtes vraiment proche, même si je dois admettre que je ne comprends pas ce que xx est censé représenter. Vous voulez absolument que les données soient correctes (ydata) et que la variable indépendante (xdata) ait la même longueur.

Je pense que le principal problème que vous utilisez en est maintenant que vos estimations initiales ne sont pas très bien, et que vous obtiendrez de bons résultats avec

result = gmodel.fit(ydata, x=xdata, a=ydata.max(), x0=xdata.mean(), sigma=xdata.std()) 

(avec xdata au lieu de ydata contrôler les valeurs initiales pour x0 et sigma).

Peut-être encore mieux serait d'ajouter certains contrôles à la plage de valeurs de paramètres, comme

params = gmodel.make_params(a=ydata.max(), 
          x0=xdata.mean(), 
          sigma=xdata.std()) 
params['x0'].min = min(xdata) 
params['x0'].max = max(xdata) 
params['sigma'].max = 5 
result = gmodel.fit(ydata, params, x=xdata) 

Enfin, en utilisant les modèles intégrés comme GaussianModel fera rapport à la fois sigma et fwhm et aussi amplitude (c'est-à-dire, intégrale du Gaussien) et height (la valeur maximale que le Gaussien prendrait).Donc, ce script:

import numpy as np 
from lmfit.models import GaussianModel 
import matplotlib.pyplot as plt 

ydata = np.array([778.,1613.8,12977.4,32990.,33165.2,13171.,2067.2,900.8,788.8,747.8]) 
xdata = np.array([639.188,639.454,639.719,639.985,640.250,640.516,640.781,641.046,641.312,641.577]) 

gmodel = GaussianModel() 
params = gmodel.make_params(amplitude=ydata.max(), 
          center=xdata.mean(), 
          sigma=xdata.std()) 
result = gmodel.fit(ydata, params, x=xdata) 

print(result.fit_report()) 
plt.scatter(xdata,ydata,label='data points') 
plt.plot(xdata, result.best_fit, 'r-') 
plt.show() 

imprimera

[[Model]] 
    Model(gaussian) 
[[Fit Statistics]] 
    # function evals = 27 
    # data points  = 10 
    # variables  = 3 
    chi-square   = 2360971.771 
    reduced chi-square = 337281.682 
    Akaike info crit = 129.720 
    Bayesian info crit = 130.628 
[[Variables]] 
    sigma:  0.27525232 +/- 0.004505 (1.64%) (init= 0.7623906) 
    center:  640.119396 +/- 0.004490 (0.00%) (init= 640.3828) 
    amplitude: 25633.2702 +/- 362.5571 (1.41%) (init= 33165.2) 
    fwhm:  0.64816968 +/- 0.010608 (1.64%) == '2.3548200*sigma' 
    height:  37152.0777 +/- 525.0717 (1.41%) == '0.3989423*amplitude/max(1.e-15, sigma)' 
[[Correlations]] (unreported correlations are < 0.100) 
    C(sigma, amplitude)   = 0.579 

et donner une assez bonne mine en forme. Pour la pratique avancée, je recommande d'essayer d'ajouter un ConstantModel() pour donner un décalage d'arrière-plan. Eh bien, et collecter plus de points de données;).

Answer 2

scipy.integrate.quad ne fait pas de convolution, comme vous semblez l'espérer. quad(function, lower_bound, upper_bound)[0] retournera une seule valeur pour l'intégrale de la fonction entre les bornes.

OTOH, curve_fit(func, ...) a besoin d'un tableau de valeurs pour la fonction du modèle, et il se plaint qu'il a obtenu un flottant, pas un ndarray.

Peut-être que vous aviez l'intention de faire curve_fit(vfunc, ...)?

Vous pourriez trouver lmfit (https://lmfit.github.io/lmfit-py/) utile. Il comporte des outils pratiques de haut niveau adaptés aux courbes et des fonctions de modèle simples comme le gaussien sont intégrées. Il dispose également d'un mécanisme permettant d'ajuster un modèle qui est la somme ou le produit de deux fonctions, et peut même créer un modèle composé de la convolution de deux fonctions. Par exemple de ceci, voir les exemples décrits dans https://lmfit.github.io/lmfit-py/model.html#composite-models-adding-or-multiplying-models.