Ajustando a distribuição log-normal em R vs. SciPy

10

Eu ajustei um modelo lognormal usando R com um conjunto de dados. Os parâmetros resultantes foram:

meanlog = 4.2991610 
sdlog = 0.5511349

Gostaria de transferir esse modelo para o Scipy, que nunca usei antes. Usando o Scipy, consegui obter uma forma e escala de 1 e 3,1626716539637488e + 90 - números muito diferentes. Eu também tentei usar o exp do meanlog e sdlog, mas continuo obtendo um gráfico bizarro.

Eu li todos os documentos que posso no scipy e ainda estou confuso sobre o significado dos parâmetros de forma e escala neste caso. Faria sentido codificar a função eu mesmo? Isso parece propenso a erros, pois sou novo no scipy.

SCIPY Lognormal (AZUL) vs. R Lognormal (VERMELHO): Scipy Lognormal (AZUL) vs. R Lognormal (VERMELHO)

Alguma idéia de qual direção tomar? A propósito, os dados se encaixam muito bem com o modelo R; portanto, se parecer com algo mais em Python, fique à vontade para compartilhar.

Obrigado!

Atualizar:

Estou executando o Scipy 0.11

Aqui está um subconjunto dos dados. A amostra real é 38k +, com uma média 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:

Estou trabalhando em uma função para capturar o 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

No entanto, isso me deu os números a seguir (tentei vários no caso de estar entendendo o significado de sdlog e 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

Alguma ideia?

Atualizar:

reexecutando com a sugestão "UPQuark":

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

A forma do gráfico é muito semelhante, no entanto, com o pico acontecendo em torno de 21.

Lillian Milagros Carrasquillo
fonte
Esta pergunta e resposta podem ajudar: stackoverflow.com/questions/8747761/…
jbowman
Obrigado, eu achei isso e aprendi o "encaixe" com o lognormal. No entanto, minhas perguntas são: por que eu receberia distribuições tão diferentes?
Lillian Milagros Carrasquillo
Você está usando o SciPy 0.9? Além disso, você poderia postar seus dados ou um subconjunto deles?
jbowman
Atualizada! É Scipy 0,11 por sinal. Assim, os erros que eu li sobre não deve ser relevante;)
Lillian Milagros Carrasquillo

Respostas:

11

Eu lutei através do código-fonte, para chegar à seguinte interpretação da rotina lognormal scipy.

x-locescalaLognormal(σ)

onde é o parâmetro "shape". σ

A equivalência entre os parâmetros scipy e o parâmetro R é a seguinte:

loc - Sem equivalente, isso é subtraído dos seus dados para que 0 se torne o menor do intervalo dos dados.

scale - , onde é a média do log da variável. (Ao ajustar, normalmente você usaria a média da amostra do log dos dados.)expμμ

shape - o desvio padrão do log da variável.

Eu chamei lognorm.pdf(x, 0.55, 0, numpy.exp(4.29))onde estão os argumentos (x, shape, loc, scale) respectivamente e gerei os seguintes valores:

x pdf

10 0.000106

20 0,002275

30 0,006552

40 0,009979

50 0,114557

60 0,1113479

70 0,103327

80 0,008941

90 0,007494

100 0,006155

que parecem combinar muito bem com sua curva R.

jbowman
fonte
Obrigado, @JBowman, essa é exatamente a explicação que eu precisava e a saída é precisamente a minha distribuição.
Lillian Milagros Carrasquillo
8

A distribuição lognormal no SciPy se encaixa na estrutura geral de todas as distribuições no SciPy. Todos eles têm uma palavra-chave de escala e localização (o padrão é 0 e 1 se não for fornecido explicitamente). Isso permite que todas as distribuições sejam alteradas e escaladas a partir de suas especificações normalizadas, com implicações claras nas estatísticas da distribuição. As distribuições normalmente também têm um ou mais parâmetros de "forma" (embora alguns, como a distribuição normal, não precisem de parâmetros adicionais).

