Método polar

01 · Punto de partida

El giro sin trigonometría

Box-Muller es exacto pero paga un coseno y un seno por par. El método polar (Marsaglia) logra lo mismo sin trigonometría: tira puntos en el cuadrado [−1,1]², se queda con los que caen en el disco S < 1 y los estira con un solo logaritmo.

Es el mismo par gaussiano del Tema 24 por otro camino: un 21 % de los intentos se rechaza, pero cada intento aceptado evita la parte más cara del cálculo.

• ¿Por qué V₁,V₂ ~ U(−1,1) con S = V₁²+V₂² < 1?
• ¿De dónde sale el factor √(−2 ln S / S)?
• ¿Qué tasa de aceptación esperar y cómo auditarla?
• ¿Cuándo conviene polar sobre Box-Muller?

02 · Definición

Disco uniforme y estiramiento radial

Condicionar (V₁,V₂) a S < 1 da un punto uniforme en el disco: su dirección es el ángulo Θ parejo y su radio al cuadrado es uniforme, el sustituto exacto de U₂ y U₁ de Box-Muller. El factor f convierte ese radio en el R gaussiano.

03 · Anatomía

Aceptar, estirar y cachear

Bucle de aceptación

Tirar V₁ = 2U−1, V₂ = 2U−1 hasta 0 < S < 1. S = 0 se rechaza (ln 0 indefinido); S ≥ 1 cae fuera del disco. Media de intentos: 4/π ≈ 1,27.

Factor único

f = sqrt(-2 * log(s) / s) una sola vez por par aceptado; Z₀ y Z₁ comparten el factor. No recalcular el log por coordenada.

Sobrante como Box-Muller

Devolver Z₀ y cachear Z₁. Sin caché se tira la mitad del trabajo aceptado, que ya costó rechazos.

Costo real

Más uniformes por par (2,55 vs 2) pero sin senos: en CPUs sin trigonometría rápida el polar gana; con sincos vectorizado la brecha se achica (Tema 39).

04 · Dónde rinde el disco

Del cuadrado a la librería

1. 1
   Mismo ruido que Box-Muller.

   N(μ,σ) vía μ + σZ, pares correlacionados y procesos gaussianos.

2. 2
   CPUs sin trigonometría.

   Microcontroladores y GPUs donde cos/sin son caros: el disco los reemplaza.

3. 3
   Clase de rechazo.

   Primer rechazo continuo que acepta sin deformar: antesala del Tema 8 aplicado.

4. 4
   Puente a ziggurat.

   La idea “capas + rechazo barato” escala al método más rápido de todos. (Tema 39)

05 · Representación en Python

Disco, factor y tasa π/4

Implementación con caché más auditoría triple: tasa de aceptación, par incorrelado y banda 68 %.

import random
import math

random.seed(210)
N = 8000
_spare = []
intentos = 0
aceptados = 0

def polar():
    global intentos, aceptados
    if _spare:
        return _spare.pop()
    while True:
        intentos += 1
        v1 = 2 * random.random() - 1
        v2 = 2 * random.random() - 1
        s = v1 * v1 + v2 * v2
        if 0 < s < 1:
            aceptados += 1
            f = math.sqrt(-2 * math.log(s) / s)
            _spare.append(v2 * f)
            return v1 * f

m = [polar() for _ in range(N)]
media = sum(m) / N
var = sum((x - media) ** 2 for x in m) / N
print("tasa:", round(aceptados / intentos, 3), "vs 0.785")
print("media:", round(media, 3), "var:", round(var, 3))
print("banda |Z|<=1:", round(sum(1 for x in m if abs(x) <= 1) / N, 3), "vs 0.683")

Tasa ≈ 0,785, media ≈ 0, var ≈ 1 y banda ≈ 0,683. Si la tasa da ≈ 1, el rechazo no está filtrando.

06 · Exploración

Laboratorio: el disco que acepta π/4

Cada punto es un intento (V₁,V₂): los celestes caen en el disco y se aceptan; los naranjas caen fuera y se rechazan. La fracción celeste debe rondar π/4 ≈ 0,785 y los Z resultantes, la banda 0,683.

Cuadrado [−1,1]² con disco inscripto. Probá N = 500 (tasa ruidosa) y N = 12000 (π/4 clavado).

Preguntas para explorar

1. Con N = 3000, ¿tasa ≈ 0,785? ¿Banda ≈ 0,683?
2. ¿Los puntos naranjas se concentran en las esquinas? ¿Por qué ahí?
3. Bajá a N = 500. ¿La tasa se aleja de π/4? ¿Sesgo o ruido?

07 · Comprensión

Confusiones frecuentes

08 · Práctica guiada

Ejercicios con Python

Ejercicio 1: tasa π/4 auditada

N = 10000 (semilla 211): verificá tasa ≈ 0,785, medias ≈ 0 y |corr| < 0,03.

import random
import math

random.seed(211)
N = 10000
intentos = aceptados = 0
pares = []
while len(pares) < N // 2:
    intentos += 1
    v1 = 2 * random.random() - 1
    v2 = 2 * random.random() - 1
    s = v1 * v1 + v2 * v2
    if 0 < s < 1:
        aceptados += 1
        f = math.sqrt(-2 * math.log(s) / s)
        pares.append((v1 * f, v2 * f))
z0 = [a for a, _ in pares]
m0 = sum(z0) / len(z0)
print("tasa:", round(aceptados / intentos, 3), "vs 0.785")
print("media Z0:", round(m0, 3))

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

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

import random
import math

random.seed(212)
mu, sigma, N = 5.0, 3.0, 8000
m = []
while len(m) < N:
    v1 = 2 * random.random() - 1
    v2 = 2 * random.random() - 1
    s = v1 * v1 + v2 * v2
    if 0 < s < 1:
        f = math.sqrt(-2 * math.log(s) / s)
        m += [mu + sigma * v1 * f, mu + sigma * v2 * f]
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: polar vs Box-Muller en uniformes

Contá uniformes por Z de cada método con N = 20000: Box-Muller 1,0 contra polar ≈ 1,27.

import random
import math

random.seed(213)
N = 20000
# Box-Muller: 2 U por 2 Z
print("Box-Muller U/Z:", 2 * (N // 2) * 2 / N / 2, "= 1.0")
# Polar: 2 U por intento, 1.27 intentos por par
intentos = aceptados = 0
while aceptados < N // 2:
    intentos += 1
    v1 = 2 * random.random() - 1
    v2 = 2 * random.random() - 1
    if 0 < v1*v1 + v2*v2 < 1:
        aceptados += 1
print("Polar U/Z:", round(2 * intentos / N, 3), "vs 1.27")

09 · Síntesis

Ideas para recordar

En el próximo tema sumamos exponenciales: la gamma con forma α.