Encontrando el ancho máximo medio máximo de un pico.

He estado tratando de calcular la mitad máxima máxima (FWHM) del ancho azul (ver imagen). El pico verde y el pico magenta combinados forman el pico azul. He estado usando la siguiente ecuación para encontrar el FWHM de los picos verde y magenta: fwhm = 2*np.sqrt(2*(math.log(2)))*sd donde sd = desviación estándar. Creé los picos verdes y magenta y conozco la desviación estándar, por eso puedo usar esa ecuación.

Creé los picos verdes y magenta usando el siguiente código:

 def make_norm_dist(self, x, mean, sd): import numpy as np norm = [] for i in range(x.size): norm += [1.0/(sd*np.sqrt(2*np.pi))*np.exp(-(x[i] - mean)**2/(2*sd**2))] return np.array(norm) 

Si no supiera que el pico azul estaba formado por dos picos y solo tenía el pico azul en mis datos, ¿cómo encontraría el FWHM?

He estado usando este código para encontrar la cima del pico:

 peak_top = 0.0e-1000 for i in x_axis: if i > peak_top: peak_top = i 

Podría dividir el peak_top por 2 para encontrar la mitad de la altura y luego tratar de encontrar los valores de y correspondientes a la mitad de la altura, pero luego me encontraría con problemas si no hubiera valores de x que coincidieran exactamente con la altura de la mitad.

Estoy bastante seguro de que hay una solución más elegante a la que estoy intentando.

Related of "Encontrando el ancho máximo medio máximo de un pico."

Puede usar la spline para ajustar la [curva azul – pico / 2], y luego encontrar sus raíces:

 import numpy as np from scipy.interpolate import UnivariateSpline def make_norm_dist(x, mean, sd): return 1.0/(sd*np.sqrt(2*np.pi))*np.exp(-(x - mean)**2/(2*sd**2)) x = np.linspace(10, 110, 1000) green = make_norm_dist(x, 50, 10) pink = make_norm_dist(x, 60, 10) blue = green + pink # create a spline of x and blue-np.max(blue)/2 spline = UnivariateSpline(x, blue-np.max(blue)/2, s=0) r1, r2 = spline.roots() # find the roots import pylab as pl pl.plot(x, blue) pl.axvspan(r1, r2, facecolor='g', alpha=0.5) pl.show() 

Aquí está el resultado:

introduzca la descripción de la imagen aquí

Esto me funcionó en iPython (rápido y sucio, se puede reducir a 3 líneas):

 def FWHM(X,Y): half_max = max(Y) / 2. #find when function crosses line half_max (when sign of diff flips) #take the 'derivative' of signum(half_max - Y[]) d = sign(half_max - array(Y[0:-1])) - sign(half_max - array(Y[1:])) #plot(X[0:len(d)],d) #if you are interested #find the left and right most indexes left_idx = find(d > 0)[0] right_idx = find(d < 0)[-1] return X[right_idx] - X[left_idx] #return the difference (full width) 

Se pueden hacer algunas adiciones para que la resolución sea más precisa, pero en el límite de que hay muchas muestras a lo largo del eje X y los datos no son demasiado ruidosos, esto funciona muy bien.

Incluso cuando los datos no son gaussianos y un poco ruidosos, funcionó para mí (solo tomo la primera y la última vez que Half Max cruza los datos).

Si sus datos tienen ruido (y siempre lo hace en el mundo real), una solución más robusta sería ajustar un gaussiano a los datos y extraer FWHM de eso:

 import numpy as np import scipy.optimize as opt def gauss(x, p): # p[0]==mean, p[1]==stdev return 1.0/(p[1]*np.sqrt(2*np.pi))*np.exp(-(xp[0])**2/(2*p[1]**2)) # Create some sample data known_param = np.array([2.0, .7]) xmin,xmax = -1.0, 5.0 N = 1000 X = np.linspace(xmin,xmax,N) Y = gauss(X, known_param) # Add some noise Y += .10*np.random.random(N) # Renormalize to a proper PDF Y /= ((xmax-xmin)/N)*Y.sum() # Fit a guassian p0 = [0,1] # Inital guess is a normal distribution errfunc = lambda p, x, y: gauss(x, p) - y # Distance to the target function p1, success = opt.leastsq(errfunc, p0[:], args=(X, Y)) fit_mu, fit_stdev = p1 FWHM = 2*np.sqrt(2*np.log(2))*fit_stdev print "FWHM", FWHM 

introduzca la descripción de la imagen aquí

La imagen trazada puede ser generada por:

 from pylab import * plot(X,Y) plot(X, gauss(X,p1),lw=3,alpha=.5, color='r') axvspan(fit_mu-FWHM/2, fit_mu+FWHM/2, facecolor='g', alpha=0.5) show() 

Una mejor aproximación filtraría los datos ruidosos por debajo de un umbral dado antes del ajuste.

Aquí hay una pequeña función que usa el enfoque de spline.

 from scipy.interpolate import splrep, sproot, splev class MultiplePeaks(Exception): pass class NoPeaksFound(Exception): pass def fwhm(x, y, k=10): """ Determine full-with-half-maximum of a peaked set of points, x and y. Assumes that there is only one peak present in the datasset. The function uses a spline interpolation of order k. """ half_max = amax(y)/2.0 s = splrep(x, y - half_max, k=k) roots = sproot(s) if len(roots) > 2: raise MultiplePeaks("The dataset appears to have multiple peaks, and " "thus the FWHM can't be determined.") elif len(roots) < 2: raise NoPeaksFound("No proper peaks were found in the data set; likely " "the dataset is flat (eg all zeros).") else: return abs(roots[1] - roots[0])