← Volver al blog

Monte Carlo con cualquier distribución: por qué funciona y cómo se hace

25 de abril de 2026 10 min de lectura monte carlo simulation statistics distributions python weibull mixture bootstrap

Cuando se compara una simulación Monte Carlo con Wi Normal vs Lognormal con los parámetros típicos de un molino de bolas (media 14 kWh/t, σ 2 kWh/t), las curvas quedan casi superpuestas. La pregunta obvia es: "si los resultados son tan parecidos, ¿dónde está el problema de asumir normalidad?"

La respuesta tiene dos partes y este artículo aborda ambas:

  1. Por qué dos distribuciones con la misma media y σ pueden verse casi idénticas (spoiler: coeficiente de variación demasiado bajo).
  2. Por qué Monte Carlo funciona con cualquier distribución — no solo normal — y cómo se implementa esa flexibilidad en Python.

Vamos a ver tres casos donde la diferencia entre asumir normal y usar la distribución correcta no es marginal, es decisiva.

1. ¿Por qué dos distribuciones pueden verse tan similares?

Cuando una variable tiene media 14 kWh/t y desviación 2 kWh/t, su coeficiente de variación (CV = σ/μ) es solo 14 %. Para valores de CV bajos, la lognormal es prácticamente indistinguible de la normal.

El siguiente código lo muestra de forma directa:

import numpy as np
import matplotlib.pyplot as plt

rng = np.random.default_rng(42)
N = 50_000
mu = 100  # mantener la media fija

fig, axes = plt.subplots(1, 3, figsize=(15, 4))
for ax, cv in zip(axes, [0.15, 0.50, 1.00]):
    sigma = cv * mu
    # Parámetros equivalentes de la lognormal
    mu_log = np.log(mu**2 / np.sqrt(mu**2 + sigma**2))
    sigma_log = np.sqrt(np.log(1 + (sigma/mu)**2))

    norm = rng.normal(mu, sigma, N)
    logn = rng.lognormal(mu_log, sigma_log, N)

    ax.hist(norm, bins=80, alpha=0.5, density=True, label="Normal")
    ax.hist(logn, bins=80, alpha=0.5, density=True, label="Lognormal")
    ax.set_title(f"CV = {cv:.2f}")
    ax.legend()
plt.tight_layout()
plt.show()

Normal vs Lognormal a tres niveles de CV

  • CV = 0.15 → casi idénticas (el régimen típico de variables controladas de planta).
  • CV = 0.50 → se separan; la lognormal muestra cola derecha.
  • CV = 1.00 → distribuciones radicalmente distintas; la normal incluso genera valores negativos.

Conclusión operativa: la elección de distribución importa cuando la variabilidad es alta relativa a la media. Variables de proceso con CV > 30 % (consumos de reactivos con upsets, tiempos a falla, leyes de mineral) están en ese régimen y asumir normalidad produce errores significativos en los percentiles.

2. El principio: Monte Carlo no exige ninguna distribución

El método Monte Carlo solo requiere dos cosas:

  1. Poder muestrear aleatoriamente de la distribución de cada variable de entrada.
  2. Poder evaluar el modelo determinístico en cada muestra.

Eso es todo. No hay ningún teorema que obligue a usar normal. La teoría no distingue entre normal, Weibull, beta, bimodal, empírica o una mezcla personalizada. Mientras tengas un generador de muestras, Monte Carlo funciona.

La confusión viene de tres lugares:

  • La normal es la distribución más fácil de explicar en cursos introductorios.
  • El Teorema del Límite Central aparece mucho en estadística, y se generaliza en mala forma.
  • Excel históricamente ofrecía solo NORM.INV de forma simple, lo que propagó el hábito.

Python rompe esa limitación: puedes muestrear de prácticamente cualquier distribución con una línea de código.

3. Cuatro formas de muestrear de distribuciones no normales

3.1 Generadores directos de NumPy

numpy.random.Generator incluye ya una decena de distribuciones listas:

rng = np.random.default_rng(42)
N = 10_000

