python - Ajuste de una distribución de Weibull usando Scipy

numpy distribution

Question

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:

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.

¡He publicado mis datos here 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?

Answer 1

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, utilizando 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 se correspondían con mi aplicación.

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

Answer 2

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.

Nuevamente, 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 histograma de distribución. Una manera fácil de pensar en MLE es que MLE es el parámetro más probable dado los datos observados. No es necesario que el histograma 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 adaptarse bien a sus datos. Pruebe, por ejemplo, Gompertz-Logistic, que mejora la probabilidad de registro en otros 100. ¡Aclamaciones!

Answer 3

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 encontrar 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()

Answer 4

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 Weibull de dos parámetros (a menudo se usa en el análisis de datos de viento) y la escala corresponde a su factor de escala.

De 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.

Answer 5

Supongo 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 histograma se ve bien, pero no 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 trama que ahora está disponible es con un ajuste no restringido, pero utilizando valores iniciales. La trama 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()

Answer 6

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:

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()

Answer 7

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 answer de @ user333700. No pude cargar tus datos: tu here apunta a una imagen, no a datos. Así que hice una prueba en mis datos en su lugar:

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()

Answer 8

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.