A aproximação de rstan ou minha grade está incorreta: decidindo entre estimativas quantílicas conflitantes na inferência bayesiana

8

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 .Nθy

p(N,θ|y)Bin(y|N,θ)N
{N|NZNmax(y)}×(0,1)Nyiy=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 thetacomo 2 simplex. Isso é apenas por simplicidade. A quantidade de interesse é theta[1]; obviamente theta[2]é informação supérflua.

Além disso, N é um valor real ( rstansó 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 = θ NNy¯=θN

posterior sobre uma grade

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%.Nrstan

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 rstanamostragem HMC-NUTS, mas não sei como entender se uma, ambas ou nenhuma explicação está correta.

  1. 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 .+ NN+N
  2. 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.NNNrstan

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: θ

p(y|N,α,β)=(Ny)Beta(y+α,Ny+β)Beta(α,β)

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 )α=1β=1θp(N|y)N1i=1Kp(yi|N,α=1,β=1)K=#(y)p(N)=N1

Sycorax diz restabelecer Monica
fonte

Respostas:

4

Penhascos: Rstan parece estar (mais próximo) correto com base em uma abordagem que integra fora analiticamente e avalia em uma grade bastante grande.θP(N)P(yN)

Para obter o posterior de , é realmente possível integrar analiticamente: que é o comprimento de . Agora, como tem Beta anterior (aqui ) e Beta é conjugado com binomial, também segue uma distribuição Beta. Portanto, a distribuição de é beta-binomialNθ

P(yN)=P(y1N)×P(y2N,y1)×P(y3N,y1,y2)×P(yKN,y1,,yK1)
KyθBeta(1,1)θN,y1,,ykyk+1N,y1,,ykpara o qual existe uma expressão fechada das probabilidades em termos da função Gamma. Portanto, podemos avaliar calculando os parâmetros relevantes dos binômios beta e multiplicando as probabilidades beta-binomiais. O código MATLAB a seguir usa essa abordagem para calcular para e normaliza para obter a posterior. P(yN)P(N)P(y|N)N=72,,500000
%The data
y = [53 57 66 67 72];

%Initialize
maxN = 500000;
logp = zeros(1,maxN); %log prior + log likelihood
logp(1:71) = -inf; 

for N = 72:maxN
    %Prior
    logp(N) = -log(N);

    %y1 has uniform distribution
    logp(N) = logp(N) - log(N+1);    
    a = 1;
    b = 1;

    %Rest of the measurements 
    for j = 2:length(y);
        %Update beta parameters
        a = a + y(j-1);
        b = b + N - y(j-1);

        %Log predictive probability of y_j (see Wikipedia article)
        logp(N) = logp(N) + gammaln(N+1) - gammaln(y(j) + 1) - ... 
         gammaln(N - y(j) + 1) + gammaln(y(j) + a) +  ...
         gammaln(N - y(j) + b) - gammaln(N + a + b) ...
            + gammaln(a+b) - gammaln(a) - gammaln(b);

    end
end

%Get the posterior of N
pmf = exp(logp - max(logp));
pmf = pmf/sum(pmf);
cdf = cumsum(pmf);

%Evaluate quantiles of interest
disp(cdf(5000))  %0.9763
for percentile = [0.025 0.25 0.5 0.75 0.975]
    disp(find(cdf>=percentile,1,'first'))
end

O cdf em é , então acho que é suficiente, mas é possível investigar a sensibilidade para aumentar o máximo . Como o cdf em é de apenas , sua grade original perde uma massa de probabilidade de cauda bastante significativa em comparação com o objetivo de encontrar o quantil . Os quantis que recebo são 0,9990N=1000000.9990maxN=500000N = 5000 0,9763 0,975 Quantil 0,025 0,25 0,5 0,75 0,975 N 95 149 235 235 478 4750NN=50000.97630.975

Quantile0.0250.250.50.750.975N951492354784750

Isenção de responsabilidade: não testei muito o código, pode haver erros (e, obviamente, também pode haver problemas numéricos com essa abordagem). No entanto, os quantis obtidos são muito próximos dos resultados de Rstan, por isso estou bastante confiante.

Juho Kokkala
fonte
(+1) Obrigado pelo interesse! Vou pegar sua sugestão e brincar com esses resultados em R e voltar para você.
Sycorax diz Restabelecer Monica
Você poderia esclarecer por que a parte posterior da distribuição beta-binomial é a sua expressão, e não a que adicionei na parte inferior da minha pergunta? Parece-me que o posterior deve ser simplesmente o produto da probabilidade beta-binomial e do anterior. Mas os resultados parecem estar bastante errados!
Sycorax diz Restabelecer Monica
1
É importante atualizar os parâmetros da distribuição Beta à medida que os y's são processados, porque todos os y compartilham o mesmo . A equação no final da pergunta assume um separado para cada , que é um modelo diferente. θ yθθy
21414 Tommy Minka