2012-01-05 29 views
22

Quiero ajustar la distribución lognormal a mis datos, usando python scipy.stats.lognormal.fit. De acuerdo con manual, fit devuelve forma, loc, escala parámetros. Pero la distribución lognormal normalmente solo necesita two parameters: media y desviación estándar.scipy, lognormal distribution - parameters

¿Cómo interpretar los resultados de la función scipy fit? Cómo obtener mean y std.dev.?

Respuesta

20

Las distribuciones en scipy se codifican de una manera genérica wrt dos ubicación de parámetros y escala de forma que la ubicación es el parámetro (loc) que desplaza la distribución a la izquierda o la derecha, mientras que scale es el parámetro que comprime o extiende la distribución .

Para los dos parámetro de distribución logarítmica normal, el "medio" y "std dev" corresponden a log (scale) y shape (se puede dejar que loc=0).

El siguiente ejemplo ilustra cómo encajar una distribución logarítmica normal para encontrar los dos parámetros de interés:

In [56]: import numpy as np 

In [57]: from scipy import stats 

In [58]: logsample = stats.norm.rvs(loc=10, scale=3, size=1000) # logsample ~ N(mu=10, sigma=3) 

In [59]: sample = np.exp(logsample) # sample ~ lognormal(10, 3) 

In [60]: shape, loc, scale = stats.lognorm.fit(sample, floc=0) # hold location to 0 while fitting 

In [61]: shape, loc, scale 
Out[61]: (2.9212650122639419, 0, 21318.029350592606) 

In [62]: np.log(scale), shape # mu, sigma 
Out[62]: (9.9673084420467362, 2.9212650122639419) 
+4

dos observaciones: hubo un error en scipy 0,9 con floc que se ha fijado en scipy 0,10 http://projects.scipy.org/scipy/ticket/1536 y segunda debido a la parametrización genérico la distribución lognormal no tiene una parametrización habitual, por ejemplo http://projects.scipy.org/scipy/ticket/1502. – user333700

+0

gracias por los parches. No obtuve el segundo: de los comentarios en el enlace, parece que no hay ningún error en realidad ... –

+0

@JakubM. Sí, si está utilizando el último scipy (0.10), la respuesta/ejemplo dado anteriormente no se contradice con ninguno de los tickets mencionados en el comentario por parte del usuario33700. – ars

5

Me acabo de pasar algún tiempo trabajando esto y quería documentar aquí: Si desea obtener la probabilidad la densidad (en el punto x) a partir de los tres valores de retorno de lognorm.fit (les permite llamar (shape, loc, scale)), es necesario utilizar esta fórmula:

x = 1/(shape*((x-loc)/scale)*sqrt(2*pi)) * exp(-1/2*(log((x-loc)/scale)/shape)**2)/scale 

Así como una ecuación que es (es locµ, shape es σ y scale es α):

x = \frac{1}{(x-\mu)\cdot\sqrt{2\pi\sigma^2}} \cdot e^{-\frac{log(\frac{x-\mu}{\alpha})^2}{2\sigma^2}}

+0

¿Por qué dividimos por (1/alfa) al final de su fórmula? – bioslime

+1

Esta no es mi fórmula, pero la manera en que scipy funciona. No veo ninguna división por (1/α} al final, así que asumiré que estás hablando de la división por 'escala' - corrígeme si no entendí. Miré este código hace un tiempo, pero si recuerde correctamente, eso es lo que se hace con cada tipo de distribución. Pero tenga en cuenta que hay otra 'escala' al principio y se cancelan mutuamente (como puede ver en la ecuación). – Chronial

-1

Creo que esto ayudará. Estaba buscando el mismo problema durante mucho tiempo y finalmente encontré una solución para mi problema. En mi caso, estaba tratando de ajustar algunos datos a la distribución lognormal usando el módulo scipy.stats.lognorm. Sin embargo, cuando finalmente obtuve los parámetros del modelo, no pude encontrar una manera de replicar mis resultados usando la media y el estándar de los datos.

En el siguiente código, explico a partir de los parámetros de media y estándar cómo producir una muestra de datos distribuidos normalmente usando el módulo scipy.stats.norm. Utilizando esos datos, ajusto el modelo normal (norm_dist_fitted) y también creo un modelo normal usando la media y la desviación estándar (mu, sigma) extraída de los datos.

El modelo original que produce los datos, ajustado y producido por (mu-sigma) -pair se compara en un gráfico.

Fig1


En la siguiente sección del código, I utilizar los datos normales para producir una muestra lognormal-distribuido. Para hacerlo, tenga en cuenta que las muestras lognormales serán exponenciales de la muestra original. Por lo tanto, la media y la desviación estándar de la muestra exponencial serán (exp(mu) y exp(sigma)).

