Complementaria 3: Cadenas de Markov en Python II

Complementaria 3: Cadenas de Markov en Python II

#

Inventario de farmacia con cadena de Markov

Esta complementaria tiene como objetivo resolver el ejercicio de CaliPharm, que modela el sistema de inventario de una farmacia.

El modelo se basa en una cadena de Markov de tiempo discreto donde:

  • La demanda diaria sigue una distribución Poisson.

  • La reposición del medicamento depende del inventario final.

  • Hay colaboración externa si el inventario se agota por completo.

Parte A. Modelar el sistema de inventario de CaliPharm como una cadena de Markov

Variable de estado:

\( X_n \): Número de unidades del medicamento disponibles al inicio del día n.

Espacio de estados:

\( S_X = \{0, 1, 2, 3, 4, 5\} \):

Otras definiciones:

Adicionalmente, se definen las siguientes variables aleatorias:

  • La variable aleatoria \(D\) que representa la demanda diaria del medicamento, la cual se distribuye según una distribución de Poisson con parámetro \(\lambda = 3\):

\( D \sim \text{Poisson}(\lambda = 3) \)

\( \mathbb{P}(D = k) = \frac{3^k e^{-3}}{k!}, \quad k = 0, 1, 2, \ldots \)

\( \mathbb{P}(D \geq k) = \sum_{r = k}^{\infty} \frac{3^r e^{-3}}{r!} \)

  • La variable aleatoria \(Y\) que representa la cantidad de unidades que podrían llegar al inventario desde hospitales aliados cuando la farmacia entra en cuarentena por haber finalizado el día con inventario cero.

\( Y \sim \text{Binomial}(n = 5, p = 0.07) \) \( \mathbb{P}(Y = k) = \binom{5}{k}(0.07)^k(0.93)^{5 - k}, \quad k = 0, 1, 2, 3, 4, 5 \)

La política de transición es:

