¿Cómo se distribuyen las sumas y diferencias de variables aleatorias exponenciales independientes?

Question

Si nos dan X ₁ , X ₂ , X ₃ , X ₄ todas las variables aleatorias exponenciales con la misma media y dijo que otra variable aleatoria T = X ₁ + X ₂ + X ₃ - X ₄ . ¿Podemos decir entonces que T tiene básicamente una distribución Gamma con parámetros (2,)? Sé que podemos sumar exponenciales y obtener una gamma; sin embargo, me pregunto si podemos aplicar la misma lógica para la resta. Creo que es posible considerando estas variables aleatorias como llegadas a un Proceso de Poisson, sin embargo quería preguntar.

Answer 1

Porque la función característica de una distribución exponencial es

$$\phi(t) = \frac{1}{1 - it},$$

cuando $X_1,\ldots, X_4$ son variables independientes distribuidas exponencialmente (todas con la misma tasa, que sin pérdida de generalidad puede tomarse como unidad), la función característica de $X = X_1+X_2+X_3-X_4$ es

$$\phi_X(t) = \phi(t)^3\phi(-t) = \frac{1}{(1-it)^3(1+it)}.$$

Utilizando fracciones parciales, esto se descompone e invierte fácilmente (por inspección o cálculo directo de la transformada de Fourier inversa) para obtener la densidad

$$f_X(x) = \frac{1}{8} \left(e^{x} \mathcal{I}(x\le 0) + e^{-x}(1 + 2x + 2x^2)\mathcal{I}(x\gt 0)\right).$$

Esto es reconocible como una mezcla de una distribución Laplace, una Gamma $(2)$ y una distribución Gamma $(3)$ distribución.

Como comprobación, en R Generé un millón de independientes $(X_1,X_2,X_3,X_4)$ tuplas y trazamos este histograma de $X_1+X_2+X_3-X_4$ para compararlo con el gráfico de $f_X,$ como se muestra.

hist(c(1,1,1,-1) %*% matrix(rexp(1e6*4), 4), freq=FALSE, breaks=200, xlab="x", main="")
curve(1/8 * (dgamma(-x, 1) + dgamma(x, 1) + 2*dgamma(x, 2) + 4*dgamma(x, 3)))

El acuerdo es excelente.

---

Generalización

Como la pregunta se refiere a "sumas y diferencias" de forma bastante general, he pensado que podría ser interesante obtener una respuesta general. Para ello, observe que la función característica (cf) de la suma de $k$ iid menos la suma de $j$ iid las variables exponenciales serán

$$\phi_{k,j}(t) = (1 - it)^{-k}(1 + it)^{-j}.$$

El objetivo es expresarlo, si es posible, como una mezcla (combinación lineal) de Gamma cfs. Recordemos que una distribución Gamma con parámetro de forma $\alpha$ tiene $(1-it)^{-\alpha}$ para su cf. Se ve inmediatamente que $\phi_{k,j}$ es la cf de la diferencia de una Gamma $(k)$ y Gamma (independiente) $(j)$ variable. Pero eso no nos hace avanzar hacia una solución.

Este problema resulta ser puramente algebraico. Al escribir

$$a = \frac{2}{1-it},\quad b = \frac{2}{1+it},$$

observe que $1/a + 1/b = 1,$ de donde $a + b = ab.$ Nuestro objetivo, entonces, puede ser alcanzado si podemos expresar $a^k b^j$ como una combinación lineal de potencias puras (positivas) de $a$ y $b.$ (Formalmente, los cálculos se realizan en el anillo de cociente de polinomios $\mathbb{Q}(a,b)/(a+b-ab).$ No serán necesarias consideraciones analíticas de convergencia, etc.).

La idea es encontrar $a^kb^j$ recursivamente.

• El caso base es $a^k=a^k:$ ya tenemos un poder puro de $a.$

• El siguiente caso es $$a^kb = a^{k-1}(ab) = a^{k-1}(a+b) = a^k + a^{k-1}b.$$ Por lo tanto, si $k\gt 1,$ repetir este proceso, encontrando finalmente $$a^kb = a^k + a^{k-1} + \cdots + a^2 + a + b.$$

