Solución:
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:
Explicación
Necesita buenos valores iniciales de modo que el curve_fit la función converge 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, divida 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 sobre la definición de la media (incluida la respuesta del desarrollador)
Dado que a los revisores no les gustó mi edición en el 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.
Tu definición regresa:
>>> sum(x * y)
125
Devuelve la definición del desarrollador:
>>> sum(x * y) / len(x)
12.5 #for Python 3.x
La media aritmética ponderada:
>>> sum(x * y) / sum(y)
5.0
De manera similar, puede 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, debe 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)
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 max (y), 5 en el ejemplo, en lugar de 1.