Cálculo da probabilidade quando

8

Estou tentando calcular esta distribuição posterior:

(θ|-)=Eu=1npEuyEu(1-pEu)1-yEutodosθ,pEu|θEu=1npEuyEu(1-pEu)1-yEu

O problema é que o numerador, que é o produto de um monte de probabilidades é muito pequeno. (Meu é grande, cerca de 1500).nBernoulli(pEu,yEu)n

Portanto, os valores posteriores para todos são calculados como 0 (estou fazendo cálculos em R).θ

Para esclarecer, cada tem seu próprio , juntos esses formam um vetor de elementos para . Cada tem seu próprio vetor de elemento de .yEupEupEunn yθnpEu

EDIT: Adicionando um exemplo de reprodução (para o numerador)

p <- sample(seq(0,1,by=0.01), 1500, replace=T)
y <- sample(c(0,1), 1500, replace=T)
dbern(y, p) # 1500-element vector, each element is < 1
prod(dbern(y, p)) # produce 0
exp(sum(log(dbern(y, p)))) # produce 0 since the sum is very negative
Heisenberg
fonte
Você tentou calcular a soma dos logs?
Ansari
1
Há discussões relacionadas aqui . Tem alguma discussão adicional sobre alguns dos detalhes de tais cálculos.
Glen_b -instala Monica

Respostas:

7

Esse é um problema comum com o cálculo de probabilidades para todos os tipos de modelos; os tipos de coisas comumente feitas são trabalhar nos logs e usar um fator de escala comum que traga os valores para um intervalo mais razoável.

Nesse caso, sugiro:

Etapa 1: Escolha um , θ 0 razoavelmente "típico" . Divida a fórmula para o numerador e o denominador do termo geral pelo numerador para θ = θ 0 , a fim de obter algo que será muito menos provável que ocorra.θθ0 0θ=θ0 0

Etapa 2: trabalhe na escala de log, isso significa que o numerador é uma exp de somas de diferenças de logs e o denominador é uma soma de exp de somas de diferenças de logs.

NB: Se algum dos seus p's for 0 ou 1, retire-os separadamente e não faça registros desses termos; eles são fáceis de avaliar como estão!

[Em termos mais gerais, este aumento de escala e trabalhando-a-a-registo de escala-pode ser visto como tendo um conjunto de log-probabilidades, e fazendo isso: log ( Σ i e L i ) = c + log ( Σ i e l i - c ) . Uma escolha óbvia para ceuEuregistro(EueeuEu)=c+registro(EueeuEu-c)cregistro(EueeuEu)=maxEu(euEu)+registro(EueeuEu-maxEu(euEu))cpara ambos, que será cancelado. Acima, isso corresponde a obter com a maior probabilidade de log.]θ0 0

Os termos usuais no numerador tendem a ter um tamanho mais moderado e, em muitas situações, o numerador e o denominador são relativamente razoáveis.

Se houver uma variedade de tamanhos no denominador, some os menores antes de adicionar os maiores.

Se apenas alguns termos dominam fortemente, concentre sua atenção em fazer o cálculo para aqueles relativamente precisos.

Glen_b -Reinstate Monica
fonte
Mas para todos os tetas, o numerador sempre vai para 0. Como eu divido o termo geral pelo numerador? (Passo 1)
Heisenberg
1
Etapa 1 é álgebra, não cálculo de computador. Seu objetivo é fornecer a você algo na Etapa 2 para calcular que não é insuficiente. A menos que você esteja dizendo que é sempre algebricamente zero ... nesse caso, você está sem dúvida fazendo algo que não deveria.
Glen_b -Reinstate Monica
ok - vou tentar. O numerador não é exatamente 0, apenas muito pequeno que R não pode computar. Obrigado!
precisa
3
Querido Deus, você está correto! Muito obrigada. Todo mundo fica dizendo "use log.likelihood", mas somente você realmente vê o problema.
precisa
1

Tente capitalizar as propriedades do uso dos logaritmos e da soma em vez de usar o produto de números decimais. Após o somatório, use o anti-log para colocá-lo novamente em sua forma mais natural. Eu acho que algo assim deve fazer o truque

exp(Eun(yEueuog(pEu)+(1-yEu)euog(1-pEu)))gexp(EunyEueuog(pEu)+(1-yEu)euog(1-pEu))

filósofos
fonte
O numerador na sua sugestão ainda produz um 0, pois a soma dentro de exp () ainda é muito negativa (<-1000). Estou fazendo algo errado? Obrigado pela ajuda!
precisa
Bem, se qualquer valor em p é realmente 0 ou 1, automaticamente o log dele produzirá -inf e o log (1-p). Caso contrário, acho que os números estão ficando muito pequenos para serem elevados à forma original.
philchalmers
2
cexp()cregistro(p(θ|-))