Suma genérica de variables aleatorias Gamma

Question

He leído que la suma de variables aleatorias Gamma con el mismo parámetro de escala es otra variable aleatoria Gamma. También he visto el artículo de Moschopoulos que describe un método para la suma de un conjunto general de variables aleatorias Gamma. He intentado implementar el método de Moschopoulos pero aún no he tenido éxito.

¿Cómo se ve la suma de un conjunto general de variables aleatorias Gamma? Para hacer esta pregunta concreta, ¿cómo se ve para:

$\text{Gamma}(3,1) + \text{Gamma}(4,2) + \text{Gamma}(5,1)$

Si los parámetros anteriores no son particularmente reveladores, por favor sugiera otros.

Answer 1

Primero, combine cualquier suma que tenga el mismo factor de escala: una variable $\Gamma(n, \beta)$ más una variable $\Gamma(m,\beta)$ forman una variable $\Gamma(n+m,\beta)$.

Luego, observe que la función característica (cf) de $\Gamma(n, \beta)$ es $(1-i \beta t)^{-n}$, por lo tanto la cf de la suma de estas distribuciones es el producto

$$\prod_{j} \frac{1}{(1-i \beta_j t)^{n_j}}.$$

Cuando los $n_j$ son todos enteros, este producto se expande como una fracción parcial en una combinación lineal de $(1-i \beta_j t)^{-\nu}$ donde los $\nu$ son enteros entre $1$ y $n_j$. En el ejemplo con $\beta_1 = 1, n_1=8$ (de la suma de $\Gamma(3,1)$ and $\Gamma(5,1)$) y $\beta_2 = 2, n_2=4$ encontramos

$$\begin{aligned}&\frac{1}{(1-i t)^{8}}\frac{1}{(1- 2i t)^{4}} = \\ &\frac{1}{(t+i)^8}-\frac{8 i}{(t+i)^7}-\frac{40}{(t+i)^6}+\frac{160 i}{(t+i)^5}+\frac{560}{(t+i)^4}-\frac{1792 i}{(t+i)^3}\\ &-\frac{5376}{(t+i)^2}+\frac{15360 i}{t+i}+\frac{256}{(2t+i)^4}+\frac{2048 i}{(2 t+i)^3}-\frac{9216}{(2t+i)^2}-\frac{30720 i}{2t+i}. \end{aligned}$$

La inversa de obtener la cf es la Transformada Inversa de Fourier, la cual es lineal: esto significa que podemos aplicarla término por término. Cada término es reconocible como un múltiplo de la cf de una distribución Gamma y por lo tanto se puede invertir fácilmente para obtener la función de densidad de probabilidad (PDF). En el ejemplo obtenemos

$$\begin{aligned} &\frac{e^{-t} t^7}{5040}+\frac{1}{90} e^{-t} t^6+\frac{1}{3} e^{-t} t^5+\frac{20}{3} e^{-t} t^4+\frac{8}{3} e^{-\frac{t}{2}} t^3+\frac{280}{3} e^{-t} t^3\\ &-128 e^{-\frac{t}{2}} t^2+896 e^{-t} t^2+2304 e^{-\frac{t}{2}} t+5376 e^{-t} t-15360 e^{-\frac{t}{2}}+15360 e^{-t} \end{aligned}$$

para la PDF de la suma.

Esto es una mezcla finita de distribuciones Gamma teniendo factores de escala iguales a los dentro de la suma y factores de forma menores o iguales a los dentro de la suma. Excepto en casos especiales (donde podría haber alguna cancelación), el número de términos está dado por el parámetro de forma total $n_1 + n_2 + \cdots$ (asumiendo que todos los $n_j$ son diferentes).

---

Como prueba, aquí hay un histograma de $10^4$ resultados obtenidos al sumar extracciones independientes de las distribuciones $\Gamma(8,1)$ y $\Gamma(4,2)$. Superpuesta en él está el gráfico de $10^4$ veces la función previa. La concordancia es muy buena.

---

Moschopoulos lleva esta idea un paso más adelante al expandir la cf de la suma en una serie infinita de funciones características Gamma siempre que uno o más de los $n_i$ no sea un entero, y luego termina la serie infinita en un punto donde esté razonablemente bien aproximada.

