Inversión numérica

01 · Punto de partida

Cuando sabés evaluar F pero no despejarla

La normal, la gamma general, la beta y cualquier densidad definida por datos o por integral no tienen Q elemental. Pero si podés evaluar F(x) —con erf, integración numérica o suma acumulada—, podés resolver F(x) = u por métodos numéricos. La monotonía de F convierte un problema “difícil” en uno con convergencia garantizada.

Este tema es el plan B universal de la inversa: más lento por muestra que una fórmula, pero aplicable a casi todo y exacto a tolerancia pedida.

• ¿Cómo encerrar la raíz entre lo y hi que sí bracketan u?
• ¿Cuántas iteraciones de bisección dan 6 decimales?
• ¿Bisección, Newton o regula falsi con F monótona?
• ¿Cómo pasar de 40 evaluaciones por muestra a 1 con tabla?

02 · Definición

Bisección sobre una monótona

Dado u ∈ (0,1): partir de [lo,hi] con F(lo) ≤ u ≤ F(hi). Repetir: mid = (lo+hi)/2; si F(mid) < u, lo = mid; si no, hi = mid. El ancho se parte a la mitad por iteración:

Newton (x ← x − (F(x)−u)/f(x)) converge cuadráticamente cerca de la raíz pero falla en colas planas donde f ≈ 0. Regla práctica: bisección para robustez, Newton solo con buen arranque y f evaluable.

03 · Anatomía

Del bracket a la tabla reutilizable

1. Dominio y tolerancia

Fijar [lo,hi] que cubra u ∈ [tiny, 1−tiny] y tolerancia tol (ej. 1e−10 en x o en F). Truncar colas extremas es una decisión de modelo, no un bug.

2. Bisección robusta

Iterar hasta hi − lo < tol o iter máx. Devolver (lo+hi)/2. Siempre converge; el error está acotado por el ancho actual.

3. Tabla de cuantiles

Elegir malla uⱼ (uniforme en el centro, densa en colas), resolver xⱼ = Q(uⱼ) una vez, guardar. Muestrear: buscar j con uⱼ ≤ u < uⱼ₊₁ e interpolar lineal.

4. Validación

Comparar tabla vs. bisección directa en 5–10 u de prueba y validar histograma vs. f. El error de interpolación se ve en colas, no en la media.

04 · Dónde aplica

El comodín antes de métodos dedicados

1. 1
   Normal(0,1).

   F vía erf: Φ(x) = ½(1+erf(x/√2)). Bisección didáctica; en producción Box-Muller o algoritmo dedicado (Temas 23–25).

2. 2
   Gamma / beta generales.

   F por integración o función especial; tabla de cuantiles reutilizable para toda la simulación.

3. 3
   Densidades a medida.

   f definida por polinomio/trozos/datos suavizados: integrar numéricamente a F y luego invertir. (Temas 34–35)

4. 4
   Mezclas.

   Elegir componente (Tema 11) y luego invertir la componente, analítica o numéricamente según el caso.

05 · Representación en Python

Bisección y tabla para la normal

Sin librerías externas: math.erf da F, bisección da Q, y una tabla acelera el muestreo masivo.

import random
import math

def Phi(x):
    return 0.5 * (1 + math.erf(x / math.sqrt(2)))

def Q_bisect(u, lo=-10.0, hi=10.0, it=40):
    for _ in range(it):
        mid = (lo + hi) / 2
        if Phi(mid) < u:
            lo = mid
        else:
            hi = mid
    return (lo + hi) / 2

print("Q(0.975) =", round(Q_bisect(0.975), 3), "(teórico 1.960)")

# Tabla m=2000 una vez, luego interpolación lineal.
m = 2000
ug = [ (j + 0.5) / m for j in range(m) ]
xg = [Q_bisect(u, it=35) for u in ug]

