¿Cómo muestra de $c^a d^{a-1} / \Gamma(a)$?

Question

Quiero la muestra de acuerdo a una densidad de $$ f(a) \propto \frac{c^d^{- 1}}{\Gamma(un)} 1_{(1,\infty)}(a) $$ donde $c$ $d$ son estrictamente positivos. (Motivación: Esto podría ser útil para el muestreo de Gibbs, cuando la forma de un parámetro Gamma, densidad homogénea antes.)

¿Alguien sabe como muestra de esta densidad fácilmente? Tal vez es estándar y sólo algo que yo no sé?

No puedo pensar en un estúpido rechazo sampliing algoritmo que va más o menos trabajo (encontrar el modo de $a^*$$f$, de la muestra $(a,u)$ de uniforme en un gran cuadro de $[0,10a^*]\times [0,f(a^*)]$ y rechazar si $u>f(a)$), pero (i) no es del todo eficiente y (ii) $f(a^*)$ será demasiado grande para que un equipo pueda manejar fácilmente, incluso moderadamente grande $c$ $d$. (Tenga en cuenta que el modo para un gran $c$ $d$ es aproximadamente a $a=cd$.)

Gracias de antemano por cualquier ayuda!

Answer 1

Me gusta @whuber la respuesta de mucho; es probable que sea muy eficiente y tiene un hermoso análisis. Pero se requiere de un profundo conocimiento con respecto a esta distribución particular. Para situaciones en las que no tienen esta visión (así que, por diferentes distribuciones), también me gusta el enfoque que funciona para todas las distribuciones donde el PDF es dos veces diferenciable y que la segunda derivada tiene un número finito de raíces. Se requiere un poco de trabajo para establecer, pero luego te tiene un motor que funciona para la mayoría de las distribuciones que usted puede lanzar en él.

Básicamente, la idea es usar un modelo lineal por tramos límite superior para el PDF que adaptar lo que está haciendo el rechazo de muestreo. Al mismo tiempo, usted tiene un modelo lineal por tramos inferiores obligado para el PDF, lo cual evita tener que evaluar el PDF con demasiada frecuencia. Los límites superior e inferior están dadas por los acordes y las tangentes a la PDF gráfico. La división inicial en intervalos es tal que en cada intervalo, el PDF es cóncava o todos convexo; cada vez que tiene que rechazar un punto (x, y) que subdividir el intervalo en x. (También se puede hacer una subdivisión adicional en x si tuviera que calcular el PDF porque el límite inferior es realmente malo.) Esto hace que las subdivisiones ocurrir especialmente con frecuencia donde la parte superior (inferior) los límites son malos, por lo que se obtiene una buena aproximación de su PDF esencialmente de forma gratuita. Los detalles son un poco difícil de conseguir, pero he tratado de explicar la mayoría de ellos en esta serie de blog posts - especialmente el último.

Los mensajes no discutir qué hacer si el PDF es ilimitado, ya sea en el dominio o en valores; me gustaría recomendar un poco la solución obvia de hacer una transformación que los convierte en finito (que sería difícil de automatizar) o el uso de un corte. Yo elegiría el punto de corte en función del número total de puntos que esperar a generar, digamos N, y elegir el punto de corte de manera que la parte eliminada tiene menos de $1 / (10 N)$ de probabilidad. (Esto es bastante fácil si usted tiene una forma cerrada para el CDF; de lo contrario también puede ser complicado.)

Este método está implementado en Arce como el método por defecto para el usuario definido por el continuo distribuciones. (Revelación completa - yo trabajo para Maplesoft.)

---

Hice un ejemplo de ejecución, generando 10^4 puntos para c = 2, d = 3, la especificación [1, 100] como el primer rango para los valores:

