lorentzian - lmfit models python

Ajuste gaussiano para Python

Question

Estoy tratando de encajar un gaussiano para mis datos (que ya es un gaussiano aproximado). Ya tomé el consejo de los que están aquí y probé curve_fit y leastsq pero creo que me falta algo más fundamental (porque no tengo idea de cómo usar el comando). Aquí hay un vistazo al guión que tengo hasta ahora.

import pylab as plb
import matplotlib.pyplot as plt

# Read in data -- first 2 rows are header in this example. 
data = plb.loadtxt(''part 2.csv'', skiprows=2, delimiter='','')

x = data[:,2]
y = data[:,3]
mean = sum(x*y)
sigma = sum(y*(x - mean)**2)

def gauss_function(x, a, x0, sigma):
    return a*np.exp(-(x-x0)**2/(2*sigma**2))
popt, pcov = curve_fit(gauss_function, x, y, p0 = [1, mean, sigma])
plt.plot(x, gauss_function(x, *popt), label=''fit'')

# plot data

plt.plot(x, y,''b'')

# Add some axis labels

plt.legend()
plt.title(''Fig. 3 - Fit for Time Constant'')
plt.xlabel(''Time (s)'')
plt.ylabel(''Voltage (V)'')
plt.show()

Lo que obtengo de esto es una forma gaussiana, que es mi información original, y una línea recta horizontal.

Además, me gustaría trazar mi gráfica usando puntos, en lugar de tenerlos conectados. Cualquier entrada es apreciada!

Answer 1

Explicación

Necesita buenos valores de inicio, de modo que la función curve_fit converja en valores "buenos". Realmente no puedo decir por qué su ajuste no convergió (aunque la definición de su media es extraña, verifique a continuación), pero le daré una estrategia que funciona para funciones gaussianas no normalizadas como la suya.

Ejemplo

Los parámetros estimados deben estar cerca de los valores finales (use la media aritmética ponderada , dividir por la suma de todos los valores):

import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
import numpy as np

x = np.arange(10)
y = np.array([0, 1, 2, 3, 4, 5, 4, 3, 2, 1])

# weighted arithmetic mean (corrected - check the section below)
mean = sum(x * y) / sum(y)
sigma = np.sqrt(sum(y * (x - mean)**2) / sum(y))

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

popt,pcov = curve_fit(Gauss, x, y, p0=[max(y), mean, sigma])

plt.plot(x, y, ''b+:'', label=''data'')
plt.plot(x, Gauss(x, *popt), ''r-'', label=''fit'')
plt.legend()
plt.title(''Fig. 3 - Fit for Time Constant'')
plt.xlabel(''Time (s)'')
plt.ylabel(''Voltage (V)'')
plt.show()

Personalmente prefiero usar numpy.

Comente la definición de la media (incluida la respuesta del desarrollador)

Como a los revisores no les gustó mi edición en # código del desarrollador , explicaré en qué caso sugeriría un código mejorado. La media del desarrollador no corresponde a una de las definiciones normales de la media.

Su definición devuelve:

>>> sum(x * y)
125

La definición del desarrollador devuelve:

>>> sum(x * y) / len(x)
12.5 #for Python 3.x

La media aritmética ponderada:

>>> sum(x * y) / sum(y)
5.0

Del mismo modo puedes comparar las definiciones de desviación estándar ( sigma ). Compare con la figura del ajuste resultante:

Comentario para usuarios de Python 2.x

En Python 2.x, además, debes usar la nueva división para no encontrar resultados extraños o convertir los números antes de la división explícitamente:

from __future__ import division

o por ejemplo

sum(x * y) * 1. / sum(y)

Answer 2

Después de perder horas tratando de encontrar mi error, el problema es su fórmula:

sigma = suma (y * (x-mean) ** 2) / n esto es incorrecto, la fórmula correcta es la casilla cuadrada de esto !;

sqrt (suma (y * (x-mean) ** 2) / n)

Espero que esto ayude

Answer 3

Hay otra forma de realizar el ajuste, que es mediante el uso del paquete ''lmfit''. Básicamente usa el cuve_fit pero es mucho mejor en ajuste y también ofrece un ajuste complejo. Las instrucciones detalladas paso a paso se dan en el siguiente enlace. http://cars9.uchicago.edu/software/python/lmfit/model.html#model.best_fit

Answer 4

Obtienes una línea recta horizontal porque no converge.

Se logra una mejor convergencia si el primer parámetro del ajuste (p0) se pone como máximo (y), 5 en el ejemplo, en lugar de 1.

Answer 5

Aquí está el código corregido:

import pylab as plb
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
from scipy import asarray as ar,exp

x = ar(range(10))
y = ar([0,1,2,3,4,5,4,3,2,1])

n = len(x)                          #the number of data
mean = sum(x*y)/n                   #note this correction
sigma = sum(y*(x-mean)**2)/n        #note this correction

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

popt,pcov = curve_fit(gaus,x,y,p0=[1,mean,sigma])

plt.plot(x,y,''b+:'',label=''data'')
plt.plot(x,gaus(x,*popt),''ro:'',label=''fit'')
plt.legend()
plt.title(''Fig. 3 - Fit for Time Constant'')
plt.xlabel(''Time (s)'')
plt.ylabel(''Voltage (V)'')
plt.show()

resultado:

Answer 6

sigma = sum(y*(x - mean)**2)

debiera ser

sigma = np.sqrt(sum(y*(x - mean)**2))