Método del alias

01 · Punto de partida

Cuando la discreta se muestrea millones de veces

Buscar en acumuladas cuesta O(n) lineal u O(log n) binaria por muestra. Con n = 10 000 categorías y N = 10⁸, eso es prohibitivo. El alias paga O(n) una vez y luego muestrea en O(1): un entero uniforme + una comparación.

La idea visual: rebanar toda la masa en n columnas de altura exacta 1/n, cada una con a lo sumo dos valores (el propio y su alias). Muestrear es elegir columna y decidir entre sus dos habitantes.

• ¿Cómo se parte una masa arbitraria en n columnas de 1/n?
• ¿Qué guardan las tablas Prob y Alias?
• ¿Por qué cada muestra cuesta solo dos uniformes discretos?
• ¿Cuándo conviene pagar el preproceso?

02 · Definición

Columnas de 1/n con a lo sumo dos valores

Escalar pᵢ por n: qᵢ = n·pᵢ promedia 1. Los qᵢ < 1 (chicos) ceden su faltante a qⱼ > 1 (grandes). Al final cada columna i tiene Prob[i] = fracción propia y Alias[i] = donante del resto.

03 · Anatomía

Construcción de Walker en O(n)

Pilas chicas y grandes

qᵢ < 1 van a “small”, qᵢ ≥ 1 a “large”. Sacar un chico s y un grande g: Prob[s] = qₛ, Alias[s] = g; descontar (1 − qₛ) de q_g y reclasificar g.

Invariante

La masa total de pilas + columnas cerradas siempre suma n. Al vaciar pilas, cada columna mide 1/n exacto (salvo redondeo).

Memoria 2n

Dos arreglos de n (Prob flotante + Alias entero). Para n = 10⁶ son ~12 MB: aceptable frente al ahorro en N.

Muestreo sin sesgo

P(X = k) = Σᵢ (1/n)·P(devolver k|columna i) = pₖ por construcción. El redondeo se absorbe forzando Prob = 1 en restos.

04 · Cuándo pagar el preproceso

Regla N ≫ n

1. 1
   Vocabularios y grillas.

   Miles de categorías fijas (lenguaje, grillas espaciales) con miles de millones de sorteos.

2. 2
   Motores de simulación.

   Demanda/discretas empíricas fijas que se muestrean en cada réplica y paso.

3. 3
   No usar si…

   p cambia por contexto, N ≈ n, o n ≤ 10 (lineal ya es O(1) práctico).

4. 4
   Detalle en Tema 14.

   Construcción paso a paso con pilas y verificación de columnas.

05 · Representación en Python

Walker + sampler en 25 líneas

Construcción clásica y muestreo con dos uniformes. La verificación compara frecuencias contra p.

import random
from collections import deque

random.seed(100)
probs = [0.10, 0.20, 0.30, 0.25, 0.10, 0.05]
n = len(probs)

def construir_alias(probs):
    n = len(probs)
    q = [p * n for p in probs]
    prob = [0.0] * n
    alias = [0] * n
    small = deque([i for i, v in enumerate(q) if v < 1 - 1e-12])
    large = deque([i for i, v in enumerate(q) if v >= 1 - 1e-12])
    while small and large:
        s = small.popleft()
        g = large.popleft()
        prob[s] = q[s]
        alias[s] = g
        q[g] = q[g] - (1 - q[s])
        if q[g] < 1 - 1e-12:
            small.append(g)
        else:
            large.append(g)
    for i in list(small) + list(large):
        prob[i] = 1.0
        alias[i] = i
    return prob, alias

def alias_sample(prob, alias):
    n = len(prob)
    i = int(random.random() * n)
    return i if random.random() < prob[i] else alias[i]

prob, alias = construir_alias(probs)
print("Prob:", [round(v, 3) for v in prob])
print("Alias:", alias)
N = 20000
m = [alias_sample(prob, alias) + 1 for _ in range(N)]
for v in range(1, 7):
    print(v, round(m.count(v) / N, 3))

