python - probabilidad - Scipy, distribución lognormal-parámetros
poisson en python (3)
Creo que esto ayudará. Estuve 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 utilizando la media y la información estándar de y.
En el código a continuación, explico a partir de los parámetros media y estándar cómo producir una muestra de datos distribuida normalmente utilizando el módulo scipy.stats.norm. Usando esos datos, ajusté el modelo normal ( norm_dist_fitted
) y también norm_dist_fitted
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 par (mu-sigma) se compara en una gráfica.
En la siguiente sección del código, uso los datos normales para producir una muestra distribuida lognormal. Para hacerlo, observe que las muestras lognormales serán la exponencial 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)
).
Ajusté los datos producidos a un lognormal
(ya que el registro de mi muestra (exp (x)) se distribuye normalmente y sigue los supuestos del modelo lognormal.
Para producir un modelo lognormal a partir de la media y la desviación estándar de sus datos originales (x), el código será:
lognorm_dist = scipy.stats.lognorm(s=sigma, loc=0, scale=np.exp(mu))
Sin embargo, si sus datos ya están en el espacio exponencial (exp (x)), entonces tiene que usar:
muX = np.mean(np.log(x))
sigmaX = np.std(np.log(x))
scipy.stats.lognorm(s=sigmaX, loc=0, scale=muX)
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()
Quiero ajustar la distribución lognormal a mis datos, usando python scipy.stats.lognormal.fit
. De acuerdo con el manual , el fit
devuelve los parámetros de forma, ubicación y escala . Pero, la distribución lognormal normalmente necesita solo dos parámetros : media y desviación estándar.
¿Cómo interpretar los resultados de la función de fit
scipy? ¿Cómo obtener media y std.dev.?
Las distribuciones en scipy están codificadas de forma genérica con la ubicación y escala de dos parámetros, de modo que la ubicación es el parámetro ( loc
) que desplaza la distribución hacia la izquierda o hacia la derecha, mientras que la scale
es el parámetro que comprime o estira la distribución.
Para los dos parámetros de distribución lognormal, "mean" y "std dev" corresponden a log ( scale
) y shape
(puede dejar loc=0
).
Lo siguiente ilustra cómo ajustar una distribución lognormal 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)
Solo dediqué un tiempo a resolver esto y quise documentarlo aquí: si desea obtener la densidad de probabilidad (en el punto x
) de los tres valores de retorno de lognorm.fit
(llamémoslos (shape, loc, scale)
), Necesitas usar esta fórmula:
x = 1 / (shape*((x-loc)/scale)*sqrt(2*pi)) * exp(-1/2*(log((x-loc)/scale)/shape)**2) / scale
Entonces, como una ecuación que es ( loc
es µ
, la shape
es σ
y la scale
es α
):