¿Por qué es rechazo de muestreo con probabilidad 2/3 de la aceptación de Beta(2,2) no más lento que 'rbeta(N,2,2)'?

Question

Yo estaba tratando de ilustrar que el rechazo de muestreo es ineficiente cuando un enfoque alternativo que no bote las muestras. (Este post obviamente es un candidato para la migración, pero creo que también podría ser una estadística aspecto.)

Mi ejemplo de que iba a ser la simulación de $Beta(2,2)$ varia con un algoritmo de rechazo y rbeta(N,2,2). Mi código es el siguiente:

# densityatmode <- dbeta(.5,2,2)=1.5
rejectionsampling_beta22 <- function(N){
     px <- runif(N,min=0,max=1)
     py <- runif(N,min=0,max=1.5)
     ikeep <- (py < 6*px*(1-px))
     return(px[ikeep])
}

samples <- 1500000

La probabilidad de aceptación es area under density/area in box, area under density de curso igual a uno, por lo 1/1.5, o, equivalentemente, tenemos 1.5 como muchos de los candidatos de las muestras en el uso de rechazo.

library(microbenchmark)
result <- microbenchmark(rejectionsampling_beta22(1.5*samples),rbeta(samples,2,2))

El resultado, lo cual es sorprendente para mí, ya que yo sentía que el rechazo de 1/3 de los sorteos debe ser un desperdicio:

Warning message:
In microbenchmark(rejectionsampling_beta22(1.5 * samples), rbeta(samples,  :
  Could not measure a positive execution time for 33 evaluations.
> result
Unit: nanoseconds
                                    expr       min        lq         mean    median        uq       max neval cld
 rejectionsampling_beta22(1.5 * samples) 180171443 254891504 248482870.58 257005186 258866373 329533988   100  b 
                    rbeta(samples, 2, 2) 253091893 254717640 261602707.69 257391248 259902160 333055032   100   c
                                   neval         0         0        63.97         1         1       604   100 a  

Answer 1

Aquí es Cheng & Fiesta Gamma generador de código, en el que R rbeta y rgamma funciones son:

function x=gamrnd_cheng(alpha)
% Gamma(alpha,1) generator using Cheng--Feast method
% Algorithm 4.35
c1=alpha-1; c2=(alpha-1/(6*alpha))/c1; c3=2/c1; c4=1+c3;
c5=1/sqrt(alpha);
flag=0;
while flag==0;
    U1=rand; U2=rand;
    if alpha>2.5
        U1=U2+c5*(1-1.86*U1);
    end
       W=c2*U2/U1;
     flag=(U1<1)&&(U1>0)&&(((c3*U1+W+1/W)<c4)||((c3*log(U1)-log(W)+W)<1));
end
x=c1*W;

que puede ser reciclado en un Beta generador en aproximadamente el mismo costo. Utiliza dos uniformes, además de un rechazo de la condición, así que para los valores de $(\alpha,\beta)$ que eligió, es decir, por un rechazo de la probabilidad de $1/3$, el aceptar-rechazar algoritmo puede ser igual de eficiente. Sin embargo, también debe ejecutar la comparación para grandes valores no enteros de $(\alpha,\beta)$ a comprobar si es o no el Cheng & Fiesta Gamma generador sigue siendo eficaz.

Por ejemplo, Joe Whittaker Beta $\mathfrak{B}(\alpha,\beta)$ generador tiene un rechazo de la condición de la forma$$U_1^{1/\alpha}+U_2^{1/\beta}>1$$which occurs with increasing frequency as $\alfa$ and $\beta$ increase. I remember Luc Devroye mentioning that, for $\mathfrak{G}(\alpha,1)$ distributions, it is not possible to find a bound on the computing time that is independent of $\alpha$...

Answer 2

El seguimiento de Xi'an, la sugerencia, el resultado en la cuestión de hecho parece ser un artefacto de la (pequeña) parámetros escogidos. En particular, cuando la elección de mayor$\alpha$$\beta$, de tal manera que la beta de la densidad se vuelve más concentrada con un correspondientemente mayor densidad en el modo, de tal manera que un cuadro más grande debe ser construido alrededor de la beta de la densidad en el rechazo de un algoritmo, de tal manera que la probabilidad de aceptación disminuye, el resultado desaparece y rbeta es mucho más eficiente.

Aquí es un ejemplo para $\alpha=\beta=10.2$:

# densityatmode <- dbeta(.5,10.2,10.2)=3.559877
# precompute the beta density:
# gamma(20.4)/gamma(10.2)^2=1231365
rejectionsampling_betaab <- function(N){
  px <- runif(N,min=0,max=1)
  py <- runif(N,min=0,max=3.559877)
  #ikeep <- (py < dbeta(px,10.2,10.2))
  ikeep <- (py < 1231365*px^9.2*(1-px)^9.2)
  return(px[ikeep])
}

AcceptanceProbability <- 1/(3.559877)
samples <- 1e5

library(microbenchmark)
result <- microbenchmark(rejectionsampling_betaab(1/AcceptanceProbability*samples),rbeta(samples,10.2,10.2))

El resultado ahora claramente favorece rbeta:

> result
Unit: milliseconds
                                                        expr       min        lq      mean    median        uq      max neval
 rejectionsampling_betaab(1/AcceptanceProbability * samples) 417.46861 426.16041 469.10049 430.37921 485.60090 635.9086   100
                                  rbeta(samples, 10.2, 10.2)  90.39748  90.94824  94.24759  91.52765  92.51329 264.5119   100