¿Cómo muestrear rápidamente X si exp (X) ~ Gamma?

Tengo un problema de muestreo simple, donde mi bucle interno se ve así:

v = sample_gamma(k, a)

donde sample_gammamuestras de la distribución Gamma para formar una muestra de Dirichlet.

Funciona bien, pero para algunos valores de k / a, algunos de los flujos de cálculo aguas abajo.

Lo adapté para usar variables de espacio de registro:

v = log(sample_gamma(k, a))

Después de adaptar todo el resto del programa, funciona correctamente (al menos me da los mismos resultados exactos en casos de prueba). Sin embargo, es más lento que antes.

¿Hay alguna forma de muestrear directamente sin usar funciones lentas como ? Intenté buscar en Google para esto, pero ni siquiera sé si esta distribución tiene un nombre común (log-gamma?).Registro gamma ( $X, \exp(X) \sim \text{Gamma}$$\log()$

Considere una pequeña forma parámetro cerca de 0, tal como α = 1 / 100 . En el rango entre 0 y α , e - α es aproximadamente 1 , por lo que el Gamma pdf es aproximadamente x α - 1 d x / Γ ( α ) . Esto se puede integrar a un CDF aproximado, F α ( x ) = x $\alpha$$\alpha = 1/100$$\alpha$$e^{-\alpha}$$1$$x^{\alpha-1}dx / \Gamma(\alpha)$ . Al invertirlo, vemos unapotencia1/α: un gran exponente. Paraα=1/100esto causa alguna posibilidad de underflow (un valor de doble precisión de menos de10-300, más o menos). Aquí hay una gráfica de la posibilidad de obtener un flujo inferior en función del logaritmo de base diez deα:$F_\alpha(x) = \frac{x^\alpha}{\alpha \Gamma(\alpha)}$$1/\alpha$$\alpha = 1/100$$10^{-300}$$\alpha$

Una solución es explotar esta aproximación para generar variables log (Gamma): en efecto, intente generar una variante Gamma y, si es demasiado pequeña, generar su logaritmo a partir de esta distribución de potencia aproximada (como se muestra a continuación). (Haga esto repetidamente hasta que el registro esté dentro del rango de flujo inferior, de modo que sea un sustituto válido para la variante de flujo inferior original). Para el cálculo de Dirichlet, reste el máximo de todos los logaritmos de cada uno de los valores de registro: esto reescala implícitamente todos Gamma varía para que no afecte los valores de Dirichlet. Trate cualquier registro resultante que sea demasiado pequeño (digamos, menor que -100) como el registro de un cero verdadero. Exponga los otros registros. Ahora puede continuar sin desbordamiento.

Esto llevará más tiempo que antes, ¡pero al menos funcionará!

$\alpha$$C = \log(\Gamma(\alpha)) + \log(\alpha)$$\alpha$$C$

Debido a que el parámetro de escala simplemente reescala la variante, no hay problema para acomodarla en estos procedimientos. Ni siquiera lo necesita si todos los parámetros de escala son iguales.

Editar

$1/\alpha$$B(\alpha)$$\Gamma(\alpha+1)$$\left(\alpha x^{\alpha-1}\right) \left(y^{\alpha}e^{-y}dy/\Gamma(\alpha+1)\right)$$z = xy$$y \to z/x$$x$$x$$z$$\infty$$0 \le y \le 1$

$\Gamma(\alpha)$

$0 \lt \alpha \lt 1$$\Gamma(\alpha+1)$$1/\alpha$$\Gamma(\alpha)$

ra$\alpha$b$\beta$

  if (a < 1)
    {
      double u = gsl_rng_uniform_pos (r);
      return gsl_ran_gamma (r, 1.0 + a, b) * pow (u, 1.0 / a);
   }

gsl_ran_gammagsl_rng_uniform_pos$(0,1)$_pos

Por lo tanto, puedo tomar el registro de la última expresión y usar

return log(gsl_ran_gamma(r, 1.0 + a, b)) + log(u)/a;

log()pow()$1/a$$1/a$