muestras_normal       = rng.normal(loc=100, scale=15, size=N)
muestras_lognormal    = rng.lognormal(mean=4.5, sigma=0.5, size=N)
muestras_exponencial  = rng.exponential(scale=2000, size=N)
muestras_weibull      = rng.weibull(a=1.8, size=N) * 3000
muestras_gamma        = rng.gamma(shape=2, scale=1.5, size=N)
muestras_beta         = rng.beta(a=18, b=2, size=N)
muestras_triangular   = rng.triangular(left=80, mode=100, right=150, size=N)
muestras_uniform      = rng.uniform(low=0.1, high=0.9, size=N)
muestras_poisson      = rng.poisson(lam=5, size=N)

Cada llamada devuelve un vector de N muestras. Eso es todo lo que Monte Carlo necesita.

3.2 Catálogo universal de scipy.stats

SciPy expone más de 80 distribuciones con una API uniforme: .rvs(..., size=N) muestrea:

from scipy import stats

# Weibull de dos parámetros (forma k, escala λ)
muestras_w = stats.weibull_min.rvs(c=1.5, scale=2000, size=N, random_state=rng)

# Beta para variables acotadas en [0, 1]
muestras_b = stats.beta.rvs(a=5, b=2, size=N, random_state=rng)

# Distribución Pareto (colas pesadas)
muestras_p = stats.pareto.rvs(b=1.16, scale=50, size=N, random_state=rng)

# Distribución t de Student (colas más pesadas que normal)
muestras_t = stats.t.rvs(df=5, loc=100, scale=15, size=N, random_state=rng)

# Normal truncada (normal restringida a un rango físico)
muestras_nt = stats.truncnorm.rvs(a=-2, b=2, loc=100, scale=10, size=N, random_state=rng)

Misma lógica, distinta distribución. Monte Carlo es agnóstico.

3.3 Transformada inversa: cualquier distribución personalizada

Si la distribución que necesitas no está en NumPy ni SciPy, puedes construirla usando el método de la transformada inversa:

Si U ~ Uniforme(0,1) y F es la función de distribución acumulada (CDF), entonces X = F⁻¹(U) tiene la distribución deseada.

Ejemplo — distribución de Pareto escrita "a mano":

# CDF: F(x) = 1 - (x_min / x)^alpha  para x >= x_min
# Inversa: x = x_min / (1 - u)^(1/alpha)

def sample_pareto(x_min, alpha, size, rng):
    u = rng.uniform(0, 1, size)
    return x_min / (1 - u)**(1/alpha)

muestras = sample_pareto(x_min=10, alpha=1.5, size=N, rng=rng)

Esto funciona para cualquier distribución de la que sepas escribir la CDF inversa. Muy útil para distribuciones empíricas de ingeniería que no tienen nombre formal.

3.4 Bootstrap: muestrear directamente de datos históricos

Si no sabes qué distribución tienes pero dispones de datos reales, no necesitas ajustar nada — muestrea directamente de los datos:

# Datos históricos de consumo de colector (ejemplo)
datos_historicos = np.array([142, 138, 155, 161, 150, 147, 310, 145, 152, 148,
                              139, 159, 143, 280, 151, 146, 154, 260, 149, 158])

muestras = rng.choice(datos_historicos, size=N, replace=True)

Esto captura automáticamente cualquier asimetría, bimodalidad o cola pesada que exista en los datos, sin asumir ninguna forma. El único requisito es tener suficientes datos históricos representativos (típicamente n ≥ 50).

4. Caso 1: donde la diferencia SÍ es dramática

Dimensionemos un tanque de almacenamiento para un reactivo cuyo consumo histórico tiene:

  • Media = 100 g/t
  • Desviación estándar = 100 g/t (CV = 1.0)

Queremos el P99 del consumo para garantizar 99 % de cobertura de stock.

N = 100_000
mu, sigma = 100, 100

# Caso A: asumir Normal
x_norm = rng.normal(mu, sigma, N)
x_norm = np.clip(x_norm, 0, None)  # evitar negativos

# Caso B: usar Lognormal (misma media, misma desviación)
mu_log = np.log(mu**2 / np.sqrt(mu**2 + sigma**2))
sigma_log = np.sqrt(np.log(1 + (sigma/mu)**2))
x_logn = rng.lognormal(mu_log, sigma_log, N)

print(f"Normal    → media={x_norm.mean():.1f}, P95={np.percentile(x_norm,95):.1f}, "
      f"P99={np.percentile(x_norm,99):.1f}, P99.9={np.percentile(x_norm,99.9):.1f}")
