Eu tenho um modelo para obter estimativas bayesianas do tamanho da população e probabilidade de detecção em uma distribuição binomial baseada apenas no número observado de objetos observados : para . Por simplicidade, assumimos que N é fixado no mesmo valor desconhecido para cada y_i . Neste exemplo, y = 53,57,66,67,73 .
Este modelo, quando estimado em rstan
, diverge dos resultados obtidos a partir de uma aproximação da grade posterior. Estou tentando definir o porquê. (Os leitores interessados podem achar que esta pergunta é uma continuação da minha resposta aqui .)
rstan
Aproximação
Para referência, este é o código do rstan.
raftery.model <- "
data{
int I;
int y[I];
}
parameters{
real<lower=max(y)> N;
simplex[2] theta;
}
transformed parameters{
}
model{
vector[I] Pr_y;
for(i in 1:I){
Pr_y[i] <- binomial_coefficient_log(N, y[i])
+multiply_log(y[i], theta[1])
+multiply_log((N-y[i]), theta[2]);
}
increment_log_prob(sum(Pr_y));
increment_log_prob(-log(N));
}
"
raft.data <- list(y=c(53,57,66,67,72), I=5)
system.time(fit.test <- stan(model_code=raftery.model, data=raft.data,iter=10))
system.time(fit <- stan(fit=fit.test, data=raft.data,iter=10000,chains=5))
Observe que eu faço a conversão theta
como 2 simplex. Isso é apenas por simplicidade. A quantidade de interesse é theta[1]
; obviamente theta[2]
é informação supérflua.
Além disso, é um valor real ( rstan
só aceita parâmetros com valor real porque é um método de gradiente), então escrevi uma distribuição binomial com valor real.
Resultados Rstan
mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
N 1078.75 256.72 15159.79 94.44 148.28 230.61 461.63 4575.49 3487 1
theta[1] 0.29 0.00 0.19 0.01 0.14 0.27 0.42 0.67 2519 1
theta[2] 0.71 0.00 0.19 0.33 0.58 0.73 0.86 0.99 2519 1
lp__ -19.88 0.02 1.11 -22.89 -20.31 -19.54 -19.09 -18.82 3339 1
Aproximação da grade
A aproximação da grade foi produzida como abaixo. As restrições de memória me impedem de fazer uma grade mais fina no meu laptop.
theta <- seq(0+1e-10,1-1e-10, len=1e3)
N <- round(seq(72, 5000, len=1e3)); N[2]-N[1]
grid <- expand.grid(N,theta)
y <- c(53,57,66,67,72)
raftery.prob <- function(x, z=y){
N <- x[1]
theta <- x[2]
exp(sum(dbinom(z, size=N, prob=theta, log=T)))/N
}
post <- matrix(apply(grid, 1, raftery.prob), nrow=length(N), ncol=length(theta),byrow=F)
post.norm <- post/sum(post)
Eu usei a aproximação da grade para produzir essa exibição da densidade posterior. Podemos ver que o posterior é em forma de banana; esse tipo de posterior pode ser problemático para o HMC métrico euclidiano. (A severidade da forma da banana é realmente suprimida aqui, já que está na escala do log.) Se você pensar na forma da banana por um minuto, perceberá que ela deve estar na linha . (Além disso, a aproximação da grade exibida neste gráfico não é normalizada por razões de clareza - caso contrário, a banana é um pouco estreita demais para ser vista com clareza.)ˉ y = θ N
Resultados da aproximação de grade
do.call(cbind, lapply(c(0.025, .25, .5, .75, .975), function(quantile){
approx(y=N, x=cumsum(rowSums(post.norm))/sum(post.norm), xout=quantile)
}))
[,1] [,2] [,3] [,4] [,5]
x 0.025 0.25 0.5 0.75 0.975
y 92.55068 144.7091 226.7845 443.6359 2475.398
Discussão
O quantil de 97,5% para é muito maior no meu modelo do que na aproximação de grade, mas seus quantis são semelhantes à aproximação de grade. Eu interpreto isso como indicando que os dois métodos geralmente estão de acordo. Mas não sei como interpretar a discrepância no quantil de 97,5%.rstan
Eu desenvolvi várias explicações possíveis para o que pode estar explicando a divergência entre a aproximação da grade e os resultados da rstan
amostragem HMC-NUTS, mas não sei como entender se uma, ambas ou nenhuma explicação está correta.
- Rstan está errado e a grade está correta. A densidade em forma de banana é problemática
rstan
, especialmente quando se desloca para , portanto essas quantidades de cauda não são confiáveis. Podemos ver o enredo do posterior sobre a grade que a cauda é muito forte em valores maiores .+ ∞ N - Rstan está correto e a grade está errada. A grade faz duas aproximações que podem prejudicar os resultados. Primeiro, a grade é apenas um conjunto finito de pontos sobre um subespaço posterior, portanto é uma aproximação aproximada. Em segundo lugar, porque é um sub espaço finito, estamos falsamente declarando que haja 0 probabilidade posterior sobre os valores maior do que o nosso maior valor de grade para . Da mesma forma, é melhor entrar nas caudas da grade, de modo que seus quanitles de cauda estejam corretos.N
rstan
Eu precisava de mais espaço para esclarecer um ponto da resposta de Juho. Se bem entendi, podemos integrar partir do posterior para obter a distribuição beta-binomial:
No nosso caso, e porque temos um uniforme anterior em . Eu acredito que o posterior deve então ser onde porque . Mas isso parece divergir amplamente da resposta de Juho. Onde eu errei?β = 1 θ p ( N | y ) α N - 1 Π K i = 1 P ( y i | N , α = 1 , β = 1 ) K = # ( y )
fonte