Generación gamma

01 · Punto de partida

Cuánto tarda la k-ésima llegada

¿Cuánto hasta el tercer cliente, la quinta falla, la décima llamada? Si cada espera es Exp(λ), la espera acumulada es Gamma(α,λ): α cuenta cuántas exponenciales se suman y λ mantiene la tasa.

Es la generalización continua de la Erlang y la madre de la chi-cuadrado: con α entero basta sumar, con α real hace falta un método general como Marsaglia-Tsang.

• ¿Por qué sumar α exponenciales da Gamma(α,λ)?
• ¿Qué forma toma con α < 1, α = 1 y α grande?
• ¿Cómo se escala la tasa con Y/λ?
• ¿Cuándo usar suma, Marsaglia-Tsang o librería?

02 · Definición

Forma α, tasa λ y casos famosos

Soporte [0,∞). E[X] = α/λ, Var(X) = α/λ². Con α = 1 es Exp(λ); con α entero es Erlang; con λ = 1/2 y α = ν/2 es chi-cuadrado con ν grados.

03 · Anatomía

Suma, Marsaglia-Tsang y escala

Suma de exponenciales O(α)

Con α = k entero: −Σ ln(Uᵢ)/λ. Gasta k uniformes por muestra. Perfecta si k ≤ 10 y didáctica siempre.

Marsaglia-Tsang para α ≥ 1

Propuesta normal con d = α−1/3 y c = 1/√(9d): v = (1+cZ)³ aceptada con test logarítmico. O(1) por muestra, exacta y con colas reales.

Boost para α < 1

Generar G ~ Gamma(α+1,1) y devolver G·U^(1/α). Reduce el caso difícil al caso fácil con un uniforme extra.

Escala de tasa

Todo se genera con λ = 1 y se divide por λ al final. Validar media y varianza después de escalar, no antes.

04 · Dónde suma la gamma

De la espera al modelo

1. 1
   Fiabilidad y colas.

   Tiempo hasta la k-ésima falla o llegada: Erlang con α entera.

2. 2
   Chi-cuadrado.

   Gamma(ν/2,1/2): base de tests e intervalos de varianza. (Tema 38)

3. 3
   Mezclas y Bayes.

   Prior conjugada de Poisson y exponencial: tasas aleatorias con forma α.

4. 4
   Puente a normal.

   Con α grande, (X−α/λ)/√(α/λ²) ⇒ N(0,1): otra campana por suma.

05 · Representación en Python

Suma entera frente a gammavariate

Mismos α y λ: la suma didáctica contra random.gammavariate (general). Ambas deben calzar en media y varianza.

import random
import math

random.seed(220)
alpha, lam, N = 3, 1.5, 8000
inv_lam = 1.0 / lam

def gamma_erlang(k, lam):
    return -sum(math.log(1 - random.random()) for _ in range(k)) / lam

m1 = [gamma_erlang(alpha, lam) for _ in range(N)]
random.seed(220)
m2 = [random.gammavariate(alpha, inv_lam) for _ in range(N)]

def resumen(m):
    media = sum(m) / N
    var = sum((x - media) ** 2 for x in m) / N
    return round(media, 3), round(var, 3)

print("erlang:", resumen(m1), "vs (2.0, 1.333)")
print("librería:", resumen(m2), "vs (2.0, 1.333)")

Media ≈ 2,0 y var ≈ 1,333 con α = 3 y λ = 1,5. Ojo: gammavariate pide escala θ = 1/λ, no tasa.

06 · Exploración

Laboratorio: α esculpe, λ escala

Movete α con λ y N fijos. El histograma usa 20 bines hasta el cuantil 99,5 % (cola aparte); los rombos son la probabilidad teórica por bin. Con α = 1 ves exponencial; con α = 6, campana asimétrica.

20 bines hasta el cuantil 99,5 % más cola. Probá α = 0,5 (pico en 0) y α = 8 (campana).

Preguntas para explorar

1. Con α = 3 y λ = 1, ¿media y varianza teóricas?
2. Bajá α a 1. ¿Recuperás la exponencial del Tema 22? ¿Media y var?
3. Duplicá λ a 2 con α = 3. ¿La media se parte a la mitad? ¿Y la varianza?

07 · Comprensión

Confusiones frecuentes

08 · Práctica guiada

Ejercicios con Python

Ejercicio 1: tercera llegada

α = 3, λ = 2, N = 8000 (semilla 221): estimá media (≈ 1,5) y varianza (≈ 0,75).

import random
import math

random.seed(221)
alpha, lam, N = 3, 2.0, 8000
m = [-sum(math.log(1 - random.random()) for _ in range(alpha)) / lam
     for _ in range(N)]
media = sum(m) / N
print(round(media, 3), "vs 1.5")
print(round(sum((x-media)**2 for x in m)/N, 3), "vs 0.75")

Ejercicio 2: chi-cuadrado como gamma

Generá χ² con ν = 4 como Gamma(2; 0,5) con N = 8000 y verificá media ≈ 4 y var ≈ 8.

import random

random.seed(222)
N = 8000
m = [random.gammavariate(2, 2.0) for _ in range(N)]
media = sum(m) / N
print(round(media, 3), "vs 4.0")
print(round(sum((x-media)**2 for x in m)/N, 3), "vs 8.0")

Ejercicio 3: escala de tasa

Generá 3000 Gamma(4,1) y dividí por λ = 2: debe dar Gamma(4,2) con media ≈ 2.

import random

random.seed(223)
alpha, lam, N = 4, 2.0, 3000
base = [random.gammavariate(alpha, 1.0) for _ in range(N)]
m = [x / lam for x in base]
media = sum(m) / N
print(round(media, 3), "vs 2.0")
print(round(sum((x-media)**2 for x in m)/N, 3), "vs 1.0")

09 · Síntesis

Ideas para recordar

En el próximo tema abrimos el catálogo: beta, Weibull, lognormal y triangular.