I equipados los datos producidos a un lognormal (desde el registro de mi muestra (exp (x)) se distribuye normalmente y sigue los supuestos del modelo lognormal.

Para producir un modelo lognormal de la media y la desviación estándar de los datos originales (x) el código será:

lognorm_dist = scipy.stats.lognorm(s=sigma, loc=0, scale=np.exp(mu)) 

Sin embargo, si los datos ya está en el espacio exponencial (exp (x)), entonces usted tiene que utilizar:

muX = np.mean(np.log(x)) 
sigmaX = np.std(np.log(x)) 
scipy.stats.lognorm(s=sigmaX, loc=0, scale=muX) 

Fig2

import scipy 
import matplotlib.pyplot as plt 
import seaborn as sns 
import numpy as np 

mu = 10 # Mean of sample !!! Make sure your data is positive for the lognormal example 
sigma = 1.5 # Standard deviation of sample 
N = 2000 # Number of samples 

norm_dist = scipy.stats.norm(loc=mu, scale=sigma) # Create Random Process 
x = norm_dist.rvs(size=N) # Generate samples 

# Fit normal 
fitting_params = scipy.stats.norm.fit(x) 
norm_dist_fitted = scipy.stats.norm(*fitting_params) 
t = np.linspace(np.min(x), np.max(x), 100) 

# Plot normals 
f, ax = plt.subplots(1, sharex='col', figsize=(10, 5)) 
sns.distplot(x, ax=ax, norm_hist=True, kde=False, label='Data X~N(mu={0:.1f}, sigma={1:.1f})'.format(mu, sigma)) 
ax.plot(t, norm_dist_fitted.pdf(t), lw=2, color='r', 
     label='Fitted Model X~N(mu={0:.1f}, sigma={1:.1f})'.format(norm_dist_fitted.mean(), norm_dist_fitted.std())) 
ax.plot(t, norm_dist.pdf(t), lw=2, color='g', ls=':', 
     label='Original Model X~N(mu={0:.1f}, sigma={1:.1f})'.format(norm_dist.mean(), norm_dist.std())) 
ax.legend(loc='lower right') 
plt.show() 


# The lognormal model fits to a variable whose log is normal 
# We create our variable whose log is normal 'exponenciating' the previous variable 

x_exp = np.exp(x) 
mu_exp = np.exp(mu) 
sigma_exp = np.exp(sigma) 

fitting_params_lognormal = scipy.stats.lognorm.fit(x_exp, floc=0, scale=mu_exp) 
lognorm_dist_fitted = scipy.stats.lognorm(*fitting_params_lognormal) 
t = np.linspace(np.min(x_exp), np.max(x_exp), 100) 

# Here is the magic I was looking for a long long time 
lognorm_dist = scipy.stats.lognorm(s=sigma, loc=0, scale=np.exp(mu)) 

# The trick is to understand these two things: 
# 1. If the EXP of a variable is NORMAL with MU and STD -> EXP(X) ~ scipy.stats.lognorm(s=sigma, loc=0, scale=np.exp(mu)) 
# 2. If your variable (x) HAS THE FORM of a LOGNORMAL, the model will be scipy.stats.lognorm(s=sigmaX, loc=0, scale=muX) 
# with: 
# - muX = np.mean(np.log(x)) 
# - sigmaX = np.std(np.log(x)) 


# Plot lognormals 
f, ax = plt.subplots(1, sharex='col', figsize=(10, 5)) 
sns.distplot(x_exp, ax=ax, norm_hist=True, kde=False, 
      label='Data exp(X)~N(mu={0:.1f}, sigma={1:.1f})\n X~LogNorm(mu={0:.1f}, sigma={1:.1f})'.format(mu, sigma)) 
ax.plot(t, lognorm_dist_fitted.pdf(t), lw=2, color='r', 
     label='Fitted Model X~LogNorm(mu={0:.1f}, sigma={1:.1f})'.format(lognorm_dist_fitted.mean(), lognorm_dist_fitted.std())) 
ax.plot(t, lognorm_dist.pdf(t), lw=2, color='g', ls=':', 
     label='Original Model X~LogNorm(mu={0:.1f}, sigma={1:.1f})'.format(lognorm_dist.mean(), lognorm_dist.std())) 
ax.legend(loc='lower right') 
plt.show() 
+0

Proporcione un poco de una explicación para acompañar su respuesta. –

+0

@ CollinM.Barrett ¡Hecho !: DI espero que esté bien: D –