Encontrar una manera de simular números aleatorios para esta distribución

Question

Estoy intentando escribir un programa en R que simule números pseudoaleatorios a partir de una distribución con la función de distribución acumulativa:

$$F(x)= 1-\exp \left(-ax-\frac{b}{p+1}x^{p+1}\right), \quad x \geq 0$$

donde $a,b>0, p \in (0,1)$

He probado el muestreo por transformada inversa, pero la inversa no parece resoluble analíticamente. Me gustaría que me sugirieran una solución a este problema.

Answer 1

Hay una solución sencilla (y si se me permite añadir, elegante) para este ejercicio: puesto que $1-F(x)$ parece un producto de dos distribuciones de supervivencia: $$(1-F(x))=\exp\left\{-ax-\frac{b}{p+1}x^{p+1}\right\}=\underbrace{\exp\left\{-ax\right\}}_{1-F_1(x)}\underbrace{\exp\left\{-\frac{b}{p+1}x^{p+1}\right\}}_{1-F_2(x)}$$ la distribución $F$ es la distribución de $$X=\min\{X_1,X_2\}\qquad X_1\sim F_1\,,X_2\sim F_2$$ En este caso $F_1$ es el exponencial $\mathcal{E}(a)$ distribución y $F_2$ es el $1/(p+1)$ -ésima potencia de un exponencial $\mathcal{E}(b/(p+1))$ distribución.

El código R asociado no puede ser más sencillo

x=pmin(rexp(n,a),rexp(n,b/(p+1))^(1/(p+1))) #simulating an n-sample

y es definitivamente mucho más rápido que las resoluciones inversa pdf y aceptar-rechazar:

> n=1e6
> system.time(results <- Vectorize(simulate,"prob")(runif(n)))
utilisateur     système      écoulé 
    89.060       0.072      89.124 
> system.time(x <- simuF(n,1,2,3))
utilisateur     système      écoulé 
     1.080       0.020       1.103 
> system.time(x <- pmin(rexp(n,a),rexp(n,b/(p+1))^(1/(p+1))))
utilisateur     système      écoulé 
     0.160       0.000       0.163 

con un ajuste sorprendentemente perfecto:

Answer 2

Siempre se puede resolver numéricamente la transformación inversa.

A continuación, realizo una búsqueda de bisección muy sencilla. Para una probabilidad de entrada dada $q$ (Yo uso $q$ puesto que ya tiene un $p$ en su fórmula), empiezo con $x_L=0$ et $x_R=1$ . Entonces doblo $x_R$ hasta $F(x_R)>q$ . Por último, biseco iterativamente el intervalo $[x_L,x_R]$ hasta que su longitud sea inferior a $\epsilon$ y su punto medio $x_M$ satisface $F(x_M)\approx q$ .

El ECDF se adapta a sus $F$ lo suficientemente bien para mis elecciones de $a$ et $b$ y es razonablemente rápido. Probablemente se podría acelerar utilizando alguna optimización de tipo Newton en lugar de la simple búsqueda de bisección.

aa <- 2
bb <- 1
pp <- 0.1

cdf <- function(x) 1-exp(-aa*x-bb*x^(pp+1)/(pp+1))

simulate <- function(prob,epsilon=1e-5) {
    left <- 0
    right <- 1
    while ( cdf(right) < prob ) right <- 2*right

    while ( right-left>epsilon ) {
        middle <- mean(c(left,right))
        value_middle <- cdf(middle)
        if ( value_middle < prob ) left <- middle else right <- middle
    }

    mean(c(left,right))
}

set.seed(1)
results <- Vectorize(simulate,"prob")(runif(10000))
hist(results)

xx <- seq(0,max(results),by=.01)
plot(ecdf(results))
lines(xx,cdf(xx),col="red")

Answer 3

Existe una resolución algo enrevesada aunque directa por aceptación-rechazo. En primer lugar, una simple diferenciación muestra que la pdf de la distribución es $$f(x)=(a+bx^p)\exp\left\{-ax-\frac{b}{p+1}x^{p+1}\right\}$$ En segundo lugar, puesto que $$f(x)=ae^{-ax}\underbrace{e^{-bx^{p+1}/(p+1)}}_{\le 1}+bx^pe^{-bx^{p+1}/(p+1)}\underbrace{e^{-ax}}_{\le 1}$$ tenemos el límite superior $$f(x)\le g(x)=ae^{-ax}+bx^pe^{-bx^{p+1}/(p+1)}$$ En tercer lugar, considerando el segundo término en $g$ tomar el cambio de variable $\xi=x^{p+1}$ es decir, $x=\xi^{1/(p+1)}$ . Entonces $$\dfrac{\text{d}x}{\text{d}\xi}=\dfrac{1}{p+1}\xi^{\frac{1}{p+1}-1}=\dfrac{1}{p+1}\xi^{\frac{-p}{p+1}}$$ es el jacobiano del cambio de variable. Si $X$ tiene una densidad de la forma $\kappa bx^pe^{-bx^{p+1}/(p+1)}$ donde $\kappa$ es la constante de normalización, entonces $\Xi=X^{1/(p+1)}$ tiene la densidad $$\kappa b\xi^{\frac{p}{p+1}}e^{-b\xi/(p+1)}\,\dfrac{1}{p+1}\xi^{\frac{-p}{p+1}}=\kappa \dfrac{b}{p+1}e^{-b\xi/(p+1)}$$ lo que significa que (i) $\Xi$ se distribuye como una exponencial $\mathcal{E}(b/(p+1))$ y (ii) la constante $\kappa$ es igual a uno. Por lo tanto, $g(x)$ acaba siendo igual a la mezcla igualmente ponderada de un Exponencial $\mathcal{E}(a)$ distribución y la $1/(p+1)$ -ésima potencia de un exponencial $\mathcal{E}(b/(p+1))$ modulo una constante multiplicativa de $2$ para tener en cuenta los pesos: $$f(x)\le g(x)=2\left(\frac{1}{2} ae^{-ax}+\frac{1}{2} bx^pe^{-bx^{p+1}/(p+1)}\right)$$ Y $g$ es fácil de simular como una mezcla.

Una representación en R del algoritmo de aceptación-rechazo es la siguiente

simuF <- function(a,b,p){
  reepeat=TRUE
  while (reepeat){
   if (runif(1)<.5) x=rexp(1,a) else
      x=rexp(1,b/(p+1))^(1/(p+1))
   reepeat=(runif(1)>(a+b*x^p)*exp(-a*x-b*x^(p+1)/(p+1))/
      (a*exp(-a*x)+b*x^p*exp(-b*x^(p+1)/(p+1))))}
  return(x)}

y para una muestra n:

simuF <- function(n,a,b,p){
  sampl=NULL
  while (length(sampl)<n){
   x=u=sample(0:1,n,rep=TRUE)
   x[u==0]=rexp(sum(u==0),b/(p+1))^(1/(p+1))
   x[u==1]=rexp(sum(u==1),a)
   sampl=c(sampl,x[runif(n)<(a+b*x^p)*exp(-a*x-b*x^(p+1)/(p+1))/
      (a*exp(-a*x)+b*x^p*exp(-b*x^(p+1)/(p+1)))])
   }
  return(sampl[1:n])}

He aquí una ilustración para a=1, b=2, p=3: