Est-ce que rstan ou mon approximation de grille est incorrecte: décider entre des estimations quantiles contradictoires dans l'inférence bayésienne

8

J'ai un modèle pour obtenir des estimations bayésiennes de la taille de la population Net probabilité de détection dans une distribution binomiale uniquement basée sur le nombre observé d’objets observés : pour . Pour simplifier, nous supposons que N est fixé à la même valeur inconnue pour chaque y_i . Dans cet exemple, y = 53,57,66,67,73 .θy

p(N,θ|y)Bin(y|N,θ)N
{N|NZNmax(y)}×(0,1)Nyiy=53,57,66,67,73

Ce modèle, lorsqu'il est estimé en rstan, diverge des résultats obtenus à partir d'une approximation en grille de la partie postérieure. J'essaie de comprendre pourquoi. (Les lecteurs intéressés pourraient trouver que cette question fait suite à ma réponse ici .)

rstan Approximation

Pour référence, il s'agit du code 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))

Notez que je lance thetaen 2 simplex. C'est juste pour la simplicité. La quantité d'intérêt est theta[1]; theta[2]est évidemment une information superflue.

De plus, N est une valeur réelle ( rstanaccepte uniquement les paramètres à valeur réelle car c'est une méthode de gradient), j'ai donc écrit une distribution binomiale à valeur réelle.

Résultats 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

Approximation de la grille

L'approximation de la grille a été produite comme ci-dessous. Les contraintes de mémoire m'empêchent de faire une grille plus fine sur mon ordinateur portable.

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)

J'ai utilisé l'approximation de la grille pour produire cet affichage de la densité postérieure. On voit que le postérieur est en forme de banane; ce type de postérieur peut être problématique pour la métrique euclidienne HMC. (La gravité de la forme de la banane est en fait supprimée ici car est sur l'échelle logarithmique.) Si vous pensez à la forme de la banane pendant une minute, vous vous rendrez compte qu'elle doit se trouver sur la ligne . (De plus, l'approximation de la grille affichée dans ce graphique n'est pas normalisée pour des raisons de clarté - sinon la banane est un peu trop étroite pour être clairement identifiée.)Ny¯=θN

postérieur sur une grille

Résultats d'approximation de la grille

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

Discussion

Le quantile de 97,5% pour est beaucoup plus grand dans mon modèle que pour l'approximation de la grille, mais ses quantiles sont similaires à l'approximation de la grille autrement. J'interprète cela comme indiquant que les deux méthodes sont généralement d'accord. Cependant, je ne sais pas comment interpréter l'écart dans le quantile de 97,5%.Nrstan

J'ai développé plusieurs explications possibles pour ce qui pourrait expliquer la divergence entre l'approximation de la grille et les résultats de l' rstanéchantillonnage HMC-NUTS, mais je ne sais pas comment comprendre si une, les deux ou aucune explication n'est correcte.

  1. Rstan a tort et la grille est correcte. La densité en forme de banane est problématique rstan, d'autant que dérive vers , donc ces quantités de queue ne sont pas fiables. Nous pouvons voir de la parcelle de la partie postérieure sur la grille que la queue est très forte à des valeurs plus grandes .N+N
  2. Rstan a raison et la grille est fausse. La grille fait deux approximations qui peuvent saper les résultats. Premièrement, la grille n'est qu'un ensemble fini de points sur un sous-espace postérieur, c'est donc une approximation approximative. En second lieu , parce que c'est un sous - espace fini, nous avoir faussement déclaré qu'il y ait 0 probabilité a posteriori sur les valeurs plus grand que notre plus grande valeur de grille pour . De même, est meilleur pour entrer dans les queues de la grille, donc ses quanitles de queue sont corrects.NNrstan

J'avais besoin de plus d'espace pour clarifier un point de la réponse de Juho. Si je comprends bien, nous pouvons intégrer hors de la partie postérieure pour obtenir la distribution bêta-binomiale: θ

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

Dans notre cas, et car nous avons un préalable uniforme sur . Je crois que le postérieur devrait alors être où car . Mais cela semble diverger énormément de la réponse de Juho. Où ai-je mal tourné?α=1β=1θp(N|y)N1i=1Kp(yi|N,α=1,β=1)K=#(y)p(N)=N1

Sycorax dit de réintégrer Monica
la source

Réponses:

4

Cliffs: Rstan semble être (plus proche) correct sur la base d'une approche qui intègre out analytiquement et évalue dans une grille assez grande.θP(N)P(yN)

Pour obtenir la partie postérieure de , il est en fait possible d'intégrer analytiquement: où est la longueur de . Maintenant, puisque a Beta avant (ici ) et Beta est conjugué au binôme, suit également une distribution Beta. Par conséquent, la distribution de est bêta-binomialeNθ

P(yN)=P(y1N)×P(y2N,y1)×P(y3N,y1,y2)×P(yKN,y1,,yK1)
KyθBeta(1,1)θN,y1,,ykyk+1N,y1,,ykpour laquelle il existe une expression sous forme fermée des probabilités en termes de fonction Gamma. Par conséquent, nous pouvons évaluer en calculant les paramètres pertinents des bêta-binômes et en multipliant les probabilités bêta-binomiales. Le code MATLAB suivant utilise cette approche pour calculer pour et se normalise pour obtenir le postérieur. 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

Le cdf à est de , donc je suppose que c'est suffisant, mais on pourrait vouloir étudier la sensibilité à l'augmentation du maximum . Le cdf à n'est que de , votre grille d'origine manque en effet une masse de probabilité de queue assez importante par rapport à l'objectif de trouver le quantile de . Les quantiles que j'obtiens sont N=1000000,9990maxN=500000NN=50000,97630,975

Quantile0,0250,250,50,750,975N951492354784750

Avertissement: je n'ai pas beaucoup testé le code, il peut y avoir des erreurs (et évidemment il pourrait y avoir des problèmes numériques avec cette approche aussi). Cependant, les quantiles obtenus sont assez proches des résultats de Rstan, donc je suis plutôt confiant.

Juho Kokkala
la source
(+1) Merci pour l'intérêt! Je prendrai votre signal et jouerai avec ces résultats dans R et je vous répondrai.
Sycorax dit Réintégrer Monica le
Pourriez-vous s'il vous plaît clarifier pourquoi la postérieure de la distribution bêta-binomiale est votre expression, plutôt que celle que j'ai ajoutée au bas de ma question? Il me semble que le postérieur devrait simplement être le produit de la vraisemblance bêta-binomiale et de l'a priori. Mais les résultats semblent tout à fait faux!
Sycorax dit Réintégrer Monica le
1
Il est important de mettre à jour les paramètres de la distribution bêta au fur et à mesure que les y sont traités, car tous les y partagent la même . L'équation à la fin de la question suppose un séparé pour chaque , qui est un modèle différent. θθy
Tom Minka