def Q_tabla(u):
    u = min(max(u, ug[0]), ug[-1])
    # Búsqueda binaria sobre malla uniforme en u.
    import bisect
    j = bisect.bisect_left(ug, u)
    if j <= 0: return xg[0]
    if j >= m: return xg[-1]
    u0, u1 = ug[j-1], ug[j]
    x0, x1 = xg[j-1], xg[j]
    t = (u - u0) / (u1 - u0)
    return x0 + t * (x1 - x0)

random.seed(50)
N = 4000
muestra = [Q_tabla(random.random()) for _ in range(N)]
print("media:", round(sum(muestra) / N, 3), " P(X<=1):",
      round(sum(1 for x in muestra if x <= 1) / N, 3), "(teórica 0.841)")

Q(0,975) debe dar ≈ 1,96. Si tu tabla da 1,92, aumentá m o refiná colas: el error vive en los extremos.

06 · Exploración

Laboratorio: ver converger la bisección

Fijá u y aumentá las iteraciones. El intervalo [lo,hi] (banda) se achica a la mitad por paso hasta abrazar el cuantil normal. La curva es Φ(x).

Eje x ∈ [−4,4], eje y = Φ. La línea horizontal es u; la banda vertical es el bracket actual; el punto es el cuantil.

Preguntas para explorar

1. Con u = 0,975, ¿a qué valor converge Q al subir a 20 iteraciones? (Pista: 1,96.)
2. Con 1 iteración, ¿qué ancho tiene el intervalo si el bracket inicial es [−10,10]?
3. Con u = 0,99 y 4 iteraciones, ¿el error es mayor o menor que con u = 0,5? ¿Por qué la pendiente de Φ lo explica?

07 · Comprensión

Confusiones frecuentes

08 · Práctica guiada

Ejercicios con Python

Ejercicio 1: cuantil normal a mano

Implementá Q_bisect con Φ vía erf y calculá Q(0,5), Q(0,8413) y Q(0,975). Deben dar ≈ 0, 1 y 1,96.

import math

def Phi(x):
    return 0.5 * (1 + math.erf(x / math.sqrt(2)))

def Q(u, it=40):
    lo, hi = -10.0, 10.0
    for _ in range(it):
        mid = (lo + hi) / 2
        if Phi(mid) < u:
            lo = mid
        else:
            hi = mid
    return (lo + hi) / 2

for u in [0.5, 0.8413, 0.975]:
    print(u, round(Q(u), 3))

Ejercicio 2: costo de la precisión

Medí el ancho tras 10, 20 y 30 iteraciones con bracket 20. Verificá la ley ancho = 20/2ⁿ.

for n in [10, 20, 30]:
    print(n, 20 / 2 ** n)

Ejercicio 3: tabla que sí alcanza

Construí tabla m=500 y m=4000 para la normal y compará Q_tabla(0,999) contra Q_bisect(0,999). ¿Cuál se acerca más?

import math
import bisect

def Phi(x):
    return 0.5 * (1 + math.erf(x / math.sqrt(2)))

def Qb(u, it=45):
    lo, hi = -10.0, 10.0
    for _ in range(it):
        mid = (lo + hi) / 2
        if Phi(mid) < u:
            lo = mid
        else:
            hi = mid
    return (lo + hi) / 2

def make_table(m):
    ug = [(j + 0.5) / m for j in range(m)]
    return ug, [Qb(u, 35) for u in ug]

def Qtab(u, ug, xg):
    import bisect
    j = bisect.bisect_left(ug, u)
    j = min(max(j, 1), len(ug) - 1)
    u0, u1, x0, x1 = ug[j-1], ug[j], xg[j-1], xg[j]
    return x0 + (u - u0) / (u1 - u0) * (x1 - x0)

ref = Qb(0.999)
for m in [500, 4000]:
    ug, xg = make_table(m)
    print(m, round(Qtab(0.999, ug, xg), 4), "ref", round(ref, 4))

09 · Síntesis

Ideas para recordar

En el próximo tema veremos la alternativa cuando ni siquiera conviene invertir: aceptación y rechazo.