Método Box-Muller

01 · Punto de partida

De la planicie al giro gaussiano

La normal no tiene inversa cerrada, pero dos uniformes sí alcanzan si se mira en polares: un radio exponencial y un ángulo uniforme, girados con coseno y seno, dan dos N(0,1) independientes de un tiro.

Es el método exacto por defecto del Tema 23 y la referencia contra la que se mide el polar del Tema 25: misma exactitud, distinto costo trigonométrico.

• ¿Por qué R = √(−2 ln U₁) y Θ = 2πU₂?
• ¿Por qué Z₀ y Z₁ salen independientes?
• ¿Qué le cuesta al método cada par (log, raíz, senos)?
• ¿Qué pasa si U₁ = 0 o si se reutiliza mal el par?

02 · Definición

Radio exponencial y ángulo uniforme

R² = −2 ln U₁ ~ Exp(1/2): el radio al cuadrado es exponencial. Θ ~ U(0,2π) reparte la dirección pareja. El giro (cos, sin) proyecta el punto en dos coordenadas gaussianas independientes.

03 · Anatomía

Costo, bordes y sobrante

Costo por par

Un log, una raíz, un cos y un sin para dos normales (≈ 0,5 de cada uno por Z). El seno/coseno domina: por eso existe el polar (Tema 25).

Borde U₁ = 0

ln 0 = −∞ y R explota. Rechazar U₁ = 0 exacto o clampear a 1e−12; U₂ = 0 es inofensivo (ángulo 0).

Sobrante (“spare”)

Generar de a pares y devolver una por vez: la segunda queda cacheada. Sin caché se desperdicia la mitad del trabajo.

Independencia frágil

Z₀ y Z₁ del mismo par son independientes entre sí; mezclar U₁,U₂ entre pares o reutilizar Θ rompe la independencia aunque cada marginal parezca normal.

04 · Dónde gira Box-Muller

Del par a la librería

1. 1
   Ruido gaussiano.

   Sensores, errores y simulación física con N(μ,σ) vía μ + σZ.

2. 2
   Pares correlacionados.

   Combinar Z₀,Z₁ independientes para lograr correlación objetivo. (Temas 28–29)

3. 3
   Referencia exacta.

   Patrón para validar el polar, ziggurat y aproximaciones rápidas. (Temas 25, 39)

4. 4
   Colas y riesgo.

   VaR y eventos 3σ+ donde la suma-12 no llega.

05 · Representación en Python

Par exacto con sobrante

Implementación con caché: un llamado produce dos normales y guarda la segunda. Auditoría de par completo.

import random
import math

random.seed(200)
N = 8000
_spare = []

def box_muller():
    if _spare:
        return _spare.pop()
    while True:
        u1 = random.random()
        if u1 > 1e-12:
            break
    u2 = random.random()
    r = math.sqrt(-2.0 * math.log(u1))
    t = 2.0 * math.pi * u2
    _spare.append(r * math.sin(t))
    return r * math.cos(t)

m = [box_muller() for _ in range(N)]
media = sum(m) / N
var = sum((x - media) ** 2 for x in m) / N
pares = [(m[i], m[i + 1]) for i in range(0, N - 1, 2)]
pm0 = sum(a for a, _ in pares) / len(pares)
pm1 = sum(b for _, b in pares) / len(pares)
cov = sum((a - pm0) * (b - pm1) for a, b in pares) / len(pares)
print("media:", round(media, 3), "var:", round(var, 3))
print("corr aprox:", round(cov, 3), "vs 0.0")
print("banda |Z|<=1:", round(sum(1 for x in m if abs(x) <= 1) / N, 3), "vs 0.683")

Media ≈ 0, var ≈ 1, corr ≈ 0 y banda ≈ 0,683. Si la corr se despega, revisá el caché del par.

06 · Exploración

Laboratorio: nube redonda e incorrelada

Cada punto es un par (Z₀,Z₁) Box-Muller: la nube debe ser redonda (simetría radial) y la correlación ≈ 0. Los anillos marcan 1σ, 2σ y 3σ. Subí N para estabilizar las estadísticas.

Dispersión (Z₀,Z₁) en [−4,4]² con anillos 1σ–3σ. Probá N = 500 (nube rala) y N = 12000 (redonda y densa).

Preguntas para explorar

1. Con N = 3000, ¿medias ≈ 0 y corr ≈ 0? Leelas en el panel.
2. ¿La nube se alarga en diagonal o es redonda? ¿Qué dice eso de la independencia?
3. Bajá a N = 500. ¿La corr se aleja de 0? ¿Es sesgo o ruido?

07 · Comprensión

Confusiones frecuentes

08 · Práctica guiada

Ejercicios con Python

Ejercicio 1: par auditado

N = 10000 (semilla 201): verificá medias ≈ 0, varianzas ≈ 1 y |corr| < 0,03.

import random
import math

random.seed(201)
N = 10000
pares = []
for _ in range(N // 2):
    u1 = max(random.random(), 1e-12)
    u2 = random.random()
    r = math.sqrt(-2 * math.log(u1))
    t = 2 * math.pi * u2
    pares.append((r * math.cos(t), r * math.sin(t)))
z0 = [a for a, _ in pares]
z1 = [b for _, b in pares]
m0, m1 = sum(z0)/len(z0), sum(z1)/len(z1)
cov = sum((a-m0)*(b-m1) for a, b in pares) / len(pares)
print(round(m0, 3), round(m1, 3), round(cov, 3))

Ejercicio 2: de Z a N(μ,σ)

Con μ = 5, σ = 3 y N = 8000, transformá el par y verificá media ≈ 5 y banda [2,8] ≈ 0,683.

import random
import math

random.seed(202)
mu, sigma, N = 5.0, 3.0, 8000
m = []
for _ in range(N // 2):
    u1 = max(random.random(), 1e-12)
    u2 = random.random()
    r = math.sqrt(-2 * math.log(u1))
    t = 2 * math.pi * u2
    m += [mu + sigma * r * math.cos(t), mu + sigma * r * math.sin(t)]
media = sum(m) / N
print(round(media, 3), "vs 5.0")
print(round(sum(1 for x in m if abs(x-mu) <= sigma) / N, 3), "vs 0.683")

Ejercicio 3: el costo del seno

Generá 20000 normales por Box-Muller y por suma-12 y compará uniformes gastados: 1 por Z contra 12.

import random
import math
import time

random.seed(203)
N = 20000
t0 = time.perf_counter()
for _ in range(N // 2):
    u1 = max(random.random(), 1e-12)
    u2 = random.random()
    r = math.sqrt(-2 * math.log(u1))
    t = 2 * math.pi * u2
    _ = (r * math.cos(t), r * math.sin(t))
t1 = time.perf_counter()
random.seed(203)
for _ in range(N):
    _ = sum(random.random() for _ in range(12)) - 6.0
t2 = time.perf_counter()
print("U por Z: Box-Muller 1.0 vs suma-12 12.0")
print("tiempos:", round(t1-t0, 3), "vs", round(t2-t1, 3))

09 · Síntesis

Ideas para recordar

En el próximo tema evitamos la trigonometría: el método polar con rechazo en el disco.