Scipy: ajuste lognormal

Ha habido bastantes publicaciones sobre el manejo de la distribución de lognorm con Scipy, pero todavía no entiendo el truco.

El parámetro 2 lognormal normalmente se describe con los parámetros \mu y \sigma que corresponden a Scipys loc=0 y \sigma=shape , \mu=np.log(scale) .

En scipy, distribución lognormal – parámetros , podemos leer cómo generar una lognorm(\mu,\sigma) usando el exponencial de una distribución aleatoria. Ahora intentemos algo más:

UNA)

¿Cuál es el problema en la creación de un lognorm directamente:

 # lognorm(mu=10,sigma=3) # so shape=3, loc=0, scale=np.exp(10) ? x=np.linspace(0.01,20,200) sample_dist = sp.stats.lognorm.pdf(x, 3, loc=0, scale=np.exp(10)) shape, loc, scale = sp.stats.lognorm.fit(sample_dist, floc=0) print shape, loc, scale print np.log(scale), shape # mu and sigma # last line: -7.63285693379 0.140259699945 # not 10 and 3 

SEGUNDO)

Uso los valores de retorno de un ajuste para crear una distribución ajustada. Pero otra vez estoy haciendo algo mal al parecer:

 samp=sp.stats.lognorm(0.5,loc=0,scale=1).rvs(size=2000) # sample param=sp.stats.lognorm.fit(samp) # fit the sample data print param # does not coincide with shape, loc, scale above! x=np.linspace(0,4,100) pdf_fitted = sp.stats.lognorm.pdf(x, param[0], loc=param[1], scale=param[2]) # fitted distribution pdf = sp.stats.lognorm.pdf(x, 0.5, loc=0, scale=1) # original distribution plt.plot(x,pdf_fitted,'r-',x,pdf,'g-') plt.hist(samp,bins=30,normed=True,alpha=.3) 

lognorm

Hice las mismas observaciones: un ajuste libre de todos los parámetros falla la mayor parte del tiempo. Puede ayudar proporcionando una mejor estimación inicial, no es necesario arreglar el parámetro.

 samp = stats.lognorm(0.5,loc=0,scale=1).rvs(size=2000) # this is where the fit gets it initial guess from print stats.lognorm._fitstart(samp) (1.0, 0.66628696413404565, 0.28031095750445462) print stats.lognorm.fit(samp) # note that the fit failed completely as the parameters did not change at all (1.0, 0.66628696413404565, 0.28031095750445462) # fit again with a better initial guess for loc print stats.lognorm.fit(samp, loc=0) (0.50146296628099118, 0.0011019321419653122, 0.99361128537912125) 

También puede crear su propia función para calcular la estimación inicial, por ejemplo:

 def your_func(sample): # do some magic here return guess stats.lognorm._fitstart = your_func 

Me di cuenta de mis errores:

A) Las muestras que estoy dibujando deben provenir del método .rvs . Como así: sample_dist = sp.stats.lognorm.rvs(3, loc=0, scale=np.exp(10), size=2000)

B) El ajuste tiene algunos problemas. Cuando arreglamos el parámetro loc , el ajuste es mucho mejor. param=sp.stats.lognorm.fit(samp, floc=0)

Este problema se ha solucionado en nuevas versiones scipy. Después de actualizar scipy0.9 a scipy0.14, el problema desaparece.

Respondí aqui

Dejo el código aquí también para perezosos: D

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

El truco es entender estas dos cosas:

  1. Si la EXP de una variable es NORMAL con MU y STD -> EXP (X) ~ scipy.stats.lognorm (s = sigma, loc = 0, scale = np.exp (mu))
  2. Si su variable (x) TIENE LA FORMA DE LOGNORMAL, el modelo será scipy.stats.lognorm (s = sigmaX, loc = 0, scale = muX) con:
    • muX = np.mean (np.log (x))
    • sigmaX = np.std (np.log (x))