Ajuste de la distribución lognormal usando Scipy vs Matlab

Estoy tratando de ajustar una distribución lognormal usando Scipy. Ya lo he hecho con Matlab antes, pero debido a la necesidad de extender la aplicación más allá del análisis estadístico, estoy intentando reproducir los valores ajustados en Scipy.

A continuación se muestra el código de Matlab que solía ajustar a mis datos:

% Read input data (one value per line) x = []; fid = fopen(file_path, 'r'); % reading is default action for fopen disp('Reading network degree data...'); if fid == -1 disp('[ERROR] Unable to open data file.') else while ~feof(fid) [x] = [x fscanf(fid, '%f', [1])]; end c = fclose(fid); if c == 0 disp('File closed successfully.'); else disp('[ERROR] There was a problem with closing the file.'); end end [f,xx] = ecdf(x); y = 1-f; parmhat = lognfit(x); % MLE estimate mu = parmhat(1); sigma = parmhat(2); 

Y aquí está la ttwig equipada:

introduzca la descripción de la imagen aquí

Ahora aquí está mi código de Python con el objective de lograr lo mismo:

 import math from scipy import stats from statsmodels.distributions.empirical_distribution import ECDF # The same input is read as a list in Python ecdf_func = ECDF(degrees) x = ecdf_func.x ccdf = 1-ecdf_func.y # Fit data shape, loc, scale = stats.lognorm.fit(degrees, floc=0) # Parameters sigma = shape # standard deviation mu = math.log(scale) # meanlog of the distribution fit_ccdf = stats.lognorm.sf(x, [sigma], floc=1, scale=scale) 

Aquí está el ajuste usando el código de Python.

introduzca la descripción de la imagen aquí

Como puede ver, ambos conjuntos de códigos son capaces de producir buenos ajustes, al menos visualmente hablando.

El problema es que hay una gran diferencia en los parámetros estimados mu y sigma.

De Matlab: mu = 1.62 sigma = 1.29. Desde Python: mu = 2.78 sigma = 1.74.

¿Por qué hay tal diferencia?

Nota: he comprobado que ambos conjuntos de datos ajustados son exactamente iguales. Misma cantidad de puntos, misma distribución.

¡Su ayuda es muy apreciada! Gracias por adelantado.

Otra información:

 import scipy import numpy import statsmodels scipy.__version__ '0.9.0' numpy.__version__ '1.6.1' statsmodels.__version__ '0.5.0.dev-1bbd4ca' 

La versión de Matlab es R2011b.

Edición:

Como se demuestra en la respuesta a continuación, la falla está en Scipy 0.9. Puedo reproducir los resultados mu y sigma de Matlab usando Scipy 11.0.

Una manera fácil de actualizar tu Scipy es:

 pip install --upgrade Scipy 

Si no tienes pip (¡deberías!):

 sudo apt-get install pip 

Hay un error en el método de fit en scipy 0.9.0 que se ha corregido en versiones posteriores de scipy.

La salida del script a continuación debe ser:

 Explicit formula: mu = 4.99203450, sig = 0.81691086 Fit log(x) to norm: mu = 4.99203450, sig = 0.81691086 Fit x to lognorm: mu = 4.99203468, sig = 0.81691081 

pero con scipy 0.9.0, es

 Explicit formula: mu = 4.99203450, sig = 0.81691086 Fit log(x) to norm: mu = 4.99203450, sig = 0.81691086 Fit x to lognorm: mu = 4.23197270, sig = 1.11581240 

La siguiente secuencia de comandos de prueba muestra tres formas de obtener los mismos resultados:

 import numpy as np from scipy import stats def lognfit(x, ddof=0): x = np.asarray(x) logx = np.log(x) mu = logx.mean() sig = logx.std(ddof=ddof) return mu, sig # A simple data set for easy reproducibility x = np.array([50., 50, 100, 200, 200, 300, 500]) # Explicit formula my_mu, my_sig = lognfit(x) # Fit a normal distribution to log(x) norm_mu, norm_sig = stats.norm.fit(np.log(x)) # Fit the lognormal distribution lognorm_sig, _, lognorm_expmu = stats.lognorm.fit(x, floc=0) print "Explicit formula: mu = %10.8f, sig = %10.8f" % (my_mu, my_sig) print "Fit log(x) to norm: mu = %10.8f, sig = %10.8f" % (norm_mu, norm_sig) print "Fit x to lognorm: mu = %10.8f, sig = %10.8f" % (np.log(lognorm_expmu), lognorm_sig) 

Con la opción ddof=1 en el estándar. dev. Cálculo para utilizar la estimación de varianza imparcial:

 In [104]: x Out[104]: array([ 50., 50., 100., 200., 200., 300., 500.]) In [105]: lognfit(x, ddof=1) Out[105]: (4.9920345004312647, 0.88236457185021866) 

Hay una nota en la documentación lognfit de matlab que dice que cuando no se utiliza la censura, lognfit calcula sigma utilizando la raíz cuadrada del estimador imparcial de la varianza. Esto corresponde a usar ddof = 1 en el código anterior.