Métodos Monte Carlo UQ: Hipercubo latino, Sobol y quasi-Monte Carlo

Cuando la salida de su simulación está impulsada por entradas inciertas, la forma en que extrae muestras de esas entradas determina si su cuantificación de incertidumbre es eficiente o derrochadora. El estándar Monte Carlo utiliza I.I.D. Sorteos aleatorios, lo que garantiza la convergencia pero a un ritmo dolorosamente lento. Métodos de muestreo de baja discrepancia: muestreo de hipercubo latino (LHS), secuencias SOBOL y la familia cuasi-Monte Carlo (QMC) más amplia: reestructurar el proceso de muestreo para que los puntos cubran el espacio de parámetros de manera más uniforme, reduciendo drásticamente la varianza de estimadores.

En este artículo vamos más allá de la breve mención en nuestra Guía de Simulaciones de Monte Carlo para Científicas. Esta es una inmersión profunda dedicada a la teoría, la implementación práctica y los marcos de decisión para LHS, secuencias SOBOL y métodos QMC, con código Python de trabajo y una comparación de cuándo gana cada enfoque.

Comida clave

• LHS estratifica cada dimensión de forma independiente, lo que garantiza que cada contenedor de histogramas marginales se llene exactamente una vez. Ofrece una convergencia más rápida que el sencillo Monte Carlo a bajo costo.
• Las secuencias SOBOL son secuencias de baja discrepancia que logran una cobertura casi uniforme en cualquier dimensión, pero requieren codificación para un uso práctico y funcionan mejor con recuentos de muestras de poderes de dos.
• La convergencia de QMC puede alcanzar O(N⁻¹) para integrar suaves, frente a O(N⁻¹ᐟ²) para Monte Carlo estándar, una caída exponencial en el tamaño de muestra requerido para la misma precisión.
• Las secuencias de Halton permiten recuentos de muestras arbitrarios y son fáciles de generar, pero exhiben un alias estructurado en dimensiones más altas. SOBOL supera a Halton cuando necesita una invariancia rotacional a través de la codificación.
• La codificación no es opcional para la producción SOBOL: las secuencias no codificadas inducen un sesgo sistemático en integrandos periódicos o discontinuos.

El problema de muestreo: por qué el aleatorio no siempre es mejor

El desafío central en la cuantificación de incertidumbre basada en el muestreo es la discrepancia del conjunto de muestras. La discrepancia mide cómo se distribuyen los puntos uniformemente en el espacio de parámetros. En el límite de muestras infinitas, cualquier discrepancia razonable llega a cero, pero la velocidad a la que hace importa enormemente.

Monte Carlo clásico

Standard Monte Carlo extrae muestras independientemente de la distribución objetivo. Por el teorema del límite central, el estimador converge con una varianza proporcional a 1/N. Eso significa que para reducir a la mitad el error estándar, debe cuadruple el recuento de muestras. Para una simulación computacionalmente costosa, un problema térmico 3D, un solucionador de flujo multifásico o una ejecución de dinámica molecular de ciencia de materiales, cuadruplicar el recuento de muestras puede significar días de reloj de pared.

La perspectiva de discrepancia

Las secuencias de baja discrepancia (LD) se construyen con precisión para minimizar la discrepancia de estrellas d*_n del conjunto de muestras. La definición es:

D*_N = sup_{A ⊆ [0,1]^s} |(1/N) Σ 1_A(x_i) − λ(A)|

donde λ(a) es la medida de Lebesgue (volumen) de A y el supremum se toma sobre todas las subcajas alineadas con eje. En inglés sencillo: las secuencias LDS minimizan la desviación máxima entre la probabilidad empírica y la verdadera sobre cualquier región rectangular.

El sobelev–dick-tractman desigualdad (el teorema de Koks-Meyers para LDS) establece que para una función f con derivados parciales mixtos hasta el orden s, el error de integración de una regla LDS está delimitado por:

|∫ f − (1/N) Σ f(x_i)| ≤ C_s · ‖f‖_S · D*_N

donde c_s es una constante que depende solo de la dimensión y ‖s‖_s es la norma sobolev. Este límite explica por qué los métodos de QMC brillan para funciones suaves: el error decae mucho más rápido que la tasa de Monte Carlo cuando el integrando tiene derivados mixtos.

La maldición de la dimensionalidad

