Necessita de algoritmo para calcular a probabilidade relativa de que os dados sejam amostrados da distribuição normal versus distribuição normal

13

Digamos que você tenha um conjunto de valores e deseja saber se é mais provável que eles tenham sido amostrados de uma distribuição gaussiana (normal) ou amostrados de uma distribuição lognormal?

É claro que, idealmente, você saberia algo sobre a população ou sobre as fontes de erro experimental, para ter informações adicionais úteis para responder à pergunta. Mas aqui, suponha que só temos um conjunto de números e nenhuma outra informação. O que é mais provável: amostragem de uma distribuição gaussiana ou amostragem de uma distribuição lognormal? Quanto mais provável? O que eu estou esperando é um algoritmo para selecionar entre os dois modelos, e espero quantificar a probabilidade relativa de cada um.

Harvey Motulsky
fonte
1
Pode ser um exercício divertido tentar caracterizar a distribuição sobre distribuições na natureza / literatura publicada. Então, novamente - nunca será mais do que um exercício divertido. Para um tratamento sério, você pode procurar uma teoria que justifique sua escolha ou fornecer dados suficientes - visualize e teste a qualidade do ajuste de cada distribuição candidata.
21713 JohnRos
3
Se é uma questão de generalização a partir da experiência, eu diria que distribuições distorcidas positivamente são o tipo mais comum, especialmente para variáveis ​​de resposta que são de interesse central e que lognormals são mais comuns que normais. Um volume de 1962 O cientista especula editado pelo famoso estatístico IJ Good que incluiu uma peça anônima "Regras de trabalho de Bloggins", contendo a afirmação "A distribuição normal do log é mais normal que a normal". (Várias das outras regras são fortemente estatística.)
Nick Cox
Parece que interpreto sua pergunta de maneira diferente de JohnRos e ansoestevez. Para mim, sua pergunta soa como uma seleção simples de modelo , ou seja, uma questão de calcular , onde M é a distribuição normal ou log-normal e D é seus dados. Se a seleção de modelos não é o que você procura, pode esclarecer? P(MD)MD
Lucas
@ Lucas Acho que sua interpretação não é muito diferente da minha. Em ambos os casos, você precisa fazer suposições a priori .
anxoestevez
2
Por que não apenas calcular a razão de verossimilhança generalizada e alertar o usuário quando favorecer o log-normal?
Scortchi - Restabelece Monica

Respostas:

7

Você pode adivinhar o tipo de distribuição ajustando cada distribuição (normal ou normal de log) aos dados pela probabilidade máxima e comparando a probabilidade de log em cada modelo - o modelo com a maior probabilidade de log sendo o mais adequado. Por exemplo, em R:

# log likelihood of the data given the parameters (par) for 
# a normal or lognormal distribution
logl <- function(par, x, lognorm=F) {
    if(par[2]<0) { return(-Inf) }
    ifelse(lognorm,
    sum(dlnorm(x,par[1],par[2],log=T)),
    sum(dnorm(x,par[1],par[2],log=T))
    )
}

# estimate parameters of distribution of x by ML 
ml <- function(par, x, ...) {
    optim(par, logl, control=list(fnscale=-1), x=x, ...)
}

# best guess for distribution-type
# use mean,sd of x for starting parameters in ML fit of normal
# use mean,sd of log(x) for starting parameters in ML fit of lognormal
# return name of distribution type with highest log ML
best <- function(x) {
    logl_norm <- ml(c(mean(x), sd(x)), x)$value
        logl_lognorm <- ml(c(mean(log(x)), sd(log(x))), x, lognorm=T)$value
    c("Normal","Lognormal")[which.max(c(logl_norm, logl_lognorm))]
}

Agora gere números a partir de uma distribuição normal e ajuste uma distribuição normal por ML:

set.seed(1)
x = rnorm(100, 10, 2)
ml(c(10,2), x)

Produz:

$par
[1] 10.218083  1.787379

$value
[1] -199.9697
...

Compare a probabilidade de log para o ajuste de ML das distribuições normal e lognormal:

ml(c(10,2), x)$value # -199.9697
    ml(c(2,0.2), x, lognorm=T)$value # -203.1891
best(x) # Normal

Tente com uma distribuição lognormal:

best(rlnorm(100, 2.6, 0.2)) # lognormal

A atribuição não será perfeita, dependendo de n, média e sd:

> table(replicate(1000, best(rnorm(500, 10, 2))))

Lognormal    Normal 
        6       994 
> table(replicate(1000, best(rlnorm(500, 2.6, 0.2))))

Lognormal    Normal 
      999         1 
