Simulación de cadenas de Markov

01 · Punto de partida

De la matriz a las historias

Con P, π(0) y N se generan trayectorias: sortear X₀, luego fila por fila hasta N. Con R réplicas se estima lo que πP calcula: fracciones, tiempos y repartos.

La simulación no reemplaza la teoría: la verifica y la lleva a indicadores del negocio.

• ¿Cómo se sortea X₀ y cada paso con su fila?
• ¿Cuántas réplicas R y qué horizonte N?
• ¿Qué burn-in se descarta para estimar π?
• ¿Coincide lo simulado con π = πP?

02 · Definición

La receta del simulador

Cuatro pasos, siempre igual:

Servidor 3 modos (P del Tema 24) con π ≈ [0,63 – 0,29 – 0,08]: en 5000 pasos simulados las fracciones caen ahí.

03 · Anatomía

Semilla, vectorización y chequeos

Cuatro decisiones de código:

Semilla

Una por réplica (maestra + índice). Repetible y distinta a la vez.

Bucle vs vector

Bucle claro hasta miles de pasos; vectorizado (numpy) a millones.

Burn-in

Descartar el arranque antes de promediar. Sin él, π(0) sesga.

Chequeos

P válida, π(0) válida, fracciones ≈ π. Fallar uno invalida el experimento.

04 · Ejemplos

Simuladores que entregan números

La misma receta, distintos P:

1. 1
   Servidor 3 modos.

   Fracción caída ≈ 8 % para SLA.

2. 2
   Clima 30 días.

   Días de lluvia esperados = 30·π₁.

3. 3
   Ruina con freno.

   Réplicas hasta absorber: reparto y tiempos.

4. 4
   Cola con tope.

   Pérdidas por tope: fracción en K.

5. 5
   Cliente mensual.

   Retención a 12 meses: entrada de P¹².

05 · Representación en Python

El simulador genérico K×K

En Python, un simulador genérico recibe P, π(0), N y semilla y devuelve la trayectoria. Con R llamadas salen las réplicas.

import random

P = [[0.6, 0.35, 0.05],
     [0.4, 0.5, 0.1],
     [0.9, 0.0, 0.1]]

def sortear(dists, u):
    ac = 0.0
    for j, p in enumerate(dists):
        ac += p
        if u < ac:
            return j
    return len(dists) - 1

def trayectoria(pi0, N, semilla=7):
    random.seed(semilla)
    x = sortear(pi0, random.random())
    hist = [x]
    for _ in range(N):
        x = sortear(P[x], random.random())
        hist.append(x)
    return hist

print(trayectoria([1.0, 0.0, 0.0], 10))

Una historia de 11 valores desde libre: salta entre 0 y 1 con visitas raras al 2.

Estimar π con burn-in

R réplicas, descarte y fracciones:

import random

P = [[0.6, 0.35, 0.05],
     [0.4, 0.5, 0.1],
     [0.9, 0.0, 0.1]]

def trayectoria(pi0, N, semilla):
    random.seed(semilla)
    u = random.random()
    ac, x = 0.0, 2
    for j, p in enumerate(pi0):
        ac += p
        if u < ac:
            x = j
            break
    hist = [x]
    for _ in range(N):
        u, ac = random.random(), 0.0
        for j, p in enumerate(P[x]):
            ac += p
            if u < ac:
                x = j
                break
        hist.append(x)
    return hist

R, N, burn = 200, 200, 40
tot = [0, 0, 0]
for r in range(R):
    h = trayectoria([1.0, 0.0, 0.0], N, semilla=1000 + r)
    for x in h[burn:]:
        tot[x] += 1
s = sum(tot)
print([round(v / s, 3) for v in tot])

06 · Exploración

Laboratorio: estimar π simulando

Servidor 3 modos con π ≈ [0,63 – 0,29 – 0,08]. Movete R (réplicas) y N (horizonte, burn-in 20 %): las barras son la fracción simulada y los puntos la π teórica. Con poco R·N el error es visible; con mucho se pegan.

Cada barra lleva su punto teórico: si la barra lo tapa, el simulador está sano. El caído (8 %) necesita más muestra que el libre.

Preguntas para explorar

1. Con R = 20 y N = 50, ¿el error es chico? ¿Qué estado falla más (el raro)?
2. Duplicá R y N. ¿El error cae a la mitad o a ~70 %? (Pista: √2.)
3. ¿Qué pasa si el burn-in fuera 0 con inicio seguro en 0?

07 · Comprensión

Confusiones frecuentes

08 · Práctica guiada

Ejercicios con Python

Ejercicio 1: una trayectoria a mano

Generá 15 pasos 2×2 desde 0 con semilla fija e imprimí la historia.

import random

P = [[0.7, 0.3], [0.5, 0.5]]
random.seed(21)
x, hist = 0, [0]
for _ in range(15):
    u, ac = random.random(), 0.0
    for j, p in enumerate(P[x]):
        ac += p
        if u < ac:
            x = j
            break
    hist.append(x)
print(hist)

Ejercicio 2: fracción sin burn-in vs con

Sin ejecutar primero: desde [1, 0] con N = 10, ¿la fracción en 0 sobrestima o subestima π₀? Comprobalo.

import random

P = [[0.7, 0.3], [0.5, 0.5]]
random.seed(3)
x, h = 0, [0]
for _ in range(10):
    x = 1 if random.random() < (0.3 if x == 0 else 0.5) else 0
    h.append(x)
print("todo:", round(sum(1 for v in h if v == 0) / len(h), 3))
print("sin 3:", round(sum(1 for v in h[3:] if v == 0) / len(h[3:]), 3))

Ejercicio 3: muestra útil

Con R = 100, N = 200 y burn 20 %, ¿cuántos estados útiles hay? ¿Y el error típico en el caído (π = 0,08)?

import math
M = 100 * 160
print("muestra:", M)
print("se:", round(math.sqrt(0.08 * 0.92 / M), 4))

09 · Síntesis

Ideas para recordar

En el próximo tema modelaremos sistemas reales: del problema a P.