• El siguiente paso nos muestra lo que ocurrirá en general, utilizando el resultado anterior para expandir cada uno de los términos que aparecen en la multiplicación final por $b$ : $$\begin{aligned} a^k b^2 &= (a^k b)b = \left(a^k + a^{k-1} + \cdots + a^2 + a + b\right)b\\ &= a^kb + a^{k-1}b + \cdots + a^2b + ab + b^2\\ &= a^k + a^{k-1} + \cdots + a^2 + a + b\\ &+ \quad a^{k-1} + \cdots + a^2 + a + b\\ &\cdots\\ &+ \quad\quad\quad (a+b) + b^2\\ &= \left(a^k + 2a^{k-1} + \cdots + ka\right) + \left(k b + b^2\right). \end{aligned}$$

Procediendo de este modo, llegamos a una fórmula general (que se demuestra fácilmente de forma inductiva):

$$a^k b^j = \sum_{i=1}^k \binom{k+j-1-i}{j-1} a^{i}\ +\ \sum_{i=1}^j\binom{k+j-1-i}{k-1} b^{i}.$$

(Es una relación bastante bonita en sí misma...)

Al dividir ambos lados por $2^{k+j}$ y reexpresar $a$ y $b$ en términos de $t,$ obtenemos la descomposición deseada, de la que podemos leer directamente los parámetros de forma Gamma y los pesos de la mezcla:

$$\phi_{k,j}(t)=\sum_{i=1}^k \binom{k+j-1-i}{j-1}2^{i-j-k} \phi(t)^i\ +\ \sum_{i=1}^j\binom{k+j-1-i}{k-1} 2^{i-j-k} \phi(-t)^i.$$

Por ejemplo, la combinación de la pregunta corresponde a una suma de $k=3$ exponenciales menos $j=1$ exponencial, para lo cual

$$\begin{aligned} \phi_{3,1}(t) &= \binom{3-1}{0}2^{1-4}\phi(t)+\binom{3-2}{0}2^{2-4}\phi(t)^2+\binom{3-3}{0}2^{3-4}\phi(t)^3\\&\quad+ \binom{1-1}{0}2^{1-4}\phi(-t)\\ &= \frac{1}{8}\phi(t)+\frac{1}{4}\phi(t)^2+\frac{1}{2}\phi(t)^3+\frac{1}{8}\phi(-t), \end{aligned}$$

que nos dice que la distribución es una mezcla de $1/8$ exponencial, $1/4$ Gamma $(2),$ $1/2$ Gamma $(3),$ y $1/8$ distribuciones exponenciales invertidas.

Aquí hay un R implementación modelada en el código anterior. También incluye un parámetro de tasa común para todos los exponenciales. (No juega ningún papel en los cálculos porque simplemente establece una unidad de medida común para todas las distribuciones).

#
# Density of the sum of `k` and difference of `j` *iid* exponentials.
# (Equivalently, the difference of a Gamma(k) and Gamma(j) distribution.)
#
f <- Vectorize(function(x, k, j, ...) {
  g <- function(x, k, j) {
    i <- seq_len(k)
    sum(choose(k+j-1-i, j-1) * 2^(i-j-k) * dgamma(x, i, ...))
  }
  g(x, k, j) + g(-x, j, k)
}, "x")
#
# Simulated data.
#
k <- 9
j <- 4
rate <- 10
hist(c(rep(1,k), rep(-1,j)) %*% matrix(rexp(1e6 * (k+j), rate), k+j), 
     freq=FALSE, breaks=200, xlab="x", 
     main=bquote(paste(k==.(k), ", ", j==.(j), ", and rate"==.(rate))))
#
# Check by overplotting the density function.
#
curve(f(x, k, j, rate=rate), add=TRUE, lwd=2)

Los componentes positivos tendrán distribuciones Gamma con parámetros de forma $1, 2, \ldots, k$ mientras que los componentes negativos, al invertirse, tendrán distribuciones Gamma con parámetros de forma $1,2,\ldots, j.$ Como ilustración, el siguiente histograma, producido por el código anterior, también muestra los componentes, convenientemente ponderados y sombreados en color para diferenciarlos. Su suma se representa en negro, trazando el histograma.