Answer 2

La ecuación de Welch-Satterthwaite podría ser utilizada para dar una respuesta aproximada en forma de una distribución gamma. Esto tiene la buena propiedad de permitirnos tratar las distribuciones gamma como si fueran (aproximadamente) cerradas bajo la adición. Esta es la aproximación en el comúnmente usado t-test de Welch.

(La distribución gamma se puede ver como una distribución chi-cuadrado escalada, permitiendo un parámetro de forma no entero.)

He adaptado la aproximación a la parametrización $k, \theta$ de la distribución gamma:

$$ k_{sum} = { (\sum_i \theta_i k_i)^2 \over \sum_i \theta_i^2 k_i } $$

$$ \theta_{sum} = { { \sum \theta_i k_i } \over k_{sum} } $$

Dejemos $k=(3,4,5)$, $\theta=(1,2,1)$

Así obtenemos aproximadamente Gamma(10.666... ,1.5)

Vemos que el parámetro de forma $k$ ha sido más o menos totalizado, pero ligeramente menos debido a que los parámetros de escala de entrada $\theta_i$ difieren. $\theta$ está configurado de forma que la suma tenga el valor medio correcto.

Answer 3

Mostraré otra posible solución, que es bastante ampliamente aplicable y, con el software R de hoy en día, bastante fácil de implementar. ¡Esa es la aproximación de la densidad del punto de silla, que debería ser más conocida!

Para la terminología sobre la distribución gamma, seguiré https://en.wikipedia.org/wiki/Gamma_distribution con la parametrización de forma/escala, $k$ es el parámetro de forma y $\theta$ es la escala. Para la aproximación de punto de silla seguiré a Ronald W Butler: "Saddlepoint approximations with applications" (Cambridge UP). La aproximación de punto de silla se explica aquí: ¿Cómo funciona la aproximación de punto de silla? aquí mostraré cómo se usa en esta aplicación.

Sea $X$ una variable aleatoria con función de generación de momentos existente $$ M(s) = E e^{sX} $$ que debe existir para $s$ en algún intervalo abierto que contenga cero. Luego define la función generadora de cumulantes por $$ K(s) = \log M(s) $$ Se sabe que $E X = K'(0), \text{Var} (X) = K''(0)$. La ecuación de punto de silla es $$ K'(\hat{s}) = x$$ que define implícitamente $s$ como una función de $x$ (que debe estar en el rango de $X$). Escribimos esta función definida implícitamente como $\hat{s}(x)$. Notar que la ecuación de punto de silla siempre tiene exactamente una solución, porque la función cumulante es convexa.

Luego, la aproximación de punto de silla a la densidad $f$ de $X$ está dada por $$ \hat{f}(x) = \frac1{\sqrt{2\pi K''(\hat{s})}} \exp(K(\hat{s}) - \hat{s} x) $$ Esta función de densidad aproximada no garantiza integrar a 1, por lo que es la aproximación de punto de silla no normalizada. Podríamos integrar numéricamente y luego normalizar para obtener una mejor aproximación. Pero esta aproximación está garantizada de ser no negativa.

Ahora sean $X_1, X_2, \dots, X_n$ variables aleatorias gamma independientes, donde $X_i$ tiene la distribución con parámetros $(k_i, \theta_i)$. Luego la función generadora de cumulantes es $$ K(s) = -\sum_{i=1}^n k_i \ln(1-\theta_i s) $$ definida para $s<1/\max(\theta_1, \theta_2, \dots, \theta_n)$. La primera derivada es $$ K'(s) = \sum_{i=1}^n \frac{k_i \theta_i}{1-\theta_i s} $$ y la segunda derivada es $$ K''(s) = \sum_{i=1}^n \frac{k_i \theta_i^2}{(1-\theta_i s)^2}. $$ En lo siguiente daré algo de código R que calcula esto, y usaré los valores de parámetros $n=3$, $k=(1,2,3)$, $\theta=(1,2,3)$. Notar que el siguiente código R utiliza un nuevo argumento en la función uniroot presentado en R 3.1, por lo que no se ejecutará en versiones anteriores de R.

shape <- 1:3 #ki
scale <- 1:3 # thetai
    # Para este caso, obtenemos expectativa=14, varianza=36