Todos estos métodos se enfrentan a la maldición de la dimensionalidad. Para LHS, el beneficio de estratificación se compone multiplicativamente entre dimensiones. Para Sobol y Halton, la estructura que da baja discrepancia en las dimensiones de 1 a 2 comienza a descomponerse a medida que aumentan las dimensiones: las secuencias de Halton en particular desarrollan patrones de alias periódicos que distorsionan la cobertura. Comprender dónde se encuentran las fortalezas y debilidades de cada método es esencial para seleccionar el muestreador adecuado.

Muestreo de hipercubo latino: teoría e implementación

El muestreo de hipercubo latino fue introducido por McKay et al. (1989) como un compromiso entre la estratificación pura y la simplicidad del muestreo aleatorio.

Cómo funciona LHS

Para un problema de dimensión D con n muestras:

1. Divida cada dimensión marginal en n bins de intervalos iguales.
2. Para cada dimensión k, permute aleatoriamente el vector entero [1, 2, …, N] para determinar qué contenedor ocupa cada muestra.
3. Dibuje n muestras uniformes de cada dimensión, estratificadas dentro de sus contenedores asignados.
4. Permute de forma independiente los índices de muestra en todas las dimensiones, produciendo una matriz d × n.

La propiedad crucial es una muestra por bin a lo largo de cada dimensión. Esto garantiza que el histograma marginal para cada dimensión esté perfectamente estratificado.

Propiedades de convergencia

LHS tiene un factor de reducción de varianza en comparación con Monte Carlo simple. En condiciones de regularidad leve:

Var_LHS(ȳ) = Var_MC(ȳ) / N + O(N⁻²)

Esto significa que LHS converge aproximadamente a la misma tasa O(n⁻¹ᐟ²) en el error estándar pero con una constante significativamente más baja, aproximadamente un factor de n en la reducción de varianza.

Implementación de Python

import numpy as np
from scipy.stats import qmc

def lhs_sample(N, d, bounds=None):
    """
    Latin Hypercube Sampling using scipy.
    
    Parameters:
        N: number of samples
        d: number of dimensions
        bounds: list of (low, high) tuples for each dimension
    
    Returns:
        N x d array of LHS samples
    """
    if bounds is None:
        bounds = [(0, 1)] * d
    
    # Create the LHS engine
    engine = qmc.Lhs(d, scramble=False)
    
    # Generate samples in [0,1]
    samples = engine.random(N)
    
    # Apply bounds
    if bounds != [(0, 1)] * d:
        # Transform to specified bounds
        for i, (low, high) in enumerate(bounds):
            samples[:, i] = low + (high - low) * samples[:, i]
    
    return samples

# Example: 3D LHS with 1000 samples
np.random.seed(42)
d = 3
N = 1000
lhs_samples = lhs_sample(N, d)
print(f"LHS shape: {lhs_samples.shape}")
print(f"Marginal means (should be ~0.5): {lhs_samples.mean(axis=0)[:3]}")

La opción scramble=False anterior es intencional para la ilustración: en producción, debe usar scramble=True para romper la estructura de correlación artificial que introduce LHS determinista.

Cuasi-Monte Carlo: Secuencias Sobol y Scrambling

Las secuencias SOBOL fueron introducidas por I. M. Sobol en 1967 como la primera secuencia práctica de baja discrepancia. Usan un gráfico dirigido para generar dígitos en la base 2, produciendo una secuencia con discrepancia estelar provenientemente baja.

Por qué la codificación es esencial

Una secuencia SOBOL sin codificar tiene una estructura que repite cada 2^m muestras. Para integrandos periódicos o discontinuos, muy comunes en los flujos de trabajo de simulación, esto crea sesgo sistemático. La secuencia se alinea con las discontinuidades, lo que hace que el integrador se pierda características importantes.

La codificación aplica una permutación aleatoria a las representaciones binarias de los índices de secuencia, rompiendo esta periodicidad conservando la propiedad de baja discrepancia. El marco cuasi-monte carlo aleatorizado (RQMC) demuestra que las secuencias SOBOL codificadas, cuando se promedian sobre la codificación, producen estimadores imparciales con intervalos de confianza convergentes.

Recuento de muestras: potencias de dos

Las secuencias de Sobol se diseñaron para que exactamente 2^m puntos proporcionen una discrepancia óptima. Tomar n = 2^15 = 32.768 puntos produce una cobertura mucho mejor que tomar n = 30.000. Si necesita exactamente 30.000 muestras, usted:

1. Genera 32.768 puntos SOBOL y usa los primeros 30.000 (pierdes la optimización)
2. Use un sampler diferente para ese conteo específico
3. Utilice secuencias de Halton, que soportan n arbitrariamente

