Tengo un modelo para lograr estimaciones bayesianas del tamaño de la población y la probabilidad de detección en una distribución binomial basada únicamente en el número observado de objetos observados : para . Por simplicidad, asumimos que N está fijado en el mismo valor desconocido para cada y_i . En este ejemplo, y = 53,57,66,67,73 .$N$$\theta$$y$

$$p(N,\theta|y)\propto \frac{ \text{Bin}(y|N,\theta)}{N}$$ $\left\{N|N\in\mathbb{Z}\land N\ge \max(y)\right\}\times(0,1)$$N$$y_i$$y=53,57,66,67,73$

Este modelo, cuando se estima en rstan, difiere de los resultados obtenidos de una aproximación de cuadrícula de la parte posterior. Estoy tratando de precisar por qué. (Los lectores interesados ​​pueden encontrar que esta pregunta es una continuación de mi respuesta aquí ).

rstan Aproximación

Como referencia, este es el código 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))

Tenga en cuenta que lanzo thetacomo un 2-simplex. Esto es solo por simplicidad. La cantidad de interés es theta[1]; Obviamente theta[2]es información superflua.

Además, $N$ es un valor real ( rstansolo acepta parámetros con valores reales porque es un método de gradiente), así que escribí una distribución binomial con valores reales.

Resultados de 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

Aproximación de cuadrícula

La aproximación a la cuadrícula se produjo como a continuación. Las restricciones de memoria me impiden hacer una cuadrícula más fina en mi computadora portátil.

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)

Utilicé la aproximación de cuadrícula para producir esta visualización de la densidad posterior. Podemos ver que la parte posterior tiene forma de plátano; Este tipo de posterior puede ser problemático para la HMC métrica euclidiana. (La gravedad de la forma del plátano en realidad se suprime aquí ya que está en la escala logarítmica). Si piensa en la forma del plátano por un minuto, se dará cuenta de que debe estar en la línea . (Además, la aproximación a la cuadrícula que se muestra en este gráfico no está normalizada por razones de claridad; de lo contrario, el plátano es un poco demasiado estrecho para distinguirlo claramente).$N$$\bar{y}=\theta N$

Resultados de aproximación de cuadrícula

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

Discusión

El cuantil del 97.5% para es mucho más grande en mi modelo que para la aproximación a la cuadrícula, de lo contrario, sus cuantiles son similares a la aproximación a la cuadrícula. Interpreto que esto indica que los dos métodos generalmente están de acuerdo. Sin embargo, no sé cómo interpretar la discrepancia en el cuantil del 97.5%.$N$rstan

He desarrollado varias explicaciones posibles de lo que podría explicar la divergencia entre la aproximación a la cuadrícula y los resultados del rstanmuestreo HMC-NUTS, pero no estoy seguro de cómo entender si una, ambas o ninguna de las explicaciones es correcta.

1. Rstan está equivocado y la cuadrícula es correcta. La densidad en forma de plátano es problemática rstan, especialmente cuando desvía hacia , por lo que estas cantidades de cola no son confiables. Podemos ver en la trama de la parte posterior sobre la cuadrícula que la cola es muy agudo en los valores más grandes .$N$$+\infty$$N$
2. Rstan es correcto y la cuadrícula está mal. La cuadrícula hace dos aproximaciones que pueden socavar los resultados. Primero, la cuadrícula es solo un conjunto finito de puntos sobre un subespacio posterior, por lo que es una aproximación aproximada. En segundo lugar, porque es un subespacio finito, estamos declarando falsamente que haya 0 probabilidad posterior sobre los valores de más grande que nuestro valor de cuadrícula grande para . Del mismo modo, es mejor entrar en las colas de la cuadrícula, por lo que sus quanitles de cola son correctos.$N$$N$rstan

---

Necesitaba más espacio para aclarar un punto de la respuesta de Juho. Si entiendo correctamente, podemos integrar fuera de la parte posterior para obtener la distribución beta-binomial: $\theta$

$$p(y|N,\alpha,\beta)={N\choose y} \frac{\text{Beta}(y+\alpha, N-y+\beta)}{\text{Beta}(\alpha,\beta)}$$

En nuestro caso, y porque tenemos un uniforme previo en . Creo que la parte posterior debería ser donde porque . Pero esto parece divergir salvajemente de la respuesta de Juho. ¿Dónde me he equivocado?$\alpha=1$$\beta=1$$\theta$$p(N|y)\propto N^{-1}\prod_{i=1}^K p(y_i|N,\alpha=1,\beta=1)$$K=\#(y)$$p(N)=N^{-1}$

Acantilados: Rstan parece ser (más cercano a) correcto basado en un enfoque que integra fuera analíticamente y evalúa en una cuadrícula bastante grande.$\theta$$P(N)P(y\mid N)$

Para obtener la parte posterior de , en realidad es posible integrar analíticamente: donde es la longitud de . Ahora, dado que tiene Beta anterior (aquí ) y Beta se conjuga con binomial, también sigue una distribución Beta. Por lo tanto, la distribución de es Beta-binomial$N$$\theta$

$$\begin{equation} P(y \mid N) = P(y_1 \mid N)\times P(y_2 \mid N,y_1)\times P(y_3 \mid N,y_1,y_2) \times \ldots P(y_K \mid N,y_1,\ldots,y_{K-1}) \end{equation}$$ $K$$y$$\theta$$Beta(1,1)$$\theta \mid N, y_1, \ldots, y_k$$y_{k+1} \mid N, y_1, \ldots, y_k$para lo cual existe una expresión de forma cerrada de las probabilidades en términos de la función Gamma. Por lo tanto, podemos evaluar calculando los parámetros relevantes de los beta-binomios y multiplicando las probabilidades beta-binomiales. El siguiente código MATLAB utiliza este enfoque para calcular para y se normaliza para obtener el posterior. $P(y\mid N)$$P(N)P(y|N)$$N=72,\ldots,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

El cdf en es , así que supongo que es suficiente, pero uno podría investigar la sensibilidad para aumentar el máximo . El cdf en es solo por lo que su cuadrícula original realmente pierde una masa de probabilidad de cola bastante significativa en comparación con el objetivo de encontrar el cuantil . Los cuantiles que obtengo son $N=100000$$0.9990$maxN=500000$N$$N=5000$$0.9763$$0.975$

$$\begin{equation} \begin{array}{lllll} \textrm{Quantile} & 0.025 & 0.25 & 0.5 & 0.75 & 0.975 \\ N & 95 & 149 & 235 & 478 & 4750 \end{array} \end{equation}$$

Descargo de responsabilidad: no probé mucho el código, puede haber errores (y obviamente también podría haber problemas numéricos con este enfoque). Sin embargo, los cuantiles obtenidos están bastante cerca de los resultados de Rstan, por lo que estoy bastante seguro.