\( \mathbb{P}_{i \to j} = \left\{ \begin{array}{ll} P(D = i - j), & i \geq 2 , j \neq 0 \\ P(D \geq i), & i \geq 2 , j = 0 \\ P(D = 5 - j), & 1 \leq i \lt 2 , j \neq 0 \\ P(D \geq 5), & 1 \leq i \lt 2 , j = 0 \\ P(Y = j), & i = 0 , j \in \{0,1,\dots,5\} \\ 0, & \text{d.l.c.} \end{array} \right. \)

A continuación, se implementa el modelo en Python.

import numpy as np
from scipy.stats import poisson, binom

# Parámetros
demanda_lambda = 3
hospitales_n = 5
prob_hospitales = 0.07
inventario_max = 5

# Espacio de estados
estados = [i for i in range(0, inventario_max+1)]

# Matriz de transición
P = np.zeros((len(estados), len(estados)))

# Llenar la matriz de transición
for i in estados:
    for j in estados:
        if i>=2 and j!=0:
            P[i,j] = poisson.pmf(i - j, demanda_lambda)
        elif i>=2 and j==0:
            P[i,j] = poisson.sf(i-1, demanda_lambda)
        elif 1<=i and i<2 and j!=0:
            P[i,j] = poisson.pmf(inventario_max - j, demanda_lambda)
        elif 1<=i and i<2 and j==0:
            P[i,j] = poisson.sf(inventario_max-1, demanda_lambda)
        elif i==0:
            P[i,j] = binom.pmf(j, hospitales_n, prob_hospitales)
        

print("Matriz de transición P:")
print(P.round(2))

Matriz de transición P:
[[0.7  0.26 0.04 0.   0.   0.  ]
 [0.18 0.17 0.22 0.22 0.15 0.05]
 [0.8  0.15 0.05 0.   0.   0.  ]
 [0.58 0.22 0.15 0.05 0.   0.  ]
 [0.35 0.22 0.22 0.15 0.05 0.  ]
 [0.18 0.17 0.22 0.22 0.15 0.05]]

Para crear la cadena de markov utilizaremos la instancia dtmc la libreria jmarkov, ya que estamos creando una cadena de tiempo discreto.

from jmarkov.dtmc import dtmc

# Crear objeto de cadena de Markov
mc = dtmc(P)

Parte B. Responder preguntas de interés

Literal a. ¿Cuál es la probabilidad de que el inventario sea igual a 5 dentro de 2 días, si hoy hay 1 unidad?

Sea \(X_0 = 1\), queremos calcular:

\( P(X_2 = 5 \mid X_0 = 1) \)

Esta probabilidad corresponde a una transición de dos pasos, por lo tanto: \( P(X_2 = 5 \mid X_0 = 1) = [P^2]_{1,5} \)

# Probabilidad de estar en estado 5 en dos días si hoy estamos en estado 1
alpha = np.zeros(len(estados))
alpha[1] = 1  # Estado inicial = 1

dist_paso2 = mc.transient_probabilities(n=2, alpha=alpha)
print(f"La probabilidad de que el inventario sea igual a 5 en 2 días, dado que hoy hay 1 unidad, es: {dist_paso2[5]:.4f}")

La probabilidad de que el inventario sea igual a 5 en 2 días, dado que hoy hay 1 unidad, es: 0.0108
Literal b. ¿Qué pasaría si el umbral de reposición cambiara de 2 a 3? ¿Se reduciría el riesgo de quiebre?

Para responder esta pregunta podemos ajustar la política de transición para que la reposición se active cuando el inventario sea menor a 3 (en vez de menor a 2). Esto implica modificar las condiciones de la matriz de transición.

Quedaría como: \( \mathbb{P}_{i \to j} = \left\{ \begin{array}{ll} P(D = i - j), & i \geq 3 , j \neq 0 \\ P(D \geq i), & i \geq 3 , j = 0 \\ P(D = 5 - j), & 1 \leq i \lt 3 , j \neq 0 \\ P(D \geq 5), & 1 \leq i \lt 3 , j = 0 \\ P(Y = j), & i = 0 , j \in \{0,1,\dots,5\} \\ 0, & \text{d.l.c.} \end{array} \right. \)

Luego de modificar la política, podemos comparar la probabilidad de que el sistema llegue al estado 0 (sin inventario) en el largo plazo.

Para esto, utilizaremos el vector \(\pi\), que representa la distribución estacionaria de la cadena de Markov, es decir, la proporción de tiempo (o frecuencia a largo plazo) que el sistema pasa en cada uno de los estados cuando se deja evolucionar por un número suficientemente grande de pasos.

En este contexto, la probabilidad de que el sistema esté en quiebre (inventario = 0) en el largo plazo corresponde al valor \(\pi[0]\).

Por tanto, comparar \(\pi[0]\) bajo distintas políticas de reposición permite analizar cuál configuración reduce más el riesgo de quiebre permanente del inventario. Una política con menor \(\pi[0]\) indica un sistema más confiable frente a agotamientos.

# Matriz de transición ajustando el umbral de reposición de <2 a <3
P_alt = np.zeros((len(estados), len(estados)))

for i in estados:
    for j in estados:
        if i>=3 and j!=0:
            P_alt[i,j] = poisson.pmf(i - j, demanda_lambda)
        elif i>=3 and j==0:
            P_alt[i,j] = poisson.sf(i-1, demanda_lambda)
        elif 1<=i and i<3 and j!=0:
            P_alt[i,j] = poisson.pmf(inventario_max - j, demanda_lambda)
        elif 1<=i and i<3 and j==0:
            P_alt[i,j] = poisson.sf(inventario_max-1, demanda_lambda)
        elif i==0:
            P_alt[i,j] = binom.pmf(j, hospitales_n, prob_hospitales)


# Comparar probabilidad de inventario agotado a largo plazo

mc_old = dtmc(P)        # Política con umbral <2
mc_new = dtmc(P_alt)    # Política con umbral <3

pi_old = mc_old.steady_state()
pi_new = mc_new.steady_state()

print(f"π[0] con umbral <2 = {pi_old[0]:.4f}")
print(f"π[0] con umbral <3 = {pi_new[0]:.4f}")

π[0] con umbral <2 = 0.5648
π[0] con umbral <3 = 0.4759
Literal c. ¿Cuánto impacta la probabilidad de colaboración p = 0.07 en la confiabilidad del sistema? ¿Cómo cambiarían las probabilidades si esta subiera a 0.2?

Para este análisis se construye una nueva matriz de transición aumentando la probabilidad de préstamo de los hospitales. Luego se compara la probabilidad de que haya inventario disponible al día siguiente de una cuarentena (es decir, cuando hoy hay 0 unidades).

# Aumentar la probabilidad de colaboración a 0.2
nueva_prob_hospitales = 0.2

# Recalcular fila i=0 de la matriz
P_colab = P.copy()
for j in estados:
    P_colab[0,j] = binom.pmf(j, hospitales_n, nueva_prob_hospitales)

# Crear cadena con nueva p de colaboración
mc_alt = dtmc(P_colab)

# Vector inicial desde estado 0
alpha_cuarentena = np.zeros(len(estados))
alpha_cuarentena[0] = 1

# Distribución después de 1 paso (reapertura)
dist_reapertura = mc_alt.transient_probabilities(n=1, alpha=alpha_cuarentena)

# Probabilidad de tener inventario > 0
prob_recuperacion = 1 - dist_reapertura[0]

# Comparación con la política original
dist_reapertura_orig = mc.transient_probabilities(n=1, alpha=alpha_cuarentena)
prob_recuperacion_orig = 1 - dist_reapertura_orig[0]

print(f"Comparación:")
print(f"Recuperación con p = 0.07: {prob_recuperacion_orig:.4f}")
print(f"Recuperación con p = 0.20: {prob_recuperacion:.4f}")

Comparación:
Recuperación con p = 0.07: 0.3043
Recuperación con p = 0.20: 0.6723

Universidad de los Andes | Vigilada Mineducación. Reconocimiento como Universidad: Decreto 1297 del 30 de mayo de 1964. Reconocimiento personería jurídica: Resolución 28 del 23 de febrero de 1949 Minjusticia. Departamento de Ingeniería Industrial Carrera 1 Este No. 19 A 40 Bogotá, Colombia Tel. (57.1) 3324320 | (57.1) 3394949 Ext. 2880 /2881 http://industrial.uniandes.edu.co