Ajuste de la distribución logarítmica normal en R vs. SciPy

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 usado antes. Usando Scipy, pude obtener una forma y escala de 1 y 3.1626716539637488e + 90 - números muy diferentes. También intenté usar la exp del meanlog y sdlog, pero sigo obteniendo gráficos extraños.

He leído todos los documentos que puedo sobre scipy y todavía estoy confundido sobre lo que significan los parámetros de forma y escala en este caso. ¿Tendría sentido codificar la función yo mismo? Sin embargo, parece propenso a errores, ya que soy nuevo en scipy.

SCIPY Lognormal (AZUL) vs. R Lognormal (ROJO):

¿Alguna idea sobre qué dirección tomar? Los datos se ajustan muy bien con el modelo R, por cierto, por lo que si se ve como algo más en Python, siéntase libre de compartir.

¡Gracias!

Actualizar:

Estoy ejecutando Scipy 0.11

Aquí hay un subconjunto de los datos. La muestra real es 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 (probé varios en caso de que confundiera 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?

Actualizar:

volver a ejecutar con la sugerencia de "UPQuark":

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

Sin embargo, la forma del gráfico es muy similar, con el pico alrededor de 21.

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

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

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

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

loc: no es equivalente, esto se resta de sus datos para que 0 se convierta en el límite mínimo del rango de datos.

scale - , donde es la media del registro de la variante. (Al realizar el ajuste, normalmente usaría la media muestral del registro de los datos).$\exp{\mu}$$\mu$

forma: la desviación estándar del registro de la variante.

Llamé a lognorm.pdf(x, 0.55, 0, numpy.exp(4.29))donde están los argumentos (x, shape, loc, scale) respectivamente, y generé 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 tu curva R.

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

Si bien este enfoque general unifica muy bien todas las distribuciones, para lognormal puede crear cierta confusión debido a la forma en que otros paquetes definen los parámetros. Aún así, es muy simple hacer coincidir cualquier distribución lognormal si quiere decir log (la media de la distribución subyacente) y sdlog (la desviación estándar de la distribución subyacente).

Primero, asegúrese de establecer el parámetro de ubicación en 0. Luego, establezca el parámetro de forma en el valor de sdlog. Finalmente, establezca el parámetro de escala en math.exp (meanlog). Por lo tanto, rv = scipy.stats.lognorm (0.5511349, scale = math.exp (4.2991610)) creará un objeto de distribución cuyo pdf coincida exactamente con su curva generada en R. Como x = numpy.linspace (0,180,1000); plot (x, rv.pdf (x)) verificará.

Básicamente, la distribución lognormal SciPy es una generalización de la distribución lognormal estándar que coincide exactamente con el estándar cuando se establece el parámetro de ubicación en 0.

Al ajustar datos con el método .fit, también puede usar palabras clave, f0..fn, floc y fshape para mantener fijos cualquiera de los parámetros de forma, ubicación y / o escala y solo se ajusta a las otras variables. Para la distribución lognormal, esto es muy útil ya que generalmente sabe que el parámetro de ubicación debe fijarse en 0. Por lo tanto, scipy.stats.lognorm.fit (conjunto de datos, floc = 0) siempre devolverá el parámetro de ubicación como 0 y solo variará el otro Parámetros de forma y escala.

El ajuste lognormal puntiagudo devuelve forma, ubicación y 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 estimaciones razonables 1.0, 0.09, 0.86, y cuando lo traza, debe 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 lo normal.

Espero que esto ayude.

Parece que la distribución en Scipy para 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 esto: 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 usar una función de densidad lognormal en Python. Si a alguien le gustaría agregar a esto, no dude en hacerlo.

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