print(f"Lognormal → media={x_logn.mean():.1f}, P95={np.percentile(x_logn,95):.1f}, "
      f"P99={np.percentile(x_logn,99):.1f}, P99.9={np.percentile(x_logn,99.9):.1f}")

Resultado típico:

Métrica Normal Lognormal Δ (%)
Media 105.2 100.1 —
P95 266 284 +7 %
P99 333 485 +46 %
P99.9 410 810 +98 %

La media y la desviación son iguales por construcción, pero la lognormal predice eventos extremos que son el doble en el P99.9. Si dimensionas el tanque con el P99 asumiendo normal, el tanque se queda corto el 1 % del tiempo por un margen considerable — y en una operación que trabaja 8760 horas al año, eso son 87 horas de paro al año por falta de reactivo.

Con CV = 1.0 la historia cambia radicalmente: la diferencia del P99 entre normal y lognormal es casi del 50 %, no marginal.

5. Caso 2: donde la normal NO TIENE SENTIDO desde el inicio

Predigamos el tiempo a falla de los revestimientos de un molino. La distribución estándar de confiabilidad es Weibull, no normal. La normal ni siquiera es candidato válido porque:

  • Permite tiempos negativos (imposible).
  • Tiene colas simétricas (el envejecimiento de componentes es fundamentalmente asimétrico).
# Parámetros típicos de liners: Weibull forma k=2.5, escala λ=3500 h
from scipy import stats

N = 50_000
tiempo_falla = stats.weibull_min.rvs(c=2.5, scale=3500, size=N, random_state=rng)

print(f"MTTF (media): {tiempo_falla.mean():.0f} h")
print(f"Mediana:      {np.median(tiempo_falla):.0f} h")
print(f"P10:          {np.percentile(tiempo_falla, 10):.0f} h")
print(f"P90:          {np.percentile(tiempo_falla, 90):.0f} h")

# Probabilidad de sobrevivir a las 2000 horas
prob_sup_2000 = np.mean(tiempo_falla > 2000)
print(f"P(falla > 2000 h) = {prob_sup_2000:.3f}")

Esto es Monte Carlo puro, sin una normal a la vista. La lógica es idéntica a cualquier otra simulación: muestrear → propagar → analizar. Solo cambió el generador de muestras.

Podemos extender a una simulación de disponibilidad anual del molino:

def simular_disponibilidad(n_anos=10_000, horas_ano=8760, rng=None):
    rng = rng or np.random.default_rng()
    horas_disponibles = np.zeros(n_anos)

    for i in range(n_anos):
        t = 0
        horas_op = 0
        while t < horas_ano:
            tiempo_hasta_falla = stats.weibull_min.rvs(
                c=2.5, scale=3500, random_state=rng)
            tiempo_reparacion = stats.lognorm.rvs(
                s=0.4, scale=24, random_state=rng)  # ~24 h ± asimétrica

            if t + tiempo_hasta_falla >= horas_ano:
                horas_op += horas_ano - t
                break
            horas_op += tiempo_hasta_falla
            t += tiempo_hasta_falla + tiempo_reparacion

        horas_disponibles[i] = horas_op

    return horas_disponibles / horas_ano  # disponibilidad fraccionaria

disp = simular_disponibilidad(n_anos=5000, rng=rng)
print(f"Disponibilidad media: {disp.mean():.3f}")
print(f"P05 de disponibilidad: {np.percentile(disp, 5):.3f}")

Aquí conviven Weibull (tiempo a falla) y Lognormal (tiempo de reparación) en la misma simulación. Ninguna es normal. Monte Carlo no se inmuta.

6. Caso 3: la distribución no es paramétrica — es una mezcla

A veces el fenómeno real no tiene nombre estadístico conocido. Por ejemplo, un reactivo se consume:

  • 85 % del tiempo en régimen normal (modo operación estándar)
  • 15 % del tiempo en régimen upset (problemas de mineral, flotabilidad baja, compensación con más reactivo)

No hay una distribución paramétrica única que describa esto. Es una mezcla (mixture model). Monte Carlo la maneja directamente:

def muestrear_consumo_mezcla(N, rng):
    # ¿Este punto está en régimen normal o en upset?
    estado = rng.uniform(0, 1, N) < 0.85   # 85% normal

    # Régimen normal: alrededor de 150 g/t, baja dispersión
    consumo_normal = rng.normal(150, 15, N)

    # Régimen upset: alrededor de 280 g/t, más dispersión
    consumo_upset = rng.lognormal(np.log(280), 0.25, N)

    # Combinar según el estado
    return np.where(estado, consumo_normal, consumo_upset)