Hubo 23 de rechazos (en rojo), 51 puntos "en libertad", que en el tiempo entre el límite inferior y el PDF real, y 9949 puntos que fueron aceptados después de comprobar sólo desigualdades lineales. Que el 74 evaluaciones de los PDF en total, alrededor de un PDF de evaluación por 135 puntos. La relación debe obtener mejores generar más puntos, ya que la aproximación se pone mejor y mejor (y a la inversa, si se genera sólo unos pocos puntos, la proporción es de lo peor).

Answer 2

El rechazo de muestreo va a funcionar excepcionalmente bien al $c d \ge \exp(5)$ y es razonable para $c d \ge \exp(2)$.

Para simplificar las matemáticas un poco, vamos a $k = c d$, escribir $x = a$, y la nota que

$$f(x) \propto \frac{k^x}{\Gamma(x)} dx$$

para $x \ge 1$. Establecimiento $x = u^{3/2}$ da

$$f(u) \propto \frac{k^{u^{3/2}}}{\Gamma(u^{3/2})} u^{1/2} du$$

para $u \ge 1$. Al $k \ge \exp(5)$, esta distribución es muy cercano a lo Normal (y se acerca como $k$ se hace más grande). Específicamente, usted puede

1. Encontrar el modo de $f(u)$ numéricamente (utilizando, por ejemplo, Newton-Raphson).

2. Expandir $\log{f(u)}$ a de segundo orden acerca de su modo.

Este rendimientos de los parámetros de un aproximan a la distribución Normal. A la alta precisión, esta aproximación Normal domina $f(u)$, excepto en el extremo de la cola. (Al $k \lt \exp(5)$, usted puede necesitar para escalar el pdf Normal hasta un poco para asegurar la dominación.)

Después de haber hecho este trabajo preliminar para cualquier valor dado de a $k$, y habiendo estimado una constante $M \gt 1$ (como se describe a continuación), la obtención de un azar de la variable aleatoria es una cuestión de:

1. Dibujar un valor de $u$ desde la que domina la distribución Normal $g(u)$.

2. Si $u \lt 1$ o si un nuevo uniforme de la variable aleatoria $X$ supera $f(u)/(M g(u))$, volver al paso 1.

3. Set $x = u^{3/2}$.

El número esperado de las evaluaciones de las $f$ debido a las discrepancias entre las $g$ $f$ es sólo ligeramente mayor que 1. (Algunas evaluaciones adicionales se producen debido a los rechazos de variables menos de $1$, pero incluso cuando se $k$ es tan baja como $2$ la frecuencia de este tipo de casos es pequeño.)

Este gráfico muestra el logaritmos de g y f como una función de u de $k=\exp(5)$. Debido a que los gráficos están tan cerca, tenemos que inspeccionar su relación para ver lo que está pasando:

Esto muestra el registro de la relación de $\log(\exp(0.004)g(u)/f(u))$; el factor de $M = \exp(0.004)$ fue incluido para asegurar que el logaritmo es positivo a lo largo de la parte principal de la distribución; es decir, garantizar a $Mg(u) \ge f(u)$, excepto, posiblemente, en las regiones de insignificante probabilidad. Haciendo $M$ suficientemente grande puede garantizar que $M \cdot g$ domina $f$ en todos, pero la mayoría de los extremos (que no tiene prácticamente ninguna posibilidad de ser elegido en una simulación de todos modos). Sin embargo, el mayor $M$ es el más frecuentemente rechazos va a ocurrir. Como $k$ crece grande, $M$ puede ser elegido muy cerca de $1$, que se incurre en prácticamente ningún tipo de penalización.

Un enfoque similar funciona incluso para $k \gt \exp(2)$, pero bastante grande, los valores de $M$ puede ser necesario cuando se $\exp(2) \lt k \lt \exp(5)$, debido a $f(u)$ es notablemente asimétrica. Por ejemplo, con $k = \exp(2)$, para obtener una razonablemente precisas $g$ necesitamos establecer $M=1$:

El rojo superior de la curva es la gráfica de $\log(\exp(1)g(u))$, mientras que la parte inferior de la curva azul es la gráfica de $\log(f(u))$. El rechazo de muestreo de $f$ en relación al $\exp(1)g$ hará que alrededor de 2/3 de todo el juicio llega a ser rechazado, triplicando el esfuerzo: todavía no es malo. La cola derecha ($u \gt 10$ o $x \gt 10^{3/2} \sim 30$) serán representadas en el rechazo de muestreo (debido a $\exp(1)g$ ya no domina $f$), sino que la cola se compone de menos de $\exp(-20) \sim 10^{-9}$ del total de la probabilidad.

Para resumir, después de un esfuerzo inicial para calcular el modo y evaluar el término cuadrático de la alimentación de la serie de $f(u)$ todo el modo-un esfuerzo que requiere de un par de decenas de función de la evaluación en la mayoría de los--usted puede utilizar el rechazo de muestreo por un costo estimado de entre 1 y 3 (o así) de evaluaciones por la variable aleatoria. El costo multiplicador desciende rápidamente a 1 $k = c d$ aumenta más de 5.

Incluso cuando sólo uno de los sorteo de $f$ es necesario, este método es razonable. Se trata en su propia cuando muchos sorteos son necesarios para el mismo valor de $k$, para entonces, la sobrecarga de los cálculos iniciales se amortiza a lo largo de muchos empates.

---

Anexo

@Cardenal ha pedido, bastante razonablemente, para el apoyo de algunos de la mano saludando análisis en el anterior. En particular, ¿por qué debe la transformación de $x = u^{3/2}$ hacer la distribución aproximadamente Normal?

A la luz de la teoría de Box-Cox transformaciones, es natural buscar algún poder de transformación de la forma $x = u^\alpha$ (para una constante $\alpha$, esperemos que no demasiado diferente de la unidad) que va a hacer una distribución "más" de lo Normal. Recordar que todas las distribuciones Normales son simplemente caracteriza: los logaritmos de sus archivos pdf son puramente cuadrática, con cero término lineal y no los términos de orden superior. Por lo tanto, podemos tomar cualquier pdf y compararla a una distribución Normal mediante la ampliación de su logaritmo como una potencia de la serie alrededor de su (mayor) de pico. Buscamos un valor de $\alpha$ que hace (al menos) el tercer poder desaparecer, al menos aproximadamente: que es lo más que podemos esperar razonablemente que un solo gratis coeficiente será lograr. A menudo, esto funciona bien.

Pero, ¿cómo conseguir una manija en esta distribución particular? Al efectuar la transformación de energía, su pdf

$$f(u) = \frac{k^{u^{\alpha}}}{\Gamma(u^{\alpha})} u^{\alpha-1}.$$

Tome su logaritmo y el uso de Stirling asintótica de expansión de $\log(\Gamma)$:

$$\log(f(u)) \approx \log(k) u^\alpha + (\alpha - 1)\log(u) - \alpha u^\alpha \log(u) + u^\alpha - \log(2 \pi u^\alpha)/2 + c u^{-\alpha}$$

(para valores pequeños de a $c$, que es no constante). Esto funciona siempre $\alpha$ es positivo, lo que nos va a suponer a ser el caso (de lo contrario, no podemos descuidar el resto de la expansión).

Calcular su tercera derivada (que, cuando se divide por $3!$, será el coeficiente de la tercera potencia de $u$ en el poder de la serie) y aprovechar el hecho de que en el pico, la primera derivada debe ser cero. Esto simplifica la tercera derivada en gran medida, de dar (aproximadamente, porque nos están haciendo caso omiso de la derivada de $c$)

$$-\frac{1}{2} u^{-(3+\alpha)} \alpha \left(2 \alpha(2 \alpha-3) u^{2 \alpha} + (\alpha^2 - 5\alpha +6)u^\alpha + 12 c \alpha \right).$$

Al $k$ no es demasiado pequeño, $u$ hecho va a ser grande en el pico. Debido a $\alpha$ es positivo, el dominante término de esta expresión es el $2\alpha$ de la potencia, de la que podemos establecer a cero, haciendo que su coeficiente de desaparecer:

