Ajuste de la distribución log-normal en R frente a SciPy

Question

He ajustado un modelo lognormal usando R con un conjunto de datos. Los parámetros resultantes fueron:

meanlog = 4.2991610 
sdlog = 0.5511349

Me gustaría transferir este modelo a Scipy, que nunca he utilizado antes. Usando Scipy, pude obtener una forma y escala de 1 y 3.1626716539637488e+90 -- números muy diferentes. También he tratado de usar la exp de la meanlog y sdlog pero sigo obteniendo una gráfica bizarra.

He leído toda la documentación que he podido sobre scipy y sigo confundido sobre el significado de los parámetros shape y scale en este caso. ¿Tendría sentido codificar la función yo mismo? Eso parece propenso a errores, sin embargo, como soy nuevo en scipy.

SCIPY Lognormal (AZUL) frente a R Lognormal (ROJO):

¿Alguna idea sobre qué dirección tomar? Los datos se ajustan muy bien con el modelo de R, por cierto, así que si se parece a otra cosa en Python, siéntase libre de compartir.

Gracias.

Actualización:

Estoy ejecutando Scipy 0.11

Aquí está un subconjunto de los datos. La muestra real es de 38k+, con una media de 81,53627:

Subconjunto:

x
[60, 170, 137, 138, 81, 140, 78, 46, 1, 168, 138, 148, 145, 35, 82, 126, 66, 147, 88, 106, 80, 54, 83, 13, 102, 54, 134, 34]
numpy.mean(x)
99.071428571428569

Alternativamente:

Estoy trabajando en una función para capturar el pdf:

def lognoral(x, mu, sigma):
    a = 1 / (x * sigma * numpy.sqrt(2 * numpy.pi) )
    b = - (numpy.log(x) - mu) ^ 2 / (2 * sigma ^ 2)
    p = a * numpy.exp(b)
    return p

Sin embargo, esto me da los siguientes números (he probado varios por si me confundía el significado de sdlog y meanlog):

>>> lognormal(54,4.2991610, 0.5511349)
0.6994656085799437
 >>> lognormal(54,numpy.exp(4.2991610), 0.5511349)
0.9846125119455129
>>> lognormal(54,numpy.exp(4.2991610), numpy.exp(0.5511349))
0.9302407837304372

¿Alguna idea?

Actualización:

volver a correr con la sugerencia de "UPQuark":

forma, loc, escala (1.0, 50.03445923295007, 19.074457156766517)

Sin embargo, la forma del gráfico es muy similar, ya que el pico se produce en torno a los 21 años.

Answer 1

Me abrí camino a través del código fuente, para llegar a la siguiente interpretación de la rutina scipy lognormal.

$\frac{x-\text{loc}}{\text{scale}} \sim \text{Lognormal}(\sigma)$

donde $\sigma$ es el parámetro "forma".

La equivalencia entre los parámetros de scipy y los de R es la siguiente:

loc - No es equivalente, se resta de los datos para que el 0 se convierta en el mínimo del rango de los datos.

escala - $\exp{\mu}$ , donde $\mu$ es la media del logaritmo de la variable. (Al realizar el ajuste, normalmente se utiliza la media muestral del logaritmo de los datos).

shape - la desviación estándar del logaritmo de la variable.

He llamado lognorm.pdf(x, 0.55, 0, numpy.exp(4.29)) donde los argumentos son (x, shape, loc, scale) respectivamente, y generan los siguientes valores:

x pdf

10 0.000106

20 0.002275

30 0.006552

40 0.009979

50 0.114557

60 0.113479

70 0.103327

80 0.008941

90 0.007494

100 0.006155

que parecen coincidir bastante bien con su curva R.

Answer 2

La distribución lognormal en SciPy se ajusta al marco general de todo en SciPy. Todas tienen una palabra clave de escala y ubicación (que por defecto son 0 y 1 si no se proporcionan explícitamente). Esto permite que todas las distribuciones sean desplazadas y escaladas desde su especificación normalizada con claras implicaciones para las estadísticas de la distribución. Las distribuciones suelen tener también uno o más parámetros de "forma" (aunque algunas, como la distribución normal, no necesitan ningún parámetro adicional).

Aunque este planteamiento general unifica muy bien todo las distribuciones, para lognormal puede crear cierta confusión debido a la forma en que otros paquetes definen los parámetros. Aun así, es muy sencillo igualar cualquier distribución lognormal si se meanlog (la media de la distribución subyacente) y sdlog (la desviación estándar de la distribución subyacente).

En primer lugar, asegúrese de establecer el parámetro de ubicación en 0. A continuación, establezca el parámetro de forma en el valor de sdlog. Por último, establezca el parámetro de escala a math.exp(meanlog). Así, rv = scipy.stats.lognorm(0.5511349, scale=math.exp(4.2991610)) creará un objeto de distribución cuyo pdf coincide exactamente con la curva generada por R. Como x = numpy.linspace(0,180,1000); plot(x, rv.pdf(x)) lo verificará.

Básicamente, la distribución lognormal de SciPy es una generalización de la distribución lognormal estándar que coincide exactamente con la estándar al establecer el parámetro de localización en 0.

Cuando se ajustan los datos con el método .fit, también se pueden utilizar las palabras clave, f0..fn, floc y fshape para mantener fijos cualquiera de los parámetros de forma, localización y/o escala y sólo ajustar sobre las otras variables. Para la distribución lognormal esto es muy útil ya que normalmente se sabe que el parámetro de localización debe ser fijo a 0. Así, scipy.stats.lognorm.fit(dataset, floc=0) siempre devolverá el parámetro de localización como 0 y sólo variará los otros parámetros de forma y escala.

Answer 3

El ajuste lognormal de Scipy devuelve la forma, la ubicación y la escala. Acabo de ejecutar lo siguiente en una matriz de datos de precios de muestra:

shape, loc, scale = st.lognorm.fit(d_in["price"])

Esto me da unas estimaciones razonables de 1,0, 0,09, 0,86, y al trazarlo hay que tener en cuenta los tres parámetros.

El parámetro de forma es la desviación estándar de la distribución normal subyacente, y la escala es la exponencial de la media de la normal.

Espero que esto ayude.

Answer 4

Parece que la distribución en Scipy para la lognormal no es la misma que en R, o en general, no es la misma que la distribución con la que estoy familiarizado. John D Cook ha tocado este tema: http://www.johndcook.com/blog/2010/02/03/statistical-distributions-in-scipy/ http://www.johndcook.com/distributions\_scipy.html

Sin embargo, no he encontrado nada concluyente sobre cómo utilizar una función de densidad lognormal en Python. Si alguien quiere añadir algo a esto, por favor, siéntase libre.

Mi solución hasta ahora es utilizar el pdf lognormal evaluado de 0 a 180 (exclusivo), y utilizado como diccionario en el python script.