Scrambled Sobol en la práctica

import numpy as np
from scipy.stats import qmc

def scrambled_sobol(N, d, seed=None):
    """
    Generate scrambled Sobol QMC samples.
    
    Parameters:
        N: number of samples (ideally a power of 2)
        d: number of dimensions
        seed: random seed for scrambling
    
    Returns:
        N x d array of QMC samples
    """
    engine = qmc.QualifiedSobol(d, scramble=True, seed=seed)
    return engine.random(N)

# Example: 2^15 Sobol samples in 10 dimensions
np.random.seed(42)
d = 10
N = 2**15  # 32,768 - optimal for Sobol
sobol_samples = scrambled_sobol(N, d, seed=12345)
print(f"Scrambled Sobol shape: {sobol_samples.shape}")

Halton vs Sobol: Trade-offs prácticos

Las dos familias SUD más comunes son Halton y Sobol. Ambos son deterministas, ambos tienen baja discrepancia, pero se comportan de manera diferente en la práctica.

Secuencias de Halton

Las secuencias de Halton utilizan primos consecutivos como bases: 2, 3, 5, 7, 11, … esto asegura la independencia matemática entre las dimensiones. Las ventajas son:

• Recuentos de muestras arbitrarios — No hay necesidad de poderes de dos
• Fácil de generar — No se necesitan tablas polinómicas primitivas
• Rápido de calcular — Extracción de dígitos simple

Las desventajas están bien documentadas:

• Aliasing estructurado — Las proyecciones de mayor dimensión exhiben patrones periódicos
• Límites teóricos más pobres — La discrepancia crece más rápido que Sobol
• No hay teorías de codificación — La codificación es más difícil de definir para bases arbitrarias

Secuencias Sobol

Las secuencias sobol usan un polinomio primitivo fijo sobre GF(2), dando:

• Mejores límites teóricos: discrepancia más baja para la base-2
• Teoría enriquecida en la lucha — Owen Scrambling está bien definido y efectivo
• Optim en potencias de dos — garantías de optimización exactas

Las desventajas:

• Requiere tablas polinómicas primitivas — Debe calcularse previamente o almacenarse
• No flexible en el número de muestras — diseñado para n = 2 ^m
• Computacionalmente más pesado — Requiere búsquedas en tabla

el veredicto

Para la mayoría de las aplicaciones de simulación con dimensiones de 5 a 20, SOBOL con codificación supera a Halton. La razón es que la codificación rompe los patrones de alias que Halton desarrolla en dimensiones más altas. Cuando necesite recuentos de muestras arbitrarios y no pueda usar la codificación, Halton puede ser la mejor opción práctica.

Guía de comparación: MC vs LHS vs Sobol vs Halton

Guía de implementación de Python

Aquí hay un ejemplo de trabajo completo que compara los cuatro muestreadores en una función de prueba simple.

import numpy as np
from scipy.stats import qmc

def test_function(x):
    """
    Test integrand: sum of cosines along with a linear drift.
    True integral on [0,1]^d = sum of (sin(1)/1) for cos + 0.5*d for linear.
    """
    d = x.shape[-1]
    true_integral = d * np.sin(1) + 0.5 * d
    return np.cos(x) + x.mean(axis=-1, keepdims=True)

def estimate_integral(samples, func, true_value):
    """Estimate integral via QMC / MC averaging."""
    f_values = func(samples)
    estimate = f_values.mean(axis=-1).mean()
    std_error = f_values.std(ddof=1) / np.sqrt(len(f_values))
    return estimate, std_error

# Setup
d = 5  # dimension
np.random.seed(42)

# 1. Monte Carlo
N_mc = 2**15
mc_samples = np.random.uniform(0, 1, (N_mc, d))
mc_est, mc_se = estimate_integral(mc_samples, test_function, None)

# 2. Latin Hypercube (scrambled)
lhs_engine = qmc.Lhs(d, scramble=True)
lhs_samples = lhs_engine.random(N_mc)
lhs_est, lhs_se = estimate_integral(lhs_samples, test_function, None)

# 3. Scrambled Sobol QMC
sobol_engine = qmc.Sobol(d, scramble=True)
sobol_samples = sobol_engine.random(N_mc)
sobol_est, sobol_se = estimate_integral(sobol_samples, test_function, None)

# 4. Halton
halton_engine = qmc.Halton(d, scramble=False)
halton_samples = halton_engine.random(N_mc)
halton_est, halton_se = estimate_integral(halton_samples, test_function, None)