waferthin
fonte
1
Você não precisa encontrar as estimativas de parâmetro de probabilidade máxima numericamente para o normal ou para o log-normal (embora isso mostre como você generalizou a ideia para comparação de outras distribuições). Além disso, abordagem muito sensata.
Scortchi - Restabelecer Monica
Eu mal usei R ou o conceito de máxima probabilidade, então aqui está uma pergunta básica. Eu sei que não podemos comparar o AIC (ou BIC) de ajustar uma distribuição normal aos dados vs. aos logs dos dados, porque o AIC ou o BIC não seriam comparáveis. É necessário ajustar dois modelos a um conjunto de dados (sem transformações; sem exclusões externas, etc.), e a transformação dos dados alterará o AIC ou o BIC, independentemente de fazer a comparação falsa. E o ML? Essa comparação é legítima?
Harvey Motulsky
Encontramos as distribuições normal e lognormal mais adequadas aos dados e calculamos a probabilidade de observar os dados assumindo que eles eram dessas distribuições (a probabilidade ou p(X|\theta)). Não estamos transformando os dados. Imprimimos a distribuição cuja probabilidade de observação dos dados é mais alta. Essa abordagem é legítima, mas tem a desvantagem de não inferirmos a probabilidade do modelo dado os dados p(M|X), ou seja, a probabilidade de os dados serem de uma distribuição normal vs lognormal (por exemplo, p (normal) = 0,1, p (lognormal) = 0,9), diferentemente da abordagem bayesiana.
waferthin
1
@ Harvey É verdade, mas irrelevante - você perguntou sobre o ajuste de distribuições normais versus log-normais aos mesmos dados, e é isso que a whannymahoots está respondendo. Como o número de parâmetros livres é o mesmo para os dois modelos, comparar AICs ou BICs reduz a comparação de probabilidade de log.
Scortchi - Restabelecer Monica
@wannymahoots Qualquer prévia razoável para uma abordagem bayesiana nesse contexto - baseando-se em estimar as probabilidades relativas de que um usuário de software está tentando ajustar dados normais ou normais - será tão pouco informativa que dará resultados semelhantes a uma abordagem com base apenas na probabilidade.
Scortchi - Restabelecer Monica
11

M{Normal,Log-normal}X={x1,...,xN}

P(MX)P(XM)P(M).

A parte difícil é obter a probabilidade marginal ,

P(XM)=P(Xθ,M)P(θM)dθ.

p(θM)XY={logx1,...,logxNYX,

P(XM=Log-Normal)=P(YM=Normal)i|1xi|.

P(θM)P(σ2,μM=Normal)P(M)

Exemplo:

P(μ,σ2M=Normal)m0=0,v0=20,a0=1,b0=100

insira a descrição da imagem aqui

Segundo Murphy (2007) (Equação 203), a probabilidade marginal da distribuição normal é então dada por

P(XM=Normal)=|vN|12|v0|12b0a0bnaNΓ(aN)Γ(a0)1πN/22N

aN,bN,vNP(μ,σ2X,M=Normal)

vN=1/(v01+N),mN=(v01m0+ixi)/vN,aN=a0+N2,bN=b0+12(v01m02vN1mN2+ixi2).

Eu uso os mesmos hiperparâmetros para a distribuição log-normal,

P(XM=Log-normal)=P({registrox1,...,registroxN}M=Normal)Eu|1xEu|.

Para uma probabilidade anterior do log-normal de 0,1, P(M=Log-normal)=0,1e dados extraídos da seguinte distribuição log-normal,

enter image description here

o posterior se comporta assim:

enter image description here

A linha sólida mostra a probabilidade mediana posterior para diferentes desenhos de NOs pontos de dados. Observe que, para pouco ou nenhum dado, as crenças estão próximas das crenças anteriores. Para cerca de 250 pontos de dados, o algoritmo quase sempre tem certeza de que os dados foram extraídos de uma distribuição log-normal.

Ao implementar as equações, seria uma boa ideia trabalhar com densidades de log em vez de densidades. Mas, caso contrário, deve ser bem direto. Aqui está o código que eu usei para gerar os gráficos:

https://gist.github.com/lucastheis/6094631

Lucas
fonte
4

Parece que você está procurando algo bastante pragmático para ajudar analistas que provavelmente não são estatísticos profissionais e precisam de algo para levá-los a fazer o que deveriam ser técnicas exploratórias padrão, como analisar gráficos de qq, gráficos de densidade, etc.

Nesse caso, por que não fazer simplesmente um teste de normalidade (Shapiro-Wilk ou o que quer que seja) nos dados originais e um nos dados transformados em log, e se o segundo valor de p for maior, levante um sinalizador para o analista considerar usar uma transformação de log ? Como bônus, cuspa um gráfico 2 x 2 do gráfico da linha de densidade e do gráfico qqnorm dos dados brutos e transformados.

Tecnicamente, isso não responderá sua pergunta sobre a probabilidade relativa, mas me pergunto se é tudo o que você precisa.

Peter Ellis
fonte
Clever. Maybe this is enough, and avoids the need to explain likelihood calculations.... Thanks.
Harvey Motulsky