$$2 \alpha-3 = 0.$$

Por eso $\alpha = 3/2$ funciona tan bien: con esta opción, el coeficiente de la cúbico plazo alrededor de la cima se comporta como $u^{-3}$, que está cerca de a $\exp(-2 k)$. Una vez $k$ supera los 10 o así, que prácticamente se puede olvidarse de él, y es razonablemente pequeña, incluso para $k$ a 2. Los poderes superiores, a partir del cuarto, jugar menos y menos de un papel como $k$ se hace grande, debido a que sus coeficientes de crecer proporcionalmente más pequeñas, también. Por cierto, los mismos cálculos (basados en la segunda derivada de $log(f(u))$ en su pico) muestran la desviación estándar de esta aproximación Normal es ligeramente menor que $\frac{2}{3}\exp(k/6)$, con el error proporcional a $\exp(-k/2)$.

Answer 3

Usted podría hacerlo por numéricamente la ejecución de la inversión método, que dice que si conecta uniforme(0,1) variables aleatorias en la inversa de la CDF, consigue un empate de la distribución. He incluido algunas R siguiente código que hace esto, y de las pocas comprobaciones que he hecho, se está trabajando bien, pero es un poco descuidado y estoy seguro que se puede optimizar.

Si usted no está familiarizado con R, lgamma() es el logaritmo de la función gamma; integrar() calcula el definitivo 1-D integral; uniroot() calcula una raíz de una función utilizando 1-D interseccion.

# density. using the log-gamma gives a more numerically stable return for 
# the subsequent numerical integration (will not work without this trick)
f = function(x,c,d) exp( x*log(c) + (x-1)*log(d) - lgamma(x) )

# brute force calculation of the CDF, calculating the normalizing constant numerically
F = function(x,c,d) 
{
   g = function(x) f(x,c,d)
   return( integrate(g,1,x)$val/integrate(g,1,Inf)$val )
}

# Using bisection to find where the CDF equals p, to give the inverse CDF. This works 
# since the density given in the problem corresponds to a continuous CDF. 
F_1 = function(p,c,d) 
{
   Q = function(x) F(x,c,d)-p
   return( uniroot(Q, c(1+1e-10, 1e4))$root )
}

# plug uniform(0,1)'s into the inverse CDF. Testing for c=3, d=4. 
G = function(x) F_1(x,3,4)
z = sapply(runif(1000),G)

# simulated mean
mean(z)
[1] 13.10915

# exact mean
g = function(x) f(x,3,4)
nc = integrate(g,1,Inf)$val
h = function(x) f(x,3,4)*x/nc
integrate(h,1,Inf)$val
[1] 13.00002 

# simulated second moment
mean(z^2)
[1] 183.0266

# exact second moment
g = function(x) f(x,3,4)
nc = integrate(g,1,Inf)$val
h = function(x) f(x,3,4)*(x^2)/nc
integrate(h,1,Inf)$val
[1] 181.0003

# estimated density from the sample
plot(density(z))

# true density 
s = seq(1,25,length=1000)
plot(s, f(s,3,4), type="l", lwd=3)

El principal arbitraria cosa que voy a hacer aquí es, asumiendo que $(1,10000)$ es suficiente soporte para la interseccion - yo era perezoso acerca de esto y podría haber una manera más eficiente para elegir este soporte. Para valores muy grandes, el cálculo numérico de la CDF (por ejemplo, $> 100000$) falla, de manera que el soporte debe estar por debajo de esta. El CDF es efectivamente igual a 1 en los puntos (a menos que $c, d$ son muy grandes), así que algo probablemente podría ser incluido que evitar el error de cálculo de la CDF para grandes valores de entrada.

Edit: Al $cd$ es muy grande, una numérica de un problema con este método. Como whuber señala en los comentarios, una vez que esto ha ocurrido, la distribución es esencialmente degenerar a su modo, lo que es un simple problema de muestreo.