print(f"MC:      {mc_est:.6f} ± {mc_se:.6f}")
print(f"LHS:     {lhs_est:.6f} ± {lhs_se:.6f}")
print(f"Sobol:   {sobol_est:.6f} ± {sobol_se:.6f}")
print(f"Halton:  {halton_est:.6f} ± {halton_se:.6f}")

La producción muestra constantemente que los métodos QMC (SOBOL y LHS) producen estimaciones con una variación sustancialmente menor que la de Monte Carlo estándar, incluso con recuentos de muestras moderados.

Uso de límites para parámetros físicos

En los flujos de trabajo de simulación reales, sus parámetros tienen límites físicos. El módulo qmc maneja esto a través del parámetro expand:

from scipy.stats import qmc
import numpy as np

# Physical bounds for a 3D problem
bounds = [
    (300, 400),   # Temperature in Kelvin
    (1e-6, 1e-4), # Thermal conductivity, W/m·K
    (0.1, 0.9)    # Porosity, dimensionless
]

d = len(bounds)
engine = qmc.Sobol(d, scramble=True)

# Generate samples in [0,1], then expand to bounds
samples = engine.random(2**12)
samples = qmc.expand(samples, [b for b in bounds])

print(f"Sampled bounds check:")
for i, (low, high) in enumerate(bounds):
    print(f"  Dim {i}: min={samples[:, i].min():.6f}, max={samples[:, i].max():.6f}")

Trampas comunes y cómo evitarlas

Escolar 1: Eliminar los primeros puntos Sobol

Los primeros puntos SOBOL (índices de 0 a 10) tienen una cobertura deficiente. Si los elimina sin tener en cuenta el cambio de índice, crea un conjunto de muestras sesgado.

Fix: Siempre comience el muestreo de SOBOL en el índice 0 y deje que la secuencia completa se llene de forma natural. Si necesita una submuestra, use el medio o el final de la secuencia, no un truncamiento arbitrario.

Escollo 2: No-Potencias-de-2 Sobol

Tomar n = 50.000 puntos SOBOL desecha la propiedad de discrepancia óptima. La secuencia se diseñó para que n = 2^m proporcione una discrepancia baja exacta.

Fix: use n = 2^m y trunca (aceptable si solo necesita un subconjunto) o cambie a LHS o Halton para recuentos arbitrarios.

Escolar 3: Sobol sin codificar en producción

Las secuencias SOBOL sin codificar tienen una estructura periódica. Para las funciones periódicas (muy comunes en la propagación de la incertidumbre), esto crea un sesgo sistemático que no se promedia con más muestras.

FIX: Utilice siempre scramble=True. El algoritmo de codificación de Owen utilizado en scipy.stats.qmc produce estimaciones casi imparciales.

Escolar 4: Malinterpretando los intervalos de confianza de QMC

A diferencia de MC, QMC no produce intervalos de confianza estadísticamente válidos de forma predeterminada. Necesita RQMC (QMC aleatorizado) con múltiples codificaciones para construir intervalos válidos.

Fix: Ejecute K Scrambles de la misma secuencia QMC y use la varianza de muestra a través de scrambles para estimar los intervalos de confianza.

Escolar 5: Aplicar LHS a insumos altamente correlacionados

LHS asume la independencia de la dimensión. Si sus parámetros están correlacionados (por ejemplo, las propiedades del material derivadas de una distribución conjunta), LHS en los marginales ignora la estructura de correlación.

Fix: Use LHS basado en cópula o transforme en coordenadas independientes antes de aplicar LHS.

Lo que recomendamos

Seleccionar el método de muestreo correcto no es un problema de talla única. Aquí está nuestro marco de decisión práctica:

Cuándo utilizar el muestreo de hipercubo latino

• Necesita recuentos de muestras arbitrarios y desea beneficios de estratificación
• Tus dimensiones son ≤ 15 y no puedes garantizar los poderes de dos
• Necesitas una solución rápida y ligera de implementación — LHS es fácil de codificar desde cero
• Su función de simulación es moderadamente suave — LHS maneja las discontinuidades mejor que el SOBOL sin codificar

Cuándo usar SOBOL QMC revuelto

• Su integrando es suave (diferenciable, no hay discontinuidades fuertes)
• Puede usar n = 2^m Recuentos de muestras
• Necesitas la mejor tasa de convergencia teórica — O(n⁻¹) Beats O(n⁻¹ᐟ²)
• Está ejecutando campañas de UQ de producción donde el recuento de muestras está delimitado por el presupuesto de cálculo