consumos = muestrear_consumo_mezcla(100_000, rng)

print(f"Media:  {consumos.mean():.1f} g/t")
print(f"P50:    {np.percentile(consumos, 50):.1f} g/t")
print(f"P95:    {np.percentile(consumos, 95):.1f} g/t")
print(f"P99:    {np.percentile(consumos, 99):.1f} g/t")

Distribución bimodal: normal + upsets

El histograma muestra dos modos — uno alrededor de 150 y otro alrededor de 280. Ninguna distribución paramétrica estándar captura esto; una normal forzada aquí daría resultados absurdos. Pero la simulación Monte Carlo funciona sin problema, porque el único requisito es poder generar muestras.

Variantes más sofisticadas de este mismo patrón:

  • Tres o más regímenes (normal / degradado / crítico).
  • Probabilidades de régimen que dependen de otras variables (ley de cabeza, hora del día, estación del año).
  • Bootstrapping directo desde datos históricos si no quieres asumir ninguna forma.

7. El patrón general de Monte Carlo con cualquier distribución

La receta que engloba todos los casos anteriores:

def monte_carlo_generico(N, modelo_determinista, generadores_inputs, rng):
    """
    N                   : número de simulaciones
    modelo_determinista : función(dict_inputs) → output
    generadores_inputs  : dict {nombre_variable: función(N, rng) → muestras}
    """
    inputs = {nombre: gen(N, rng) for nombre, gen in generadores_inputs.items()}
    outputs = np.zeros(N)
    for i in range(N):
        muestra_i = {k: v[i] for k, v in inputs.items()}
        outputs[i] = modelo_determinista(muestra_i)
    return outputs

Un ejemplo mezclando las cuatro formas de muestreo — normal, lognormal, Weibull y empírico:

datos_ley_historica = np.array([...])  # tu serie histórica de leyes

generadores = {
    "potencia":  lambda N, rng: rng.normal(500, 15, N),
    "wi":        lambda N, rng: rng.lognormal(np.log(14), 0.25, N),
    "t_falla":   lambda N, rng: stats.weibull_min.rvs(c=2.5, scale=3500,
                                                      size=N, random_state=rng),
    "ley":       lambda N, rng: rng.choice(datos_ley_historica, size=N, replace=True),
}

def mi_modelo(x):
    # tu lógica de negocio aquí — recuperación, tonelaje, VAN, etc.
    return x["potencia"] * x["ley"] / x["wi"]

resultados = monte_carlo_generico(50_000, mi_modelo, generadores, rng)

Ninguna variable necesita ser normal. Ese es el punto de todo el artículo.

8. Cierre

Monte Carlo no está atado a la distribución normal. Trabaja con lo que le entregues: normal, lognormal, Weibull, beta, gamma, triangular, Pareto, t de Student, mezclas, bootstrap de datos históricos, o distribuciones personalizadas construidas por transformada inversa.

La normalidad es una elección por defecto cuando no se tiene mejor información — pero en ingeniería de procesos casi siempre se tiene mejor información, ya sea por la física del fenómeno (confiabilidad → Weibull) o por los datos históricos disponibles (ajuste MLE o bootstrap directo).

El flujo robusto de Monte Carlo es:

  1. Identificar la distribución de cada variable (ajuste + validación — cubierto en el artículo anterior sobre identificación de distribuciones).
  2. Elegir el generador correcto según el caso: NumPy, SciPy, transformada inversa, o bootstrap.
  3. Propagar por el modelo determinista.
  4. Reportar percentiles, no solo la media.

Si tus datos tienen CV alta, son asimétricos, tienen outliers frecuentes, o vienen de mezclas de regímenes — asumir normalidad introduce errores que pueden ser de dos dígitos en los percentiles extremos. En diseño de planta, inventarios de reactivos o disponibilidad mecánica, eso se traduce directamente en sobredimensionamiento innecesario o subdimensionamiento peligroso.

Python ofrece todo el toolbox para hacerlo bien. No hay excusa técnica para seguir asumiendo normal por defecto.

← Volver al blog