Frecuencias ≈ tabla original. Con la justa, Prob = [1,1,1,1,1,1]: cada columna es pura.

06 · Exploración

Laboratorio: alias que calza en O(1)

Elegí distribución y N. El muestreo usa alias real en el navegador (dos uniformes por muestra). Las tablas Prob/Alias se muestran debajo; el histograma debe pegarse a los rombos teóricos.

Seis barras vs rombos teóricos. Cambiá a justa para ver Prob todo en 1.

Preguntas para explorar

1. Con cargada, ¿qué columnas tienen alias? Leé la línea Prob/Alias.
2. Con justa, ¿por qué Prob = 1 en todas? ¿Qué significa “columna pura”?
3. Si N = 500, ¿el error máx. es del método o del ruido? ¿Cómo lo distinguís?

07 · Comprensión

Confusiones frecuentes

08 · Práctica guiada

Ejercicios con Python

Ejercicio 1: verificar columnas

Con la tabla del dado cargado, comprobá que cada pᵢ se recupera desde (Prob, Alias) con la fórmula de la buena práctica.

probs = [0.10, 0.20, 0.30, 0.25, 0.10, 0.05]
prob = [0.6, 1.0, 1.0, 1.0, 0.6, 0.3]
alias = [2, 1, 2, 3, 3, 3]  # 0-based: [2,1,2,3,3,3]-1
n = 6
for k in range(n):
    pk = prob[k] / n + sum((1 - prob[j]) / n
                           for j in range(n) if alias[j] == k and prob[j] < 1)
    print(k + 1, round(pk, 3), "vs", probs[k])

Ejercicio 2: justa trivial

Construí alias para 6 × 1/6 y verificá Prob todo 1. Muestreá 6000 y compará contra randint.

import random

random.seed(101)
prob = [1.0] * 6
alias = list(range(6))
m = [(int(random.random() * 6) if random.random() < prob[int(random.random()*6)] else 0) for _ in range(6)]
# Muestreo correcto simple:
random.seed(101)
m = []
for _ in range(6000):
    i = int(random.random() * 6)
    m.append(i if random.random() < prob[i] else alias[i])
print([round(m.count(v) / len(m), 3) for v in range(6)])

Ejercicio 3: medir la ganancia

Compará tiempo de 200 000 muestras por bisect vs alias con n = 100 categorías. Alias debe ganar claro.

import random
import bisect
import time
from collections import deque

random.seed(102)
n = 100
probs = [random.random() for _ in range(n)]
s = sum(probs)
probs = [p / s for p in probs]
acum = []
a = 0.0
for p in probs:
    a += p
    acum.append(a)
acum[-1] = 1.0

def construir(probs):
    n = len(probs)
    q = [p * n for p in probs]
    prob = [0.0] * n
    alias = [0] * n
    small = deque([i for i, v in enumerate(q) if v < 1])
    large = deque([i for i, v in enumerate(q) if v >= 1])
    while small and large:
        st = small.popleft()
        lg = large.popleft()
        prob[st] = q[st]
        alias[st] = lg
        q[lg] = q[lg] - (1 - q[st])
        (small if q[lg] < 1 else large).append(lg)
    for i in list(small) + list(large):
        prob[i] = 1.0
    return prob, alias

prob, alias = construir(probs)
N = 200000
t0 = time.perf_counter()
for _ in range(N):
    _ = bisect.bisect_left(acum, random.random())
t1 = time.perf_counter()
for _ in range(N):
    i = int(random.random() * n)
    _ = i if random.random() < prob[i] else alias[i]
t2 = time.perf_counter()
print(f"bisect: {t1-t0:.3f}s  alias: {t2-t1:.3f}s")

09 · Síntesis

Ideas para recordar

En el próximo tema construiremos esas tablas a mano: construcción de tablas de alias.