Ajuste de una distribución de Weibull usando Scipy

Estoy tratando de recrear la distribución de máxima probabilidad, ya puedo hacer esto en Matlab y R, pero ahora quiero usar scipy. En particular, me gustaría estimar los parámetros de distribución de Weibull para mi conjunto de datos.

He intentado esto:

import scipy.stats as s import numpy as np import matplotlib.pyplot as plt def weib(x,n,a): return (a / n) * (x / n)**(a - 1) * np.exp(-(x / n)**a) data = np.loadtxt("stack_data.csv") (loc, scale) = s.exponweib.fit_loc_scale(data, 1, 1) print loc, scale x = np.linspace(data.min(), data.max(), 1000) plt.plot(x, weib(x, loc, scale)) plt.hist(data, data.max(), normed=True) plt.show() 

Y consigue esto:

 (2.5827280639441961, 3.4955032285727947) 

Y una distribución que se parece a esto:

Distribución de Weibull utilizando Scipy.

    He estado usando el exponweib después de leer este http://www.johndcook.com/distributions_scipy.html . También he probado las otras funciones de Weibull en scipy (¡por si acaso!).

    En Matlab (usando la herramienta de ajuste de distribución – vea la captura de pantalla) y en R (usando tanto la función de biblioteca MASS fitdistr como el paquete GAMLSS) obtengo los parámetros a (loc) y b (escala) más como 1.58463497 5.93030013. Creo que los tres métodos utilizan el método de máxima verosimilitud para el ajuste de la distribución.

    Distribución de Weibull utilizando Matlab.

    ¡He publicado mis datos aquí si te gustaría tener una oportunidad! Y para completar estoy usando Python 2.7.5, Scipy 0.12.0, R 2.15.2 y Matlab 2012b.

    ¿Por qué estoy obteniendo un resultado diferente?

    Mi conjetura es que desea estimar el parámetro de forma y la escala de la distribución de Weibull mientras mantiene la ubicación fija. La fijación de loc supone que los valores de sus datos y de la distribución son positivos con un límite inferior en cero.

    floc=0 mantiene la ubicación fija en cero, f0=1 mantiene el primer parámetro de forma de la weibull exponencial fija en uno.

     >>> stats.exponweib.fit(data, floc=0, f0=1) [1, 1.8553346917584836, 0, 6.8820748596850905] >>> stats.weibull_min.fit(data, floc=0) [1.8553346917584836, 0, 6.8820748596850549] 

    El ajuste comparado con el histogtwig se ve bien, pero no es muy bueno. Las estimaciones de los parámetros son un poco más altas que las que usted menciona que son de R y matlab.

    Actualizar

    Lo más cercano que puedo llegar a la ttwig que ahora está disponible es con un ajuste no restringido, pero usando valores iniciales. La ttwig es aún menos pico. Tenga en cuenta que los valores de ajuste que no tienen una f al frente se utilizan como valores de inicio.

     >>> from scipy import stats >>> import matplotlib.pyplot as plt >>> plt.plot(data, stats.exponweib.pdf(data, *stats.exponweib.fit(data, 1, 1, scale=02, loc=0))) >>> _ = plt.hist(data, bins=np.linspace(0, 16, 33), normed=True, alpha=0.5); >>> plt.show() 

    ajuste exponweib

    Es fácil verificar qué resultado es el verdadero MLE, solo necesita una función simple para calcular la probabilidad de registro:

     >>> def wb2LL(p, x): #log-likelihood return sum(log(stats.weibull_min.pdf(x, p[1], 0., p[0]))) >>> adata=loadtxt('/home/user/stack_data.csv') >>> wb2LL(array([6.8820748596850905, 1.8553346917584836]), adata) -8290.1227946678173 >>> wb2LL(array([5.93030013, 1.57463497]), adata) -8410.3327470347667 

    El resultado del método de fit de exponweib y R fitdistr (@Warren) es mejor y tiene mayor probabilidad de registro. Es más probable que sea el verdadero MLE. No es sorprendente que el resultado de GAMLSS sea diferente. Es un modelo estadístico completamente diferente: modelo aditivo generalizado.

    ¿Todavía no está convencido? Podemos dibujar un gráfico de límite de confianza 2D alrededor de MLE, ver el libro de Meeker y Escobar para más detalles. Región de confianza multidimensional

    De nuevo, esto verifica que la array([6.8820748596850905, 1.8553346917584836]) es la respuesta correcta ya que la probabilidad de logl es menor que en cualquier otro punto en el espacio de parámetros. Nota:

     >>> log(array([6.8820748596850905, 1.8553346917584836])) array([ 1.92892018, 0.61806511]) 

    Es posible que BTW1, ajuste MLE no se ajuste bien al histogtwig de distribución. Una manera fácil de pensar acerca de MLE es que MLE es el parámetro más probable dado los datos observados. No es necesario que el histogtwig se ajuste bien a la vista, eso será algo que minimice el error cuadrático medio.

    BTW2, sus datos parecen ser leptokurtic y sesgados a la izquierda, lo que significa que la distribución de Weibull puede no ajustarse bien a sus datos. Pruebe, por ejemplo, Gompertz-Logistic, que mejora la probabilidad de registro en otros 100. introduzca la descripción de la imagen aquíintroduzca la descripción de la imagen aquí ¡Aclamaciones!

    Tenía curiosidad acerca de su pregunta y, a pesar de que no es una respuesta, compara el resultado de Matlab con su resultado y con el resultado utilizando leastsq , que mostró la mejor correlación con los datos dados:

    introduzca la descripción de la imagen aquí

    El código es el siguiente:

     import scipy.stats as s import numpy as np import matplotlib.pyplot as plt import numpy.random as mtrand from scipy.integrate import quad from scipy.optimize import leastsq ## my distribution (Inverse Normal with shape parameter mu=1.0) def weib(x,n,a): return (a / n) * (x / n)**(a-1) * np.exp(-(x/n)**a) def residuals(p,x,y): integral = quad( weib, 0, 16, args=(p[0],p[1]) )[0] penalization = abs(1.-integral)*100000 return y - weib(x, p[0],p[1]) + penalization # data = np.loadtxt("stack_data.csv") x = np.linspace(data.min(), data.max(), 100) n, bins, patches = plt.hist(data,bins=x, normed=True) binsm = (bins[1:]+bins[:-1])/2 popt, pcov = leastsq(func=residuals, x0=(1.,1.), args=(binsm,n)) loc, scale = 1.58463497, 5.93030013 plt.plot(binsm,n) plt.plot(x, weib(x, loc, scale), label='weib matlab, loc=%1.3f, scale=%1.3f' % (loc, scale), lw=4.) loc, scale = s.exponweib.fit_loc_scale(data, 1, 1) plt.plot(x, weib(x, loc, scale), label='weib stack, loc=%1.3f, scale=%1.3f' % (loc, scale), lw=4.) plt.plot(x, weib(x,*popt), label='weib leastsq, loc=%1.3f, scale=%1.3f' % tuple(popt), lw=4.) plt.legend(loc='upper right') plt.show() 

    Sé que es una publicación antigua, pero acabo de enfrentar un problema similar y este hilo me ayudó a resolverlo. Pensé que mi solución podría ser útil para otros como yo:

     # Fit Weibull function, some explanation below params = stats.exponweib.fit(data, floc=0, f0=1) shape = params[1] scale = params[3] print 'shape:',shape print 'scale:',scale #### Plotting # Histogram first values,bins,hist = plt.hist(data,bins=51,range=(0,25),normed=True) center = (bins[:-1] + bins[1:]) / 2. # Using all params and the stats function plt.plot(center,stats.exponweib.pdf(center,*params),lw=4,label='scipy') # Using my own Weibull function as a check def weibull(u,shape,scale): '''Weibull distribution for wind speed u with shape parameter k and scale parameter A''' return (shape / scale) * (u / scale)**(shape-1) * np.exp(-(u/scale)**shape) plt.plot(center,weibull(center,shape,scale),label='Wind analysis',lw=2) plt.legend() 

    Alguna información extra que me ayudó a entender:

    La función de Scipy Weibull puede tomar cuatro parámetros de entrada: (a, c), loc y escala. Desea corregir la ubicación y el primer parámetro de forma (a), esto se hace con floc = 0, f0 = 1. El ajuste luego le dará el parámetro c y la escala, donde c corresponde al parámetro de forma de la distribución de Weibull de dos parámetros (a menudo utilizado en el análisis de datos de viento) y la escala corresponde a su factor de escala.

    De los documentos:

     exponweib.pdf(x, a, c) = a * c * (1-exp(-x**c))**(a-1) * exp(-x**c)*x**(c-1) 

    Si a es 1, entonces

     exponweib.pdf(x, a, c) = c * (1-exp(-x**c))**(0) * exp(-x**c)*x**(c-1) = c * (1) * exp(-x**c)*x**(c-1) = c * x **(c-1) * exp(-x**c) 

    A partir de esto, la relación con la función de Weibull del ‘análisis del viento’ debería ser más clara.

    Tuve el mismo problema, pero encontré que establecer loc=0 en exponweib.fit preparó la bomba para la optimización. Eso era todo lo que se necesitaba de la respuesta de @ user333700. No pude cargar tus datos: tu enlace de datos apunta a una imagen, no a datos. Así que hice una prueba en mis datos en su lugar:

    Parcela de distribución ajustada a datos problemáticos (bimodales?)

     import scipy.stats as ss import matplotlib.pyplot as plt import numpy as np N=30 counts, bins = np.histogram(x, bins=N) bin_width = bins[1]-bins[0] total_count = float(sum(counts)) f, ax = plt.subplots(1, 1) f.suptitle(query_uri) ax.bar(bins[:-1]+bin_width/2., counts, align='center', width=.85*bin_width) ax.grid('on') def fit_pdf(x, name='lognorm', color='r'): dist = getattr(ss, name) # params = shape, loc, scale # dist = ss.gamma # 3 params params = dist.fit(x, loc=0) # 1-day lag minimum for shipping y = dist.pdf(bins, *params)*total_count*bin_width sqerror_sum = np.log(sum(ci*(yi - ci)**2. for (ci, yi) in zip(counts, y))) ax.plot(bins, y, color, lw=3, alpha=0.6, label='%s err=%3.2f' % (name, sqerror_sum)) return y colors = ['r-', 'g-', 'r:', 'g:'] for name, color in zip(['exponweib', 't', 'gamma'], colors): # 'lognorm', 'erlang', 'chi2', 'weibull_min', y = fit_pdf(x, name=name, color=color) ax.legend(loc='best', frameon=False) plt.show() 

    Ha habido algunas respuestas a esto ya aquí y en otros lugares. Me gusta en la distribución de Weibull y los datos en la misma figura (con numpy y scipy)

    Todavía me tomó un tiempo dar con un ejemplo de juguete limpio, así que pensé que sería útil publicar.

     from scipy import stats import matplotlib.pyplot as plt #input for pseudo data N = 10000 Kappa_in = 1.8 Lambda_in = 10 a_in = 1 loc_in = 0 #Generate data from given input data = stats.exponweib.rvs(a=a_in,c=Kappa_in, loc=loc_in, scale=Lambda_in, size = N) #The a and loc are fixed in the fit since it is standard to assume they are known a_out, Kappa_out, loc_out, Lambda_out = stats.exponweib.fit(data, f0=a_in,floc=loc_in) #Plot bins = range(51) fig = plt.figure() ax = fig.add_subplot(1, 1, 1) ax.plot(bins, stats.exponweib.pdf(bins, a=a_out,c=Kappa_out,loc=loc_out,scale = Lambda_out)) ax.hist(data, bins = bins , normed=True, alpha=0.5) ax.annotate("Shape: $k = %.2f$ \n Scale: $\lambda = %.2f$"%(Kappa_out,Lambda_out), xy=(0.7, 0.85), xycoords=ax.transAxes) plt.show() 

    el orden de ubicación y escala está desordenado en el código:

     plt.plot(x, weib(x, scale, loc)) 

    el parámetro de escala debe venir primero.

    En la función de ajuste, hay 3 parámetros a considerar:

    1. Los parámetros de forma: en este caso, tenemos dos parámetros de forma, que pueden fijarse de acuerdo con f0 y f1. (Pruébalo por ti mismo!). Generalmente, el nombre del parámetro se denota por f% d, donde d es el número de forma.

    2. El parámetro de ubicación: usa floc para arreglar esto. Si no arregla floc, la media de los datos se genera como loc.

    3. El parámetro de escala: Use fscale para arreglar esto.

    El retorno de cualquier ajuste sale en este orden.

    Siguiendo las líneas de @ Peter9192, encontré la mejor opción para un CDF Weibull de ~ 20-30 muestras de datos utilizando lo siguiente: _,gamma,_alpha=scipy.stats.exponweib.fit(data,floc=0,f0=1)

    La fórmula para CDF es:

    1-np.exp(-np.power(x/alpha,gamma)) Los valores para los datos que estimé usando un método de estimador de KM, correspondientes a una distribución de Weibull me dieron buenos valores.

    Para corregir a como 1, no encontré loc = 0, scale = 1 como el mejor método como se puede ver claramente en los 4 valores de parámetros devueltos. En segundo lugar, usando gamma, alfa de ella no dio la media Weibull correcta.

    Por último, confirmé qué método funciona mejor calculando la media de la distribución de Weibull usando: Mean=alpha*scipy.special.gamma(1+(1/gamma)) Los valores que obtuve correspondían a mi aplicación.

    Puede consultar la media y las fórmulas de CDF aquí para referencia: https://en.m.wikipedia.org/wiki/Weibull_distribution