Calcul de la probabilité lorsque

8

J'essaie de calculer cette distribution postérieure:

(θ|)=i=1npiyi(1pi)1yiallθ,pi|θi=1npiyi(1pi)1yi

Le problème est que le numérateur, qui est le produit d'un tas de Bernoulli(pi,yi)les probabilités sont trop faibles. (Man est grande, environ 1500).

Par conséquent, les valeurs postérieures pour tous θ tous sont calculés à 0 (je fais des calculs dans R).

Pour clarifier, chaque yi a sa propre pi, ensemble ces pifait un vecteur de n éléments pour n y's. Chaqueθ a sa propre n-élément vecteur de pje.

EDIT: Ajout d'un exemple de reproduction (pour le numérateur)

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
la source
Avez-vous plutôt essayé de calculer la somme des journaux?
Ansari
1
Il y a une discussion connexe ici . Il contient des informations supplémentaires sur certains détails de ces calculs.
Glen_b -Reinstate Monica

Réponses:

7

Il s'agit d'un problème courant avec le calcul des probabilités pour toutes sortes de modèles; les tâches généralement effectuées consistent à travailler sur les journaux et à utiliser un facteur d'échelle commun qui ramène les valeurs dans une plage plus raisonnable.

Dans ce cas, je suggère:

Étape 1: Choisissez un assez "typique" θ, θ0. Divisez la formule pour le numérateur et le dénominateur du terme général par le numérateur pourθ=θ0, afin d'obtenir quelque chose qui sera beaucoup moins susceptible de déborder.

Étape 2: travailler sur l'échelle logarithmique, cela signifie que le numérateur est une exp de sommes de différences de journaux, et le dénominateur est une somme d'exp de sommes de différences de journaux.

NB: Si l'un de vos p vaut 0 ou 1, retirez-les séparément et ne prenez pas les journaux de ces termes; ils sont faciles à évaluer tels quels!

[En termes plus généraux, cette mise à l'échelle et travail sur l'échelle du journal peut être considérée comme prenant un ensemble de log-vraisemblances, lje et ce faisant: Journal(jeelje)=c+Journal(jeelje-c). Un choix évident pourc est de faire le plus grand terme 0, ce qui nous laisse: Journal(jeelje)=maxje(lje)+Journal(jeelje-maxje(lje)). Notez que lorsque vous avez un numérateur et un dénominateur, vous pouvez utiliser le mêmecpour les deux, qui seront alors annulés. Dans ce qui précède, cela correspond à prendre le avec la probabilité de log la plus élevée.]θ0

Les termes usuels du numérateur auront tendance à être de taille plus modérée, et donc dans de nombreuses situations, le numérateur et le dénominateur sont tous deux relativement raisonnables.

S'il y a une gamme de tailles dans le dénominateur, additionnez les plus petites avant d'ajouter les plus grandes.

Si seuls quelques termes dominent fortement, vous devez concentrer votre attention sur la précision du calcul de ceux-ci.

Glen_b -Reinstate Monica
la source
Mais pour tout thêta, le numérateur passe toujours à 0. Comment puis-je diviser le terme général par le numérateur? (Étape 1)
Heisenberg
1
L'étape 1 est l' algèbre et non le calcul informatique. Son but est de vous donner quelque chose à l'étape 2 pour calculer qui ne déborde pas. À moins que vous ne disiez que c'est toujours algébriquement nul ... auquel cas vous faites sans doute quelque chose que vous ne devriez pas faire.
Glen_b -Reinstate Monica
ok - je vais essayer. Le numérateur n'est pas exactement 0, seulement très petit que R ne peut pas calculer. Merci!
Heisenberg
3
Cher Dieu, tu as raison! Merci beaucoup. Tout le monde dit "use log.likelihood" mais vous seul voyez vraiment le problème.
Heisenberg
1

Essayez de tirer parti des propriétés de l'utilisation des logarithmes et de la somme plutôt que de prendre le produit de nombres décimaux. Après la sommation, utilisez simplement l'anti-log pour le remettre dans votre forme la plus naturelle. Je pense que quelque chose comme ça devrait faire l'affaire

eXp(jen(yjelog(pje)+(1-yje)log(1-pje)))geXp(jenyjelog(pje)+(1-yje)log(1-pje))

philchalmers
la source
Le numérateur dans votre suggestion produit toujours un 0 car la somme dans exp () est toujours très négative (<-1000). Suis-je en train de faire quelque chose de mal? Merci de votre aide!
Heisenberg
Eh bien, si une valeur de p est en fait 0 ou 1, alors automatiquement son journal produira -inf et donc log (1-p). Sinon, je pense que les chiffres deviennent trop petits pour être ramenés à la forme originale.
philchalmers
2
Notez que vous pouvez ajouter et soustraire n'importe quelle constante c des termes à l'intérieur du exp()l'expression ci-dessus sans changer le résultat. réglagec égal au négatif de la valeur maximale de Journal(p(θ|-))fournit la meilleure précision numérique
probabilitéislogique