Cuándo usar Halton

• Necesita recuentos de muestras arbitrarios y no puede usar LHS por alguna razón
• Sus dimensiones son ≤ 10 — Halton funciona bien en dimensiones bajas
• Estás prototipando y necesitas una secuencia rápida sin tablas

Cuándo usar Monte Carlo estándar

• Tu función tiene fuertes discontinuidades y no se puede codificar
• Sus dimensiones superan los 25 — Todos los métodos LDS se degradan y gana la simplicidad de MC
• Necesitas intervalos de confianza estadísticamente válidos sin tener que hacer una sobrecarga
• Su función es extremadamente irregular — La estructura SOBOL/LHS puede agregar sesgo

Una regla práctica práctica

Para la mayoría de los flujos de trabajo de simulación científica en 3 a 20 dimensiones con modelos avanzados suaves, Scrambled Sobol QMC con n = 2^12 a 2^15 muestras es el punto óptimo. Le brinda una reducción de varianza de 10 a 100 × sobre el Monte Carlo estándar en un recuento de muestras que se completa en horas en lugar de días en un solo nodo.

Resumen y próximos pasos

Los métodos de Monte Carlo para simulaciones científicas no se tratan solo de dibujar números aleatorios. La elección de la estrategia de muestreo (hipercubo latino, secuencias SOBOL o Halton) determina si su campaña de UQ es eficiente o derrochadora.

Las ideas clave son:

• Asuntos de discrepancia: las secuencias de baja discrepancia cubren el espacio de parámetros mucho más uniformemente que los puntos aleatorios
• La codificación es esencial para la producción SOBOL para evitar el sesgo sistemático
• Los poderes de dos son óptimos para SOBOL: diseñe su presupuesto de muestra en consecuencia
• Tasas de convergencia de QMC de O(n⁻¹) para funciones fluidas, lo convierten en la opción predeterminada para la mayoría de los flujos de trabajo de simulación

Nuestra Monte Carlo Methods for Scientific Simulations guía proporciona una introducción más amplia a los enfoques de Monte Carlo. Este artículo se enfoca en los métodos de baja discrepancia que le brindan la mayor eficiencia por muestra.

Para una implementación práctica, integre los patrones de código QMC anteriores en los scripts de campaña de simulación. Comience con n = 2^12, evalúe la convergencia monitoreando la varianza del estimador y amplíe hasta 2^15 si su presupuesto lo permite. Los resultados hablarán por sí mismos.

Guías relacionadas

• Cuantificación de incertidumbre y análisis de sensibilidad En Simulación Científica — Contexto UQ más amplio que incluye índices SOBOL y sensibilidad basada en la varianza
• Métodos Monte Carlo para simulaciones científicas: Guía de Python — Cobertura más amplia de Monte Carlo que incluye muestreo de importancia y métodos de adaptación
• verificación vs Validación en simulaciones científicas: una guía práctica — Cómo verificar que los resultados de su UQ sean correctos numéricamente
• Barridos de parámetros reproducibles: diseño de campañas de simulación para publicación — Documentación y mejores prácticas de reproducibilidad para campañas de muestreo

referencias

1. Kucherenko, M. et al. «Uso de las secuencias latinas de hipercubo y sobol para la cuantificación de la incertidumbre». arxiv:1505.02350, 2015. https://arxiv.org/abs/1505.02350
2. Chrisman, E. «Comparación del Día de Pi: Muestreo de Monte Carlo vs. Hipercubo Latino vs. Sobol de Muestreo». Analytica, 2022. https://analytica.com/blog/pi-day-comparison-monte-carlo-vs-latin-hypercube-vs-sobol-sampling/
3. escipia. «Referencia de muestreo cuasi-monte carlo». https://docs.scipy.org/doc/scipy/reference/stats.qmc.html
4. escipia. «Tutorial cuasi-Monte Carlo». https://scipy.github.io/devdocs/tutorial/stats/quasi_monte_carlo.html
5. Software QMC. «Documentación de QMCSoftware». https://qmcsoftware.github.io/qmcsoftware/
6. Diccionario Helmholtz-Uq. «Métodos cuasi-Monte Carlo». https://diccionario.helmholtz-uq.de/content/quasi_montecarlo_methods.html
7. ndcbe. «Cuaderno de cursos de cuantificación de incertidumbre basado en muestreo». https://ndcbe.github.io/cbe67701-uncertainty-quantification/07.01-sampling-based-uncertainty-quantification.html