Prodotto di probabilità troppo piccolo - R dà solo 0

Question

Sto provando a calcolare questa distribuzione posteriore in R. Il problema è che il numeratore, che è il prodotto di un gruppo di dbern (p_i, y_i) < 1, è troppo piccolo. (Il mio n è circa 1500). Quindi, R sputa 0 e anche i valori posteriori per tutti \ theta sono 0.

Per chiarire, ogni y_i ha il proprio p_i, insieme questi p_i fanno un vettore di n elementi per n y's. Ogni theta ha il proprio vettore n-elemento di p_i.

Un esempio riproducibile (del numeratore)

p <- sample(seq(0.001,0.999,by=0.01), 1500, replace=T) 
y <- sample(c(0,1), 1500, replace=T) 
dbern(y, p) # 1500-element vector, each element < 1 
prod(dbern(y, p)) # produces 0 
exp(sum(log(dbern(y, p)))) # produces 0 

EDIT (contesto): Sto facendo un analisi bayesiana punto di cambio (jstor.org/stable/25791783 - occidentale e Kleykamp 2004). A differenza del continuo y nel foglio, il mio y è binario, quindi sto usando il metodo di aumento dei dati in Albert e Chib (1993). Con questo metodo, la probabilità di y è Bernoulli, con p = cdf-normal (x'B).

Quindi, come p dipende da theta? È perché theta è il punto di svolta. Una delle x è un tempo fittizio - se theta = 10, ad esempio, il tempo fittizio = 1 per tutte le osservazioni dopo il giorno 10 e = 0 per tutte le osservazioni prima del giorno 10.

Quindi, p dipende da x, x dipende da theta - quindi, p dipende da theta.

Ho bisogno della quantità di cui sopra perché è il pieno condizionale di theta nel campionamento Gibbs.

Answer 1

Ho fatto questa domanda sulla Croce Convalidato anche, e glen_b mi ha dato questa soluzione (testato) qui di seguito:

---

Questo è un problema comune con il calcolo delle probabilità per tutti i tipi di modelli; il tipo di cose che sono comunemente fatte sono lavorare sui log e usare un fattore di ridimensionamento comune che porta i valori in un intervallo più ragionevole.

In questo caso, suggerirei:

Fase 1: Scegli un abbastanza "tipica" θ, θ0. Dividere la formula per il numeratore e il denominatore del termine generale dal numeratore per θ = θ0, in modo da ottenere qualcosa che sarà molto meno probabilità di underflow.

Passo 2: lavorare sulla scala di registro, questo significa che il numeratore è un exp di somme di differenze di log e il denominatore è una somma di exp di somme di differenze di log.

NB: Se uno dei tuoi p è 0 o 1, estrarli separatamente e non prendere i registri di tali termini; sono facili da valutare così come sono!

I termini usuali nel numeratore tenderanno ad essere di dimensioni più moderate, e quindi in molte situazioni il numeratore e il denominatore sono entrambi ragionevoli.

Se al denominatore è presente una gamma di dimensioni, aggiungere quelle più piccole prima di aggiungere quelle più grandi.

Se uno o pochi termini dominano pesantemente, si dovrebbe concentrare l'attenzione sul rendere il calcolo per quelli relativamente accurati.

Answer 2

Un modo per affrontare problemi di precisione come questo è lavorare nello spazio di log. Comunque questo introduce un log-sum-product dal denominatore, che generalmente può essere doloroso.

Se si sta calcolando il posteriore ai fini dell'ottimizzazione, tenere presente che potrebbe essere possibile eliminare completamente il denominatore: non è necessario normalizzarlo per trovare uno argmax.

Answer 3

Beh mi sono imbattuto il tuo esempio, si sta (come previsto) ottenendo 0 perché c'è una 0

p <- sample(seq(0,1,by=0.01), 1500, replace=T) 
y <- sample(c(0,1), replace=T) 
x <- dbern(y, p) 
any(x == 0) 
## [1] TRUE