python - pip install numpy scipy matplotlib pandas

Gaussian Fit en un conjunto de datos ruidoso e "interesante"

Question

Tengo algunos datos (difracción de rayos X) que se ven así:

Quiero ajustar un gaussiano a este conjunto de datos para obtener el FWHM de la parte ''más ancha''. El doble pico alrededor de los 7 grados theta no es información importante y proviene de fuentes no deseadas.

Para ser más claro, quiero algo como esto (que hice en pintura :)):

Intenté escribir algo en Python con el siguiente código:

import math
from pylab import *
import numpy as np
import scipy as sp
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
data2=np.loadtxt(''FWHM.spc'')    
x2,y2=data2[:,0],data2[:,7]
plt.title(''Full Width Half Max of 002 Peak'')
plt.plot(x2, y2, color=''b'')
plt.xlabel(''$//theta$'', fontsize=10)
plt.ylabel(''Intensity'', fontsize=10)
plt.xlim([3,11])
plt.xticks(np.arange(3, 12, 1), fontsize=10)
plt.yticks(fontsize=10)

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

mean = sum(x2*y2)/sum(y2)
sigma2 = sqrt(abs(sum((x2-mean)**2*y2)/sum(y2)))
popt, pcov = curve_fit(func, x2, y2, p0 = [1, mean, sigma2])
ym = func(x2, popt[0], popt[1], popt[2])
plt.plot(x2, ym, c=''r'', label=''Best fit'')
FWHM = round(2*np.sqrt(2*np.log(2))*popt[2],4)
axvspan(popt[1]-FWHM/2, popt[1]+FWHM/2, facecolor=''g'', alpha=0.3, label=''FWHM = %s''%(FWHM))
plt.legend(fontsize=10)
plt.show()

y obtengo el siguiente resultado:

Claramente, esto está muy lejos de lo que se desea. ¿Alguien tiene algunos consejos sobre cómo puedo lograr esto?

(He incluido los datos aquí: https://justpaste.it/109qp )

Answer 1

Como se menciona en los comentarios de OP, una de las formas de restringir una señal en presencia de datos no deseados es modelarla junto con la señal deseada. Por supuesto, este enfoque es válido solo cuando hay un modelo válido disponible para esos datos contaminantes. Para los datos que proporciona, uno puede considerar un modelo compuesto que suma los siguientes componentes:

1. Una línea de base lineal, porque todos sus puntos se compensan constantemente desde cero.
2. Dos componentes Gaussianos estrechos que modelarán la característica de doble pico en la parte central de tu espectro.
3. Un componente gaussiano estrecho. Este es el que en realidad estás tratando de restringir.

Los cuatro componentes (doble pico cuenta dos veces) se pueden ajustar de forma simultánea una vez que se pasa una adivinanza inicial razonable a curve_fit :

def composite_spectrum(x, # data
                       a, b, # linear baseline
                       a1, x01, sigma1, # 1st line
                       a2, x02, sigma2, # 2nd line
                       a3, x03, sigma3 ): # 3rd line
    return (x*a + b + func(x, a1, x01, sigma1)
                    + func(x, a2, x02, sigma2)
                    + func(x, a3, x03, sigma3))

guess = [1, 200, 1000, 7, 0.05, 1000, 6.85, 0.05, 400, 7, 0.6]

popt, pcov = curve_fit(composite_spectrum, x2, y2, p0 = guess)
plt.plot(x2, composite_spectrum(x2, *popt), ''k'', label=''Total fit'')
plt.plot(x2, func(x2, *popt[-3:])+x2*popt[0]+popt[1], c=''r'', label=''Broad component'')
FWHM = round(2*np.sqrt(2*np.log(2))*popt[10],4)
plt.axvspan(popt[9]-FWHM/2, popt[9]+FWHM/2, facecolor=''g'', alpha=0.3, label=''FWHM = %s''%(FWHM))
plt.legend(fontsize=10)
plt.show()

En un caso en el que las fuentes no deseadas no se pueden modelar correctamente, la discontinuidad no deseada se puede enmascarar como lo sugiere Mad Physicist . Para el caso más simple, puedes simplemente enmascarar todo dentro de [6.5; 7.4] [6.5; 7.4] intervalo.