Embora essa abordagem geral unifique bem todas as distribuições, para lognormal, ela pode criar alguma confusão devido à maneira como outros pacotes definem os parâmetros. Ainda assim, é muito simples corresponder a qualquer distribuição normal de log se você quiser log (a média da distribuição subjacente) e sdlog (o desvio padrão da distribuição subjacente).

Primeiro, verifique se você definiu o parâmetro location como 0. Em seguida, defina o parâmetro shape para o valor de sdlog. Por fim, defina o parâmetro scale como math.exp (meanlog). Portanto, rv = scipy.stats.lognorm (0.5511349, scale = math.exp (4.2991610)) criará um objeto de distribuição cujo pdf corresponda exatamente à sua curva gerada por R exatamente. Como x = numpy.linspace (0,180,1000); plot (x, rv.pdf (x)) será verificado.

Basicamente, a distribuição lognormal SciPy é uma generalização da distribuição lognormal padrão que corresponde exatamente ao padrão ao definir o parâmetro de localização como 0.

Ao ajustar os dados com o método .fit, você também pode usar as palavras-chave, f0..fn, floc e fshape para manter fixos os parâmetros de forma, localização e / ou escala e ajustar-se apenas às outras variáveis. Para a distribuição lognormal, isso é muito útil, pois geralmente você sabe que o parâmetro location deve ser fixado em 0. Portanto, scipy.stats.lognorm.fit (conjunto de dados, floc = 0) sempre retornará o parâmetro location como 0 e varia apenas o outro parâmetros de forma e escala.

Travis Oliphant
fonte
3

O ajuste normal do log da Scipy retorna forma, localização e escala. Acabei de executar o seguinte em uma matriz de dados de preços de amostra:

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

Isso me dá estimativas razoáveis ​​de 1,0, 0,09, 0,86 e, quando você o traça, deve levar em consideração todos os três parâmetros.

O parâmetro shape é o desvio padrão da distribuição normal subjacente e a escala é o exponencial da média do normal.

Espero que isto ajude.

upquark
fonte
Obrigado por responder! Depois de ter esses valores (loc, scale, shape), tento encontrar o pdf (x) para cada x que me interessa (aqui são valores de 0 a 180, exclusivos). scipy.stats.lognorm.pdf (i, loc, scale, shape) No entanto, plotando esses, recebo o gráfico acima.
Lillian Milagros Carrasquillo
OK, vi você mencionando apenas forma e escala, por isso mencionei que há três parâmetros retornados por padrão de fit (). Você também disse que está confuso sobre o significado dos parâmetros de forma e escala, e tentei abordar isso. Eu nunca tive o ajuste lognormal retornar valores absurdos como no seu caso, qual é o parâmetro location?
upquark
Apenas atualizei a pergunta para responder a isso. Obrigado por pensar sobre isso.
Lillian Milagros Carrasquillo
Chame scipy.stats.lognorm.pdf (x, shape, loc, scale) em vez de scipy.stats.lognorm.pdf (i, loc, scale, shape).
25712 upquark
Obrigado, upquark, eu fiz isso também com resultados semelhantes. Toda a forma do gráfico continua sendo muito diferente dos resultados esperados, resultando em R. Parece uma distribuição totalmente diferente daquela em R, na verdade.
Lillian Milagros Carrasquillo
1

Parece que a distribuição no Scipy para o lognormal não é a mesma que em R ou, geralmente, não é a mesma que a distribuição com a qual estou familiarizado. John D Cook tocou nisto: http://www.johndcook.com/blog/2010/02/03/statistical-distributions-in-scipy/ http://www.johndcook.com/distributions_scipy.html

No entanto, não encontrei nada conclusivo sobre como usar uma função de densidade lognormal no Python. Se alguém quiser adicionar isso, sinta-se à vontade.

Minha solução até agora é usar o pdf lognormal avaliado de 0 a 180 (exclusivo) e usado como um dicionário no script python.

Lillian Milagros Carrasquillo
fonte