make_cumgenfun  <-  function(shape, scale) {
    # retornamos una lista(shape, scale, K, K', K'')
    n  <-  length(shape)
    m <-   length(scale)
    stopifnot( n == m, shape > 0, scale > 0 )
    return( list( shape=shape,  scale=scale, 
          Vectorize(function(s) {-sum(shape * log(1-scale * s) ) }),
          Vectorize(function(s) {sum((shape*scale)/(1-s*scale))}) ,
          Vectorize(function(s) { sum(shape*scale*scale/(1-s*scale))
      }))    )
    }

solve_speq  <-  function(x, cumgenfun) {
       # ¡Devuelve punto de silla!
       shape <- cumgenfun[[1]]
       scale <- cumgenfun[[2]]
       Kd  <-   cumgenfun[[4]]
       uniroot(function(s) Kd(s)-x,lower=-100,
                      upper = 0.3333, 
                      extendInt = "upX")$root
    }

make_fhat <-  function(shape,  scale) {
    cgf1  <-  make_cumgenfun(shape, scale)
    K  <-  cgf1[[3]]
    Kd <-  cgf1[[4]]
    Kdd <- cgf1[[5]]
    # Función para encontrar fhat para un x específico:
    fhat0  <- function(x) {
    # Resolver la ecuación de punto de silla:
    s  <-  solve_speq(x, cgf1)
    # Calcular el valor de densidad de punto de silla:
    (1/sqrt(2*pi*Kdd(s)))*exp(K(s)-s*x)
        }
        # Devolver una versión vectorizada:
   return(Vectorize(fhat0))
} #end make_fhat

fhat  <-  make_fhat(shape, scale)
plot(fhat, from=0.01,  to=40, col="red", 
         main="aproximación de punto de silla no normalizada\n
               a la suma de tres variables gamma")

resultando en el siguiente gráfico:

Dejaré la aproximación de punto de silla normalizada como ejercicio.

Answer 4

Una solución exacta para la convolución (es decir, la suma) de $n$ distribuciones gamma se da como Eq. (1) en el pdf vinculado por DiSalvo. Dado que esto es un poco largo, llevará algún tiempo copiarlo aquí. Para solo dos distribuciones gamma, su suma exacta en forma cerrada está especificada por la Eq. (2) de DiSalvo y sin pesos por la Eq. (5) de Wesolowski et al., que también aparece en el sitio de CV como respuesta a esa pregunta. Es decir, $$\mathrm{G}\mathrm{D}\mathrm{C}\left(\mathrm{a}\kern0.1em ,\mathrm{b}\kern0.1em ,\alpha, \beta; \tau \right)=\left\{\begin{array}{cc}\hfill \frac{{\mathrm{b}}^{\mathrm{a}}{\beta}^{\alpha }}{\Gamma \left(\mathrm{a}+\alpha \right)}{e}^{-\mathrm{b}\tau }{\tau^{\mathrm{a}+\alpha-1}}{}_1F_1\left[\alpha, \mathrm{a}+\alpha, \left(\mathrm{b}-\beta \right)\tau \right],\hfill & \hfill \tau >0\hfill \\ {}\hfill \kern2em 0\kern6.6em ,\hfill \kern5.4em \tau \kern0.30em \le \kern0.30em 0\hfill \end{array}\right.,$$ donde la notación en las preguntas anteriores; $Gamma(a,b) \rightarrow \Gamma(a,1/b)$, aquí. Es decir, $b$ y $\beta$ son constantes de tasa aquí y no escalares de tiempo.

Answer 5

Según Ansari et al. 2012, la PDF y CDF de variables aleatorias gamma independientes con diferentes distribuciones pueden expresarse en términos de la función H de Fox (función H-bar). El documento referenciado a continuación también contiene una implementación de esta función en el lenguaje Wolfram.

REFERENCIAS:

• Ansari, I.S., Yilmaz, F., et al., "Nuevos Resultados sobre la Suma de Variables Aleatorias Gamma con Aplicación al Rendimiento de Sistemas de Comunicación Inalámbrica sobre Canales de Desvanecimiento Nakagami-m", arXiv:1202.2576 [cs.IT], 2012.