Skip to content

Commit 3a7f4f5

Browse files
mgabyurangomgabyurango
committed
Update documentation
1 parent a711ae0 commit 3a7f4f5

File tree

3 files changed

+572
-14
lines changed

3 files changed

+572
-14
lines changed

‎_sources/content/complementarias/Complementaria3.ipynb‎

Lines changed: 323 additions & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -7,27 +7,339 @@
77
"### <h3 style=\"color: #ADD8E6;\">Complementaria 3: Cadenas de Markov en Python II</h3>"
88
]
99
},
10+
{
11+
"cell_type": "markdown",
12+
"metadata": {},
13+
"source": [
14+
"<h3 style=\"color: #ADD8E6;\">Inventario de farmacia con cadena de Markov</h3>\n",
15+
"\n",
16+
"Esta complementaria tiene como objetivo resolver el ejercicio de CaliPharm, que modela el sistema de inventario de una farmacia.\n",
17+
"\n",
18+
"El modelo se basa en una cadena de Markov de tiempo discreto donde:\n",
19+
"- La demanda diaria sigue una distribución Poisson.\n",
20+
"- La reposición del medicamento depende del inventario final.\n",
21+
"- Hay colaboración externa si el inventario se agota por completo.\n"
22+
]
23+
},
24+
{
25+
"cell_type": "markdown",
26+
"metadata": {},
27+
"source": [
28+
"<h4 style=\"color: #ADD8E6;\">Parte A. Modelar el sistema de inventario de CaliPharm como una cadena de Markov</h4>\n",
29+
"\n",
30+
"**Variable de estado:**\n",
31+
"\n",
32+
"$ X_n $: Número de unidades del medicamento disponibles al inicio del día n.\n",
33+
"\n",
34+
"**Espacio de estados:** \n",
35+
"\n",
36+
"$ S_X = \\{0, 1, 2, 3, 4, 5\\} $:\n",
37+
"\n",
38+
"**Otras definiciones:** \n",
39+
"\n",
40+
"Adicionalmente, se definen las siguientes variables aleatorias:\n",
41+
"\n",
42+
"- 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$:\n",
43+
"$$\n",
44+
"D \\sim \\text{Poisson}(\\lambda = 3)\n",
45+
"$$\n",
46+
"\n",
47+
"$$\n",
48+
"\\mathbb{P}(D = k) = \\frac{3^k e^{-3}}{k!}, \\quad k = 0, 1, 2, \\ldots\n",
49+
"$$\n",
50+
"\n",
51+
"$$\n",
52+
"\\mathbb{P}(D \\geq k) = \\sum_{r = k}^{\\infty} \\frac{3^r e^{-3}}{r!}\n",
53+
"$$\n",
54+
"\n",
55+
"\n",
56+
"\n",
57+
"- 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.\n",
58+
"\n",
59+
"$$ Y \\sim \\text{Binomial}(n = 5, p = 0.07) $$\n",
60+
"$$\n",
61+
"\\mathbb{P}(Y = k) = \\binom{5}{k}(0.07)^k(0.93)^{5 - k}, \\quad k = 0, 1, 2, 3, 4, 5\n",
62+
"$$\n",
63+
"\n",
64+
"\n",
65+
"La política de transición es:\n",
66+
"\n",
67+
"$$\n",
68+
"\\mathbb{P}_{i \\to j} = \n",
69+
"\\left\\{\n",
70+
"\\begin{array}{ll}\n",
71+
"P(D = i - j), & i \\geq 2 , j \\neq 0 \\\\\n",
72+
"P(D \\geq i), & i \\geq 2 , j = 0 \\\\\n",
73+
"P(D = 5 - j), & 1 \\leq i \\lt 2 , j \\neq 0 \\\\\n",
74+
"P(D \\geq 5), & 1 \\leq i \\lt 2 , j = 0 \\\\\n",
75+
"P(Y = j), & i = 0 , j \\in \\{0,1,\\dots,5\\} \\\\\n",
76+
"0, & \\text{d.l.c.}\n",
77+
"\\end{array}\n",
78+
"\\right.\n",
79+
"$$\n",
80+
"\n",
81+
"\n"
82+
]
83+
},
84+
{
85+
"cell_type": "markdown",
86+
"metadata": {},
87+
"source": [
88+
"A continuación, se implementa el modelo en Python."
89+
]
90+
},
1091
{
1192
"cell_type": "code",
1293
"execution_count": null,
13-
"metadata": {
14-
"vscode": {
15-
"languageId": "plaintext"
94+
"metadata": {},
95+
"outputs": [
96+
{
97+
"name": "stdout",
98+
"output_type": "stream",
99+
"text": [
100+
"Matriz de transición P:\n",
101+
"[[0.7 0.26 0.04 0. 0. 0. ]\n",
102+
" [0.18 0.17 0.22 0.22 0.15 0.05]\n",
103+
" [0.8 0.15 0.05 0. 0. 0. ]\n",
104+
" [0.58 0.22 0.15 0.05 0. 0. ]\n",
105+
" [0.35 0.22 0.22 0.15 0.05 0. ]\n",
106+
" [0.18 0.17 0.22 0.22 0.15 0.05]]\n"
107+
]
16108
}
17-
},
109+
],
110+
"source": [
111+
"import numpy as np\n",
112+
"from scipy.stats import poisson, binom\n",
113+
"\n",
114+
"# Parámetros\n",
115+
"demanda_lambda = 3\n",
116+
"hospitales_n = 5\n",
117+
"prob_hospitales = 0.07\n",
118+
"inventario_max = 5\n",
119+
"\n",
120+
"# Espacio de estados\n",
121+
"estados = [i for i in range(0, inventario_max+1)]\n",
122+
"\n",
123+
"# Matriz de transición\n",
124+
"P = np.zeros((len(estados), len(estados)))\n",
125+
"\n",
126+
"# Llenar la matriz de transición\n",
127+
"for i in estados:\n",
128+
" for j in estados:\n",
129+
" if i>=2 and j!=0:\n",
130+
" P[i,j] = poisson.pmf(i - j, demanda_lambda)\n",
131+
" elif i>=2 and j==0:\n",
132+
" P[i,j] = poisson.sf(i-1, demanda_lambda)\n",
133+
" elif 1<=i and i<2 and j!=0:\n",
134+
" P[i,j] = poisson.pmf(inventario_max - j, demanda_lambda)\n",
135+
" elif 1<=i and i<2 and j==0:\n",
136+
" P[i,j] = poisson.sf(inventario_max-1, demanda_lambda)\n",
137+
" elif i==0:\n",
138+
" P[i,j] = binom.pmf(j, hospitales_n, prob_hospitales)\n",
139+
" \n",
140+
"\n",
141+
"print(\"Matriz de transición P:\")\n",
142+
"print(P.round(2))"
143+
]
144+
},
145+
{
146+
"cell_type": "markdown",
147+
"metadata": {},
148+
"source": [
149+
"Para crear la cadena de markov utilizaremos la instancia `dtmc` la libreria `jmarkov`, ya que estamos creando una cadena de tiempo discreto."
150+
]
151+
},
152+
{
153+
"cell_type": "code",
154+
"execution_count": null,
155+
"metadata": {},
18156
"outputs": [],
19-
"source": []
157+
"source": [
158+
"from jmarkov.dtmc import dtmc\n",
159+
"\n",
160+
"# Crear objeto de cadena de Markov\n",
161+
"mc = dtmc(P)"
162+
]
163+
},
164+
{
165+
"cell_type": "markdown",
166+
"metadata": {},
167+
"source": [
168+
"<h4 style=\"color: #ADD8E6;\">Parte B. Responder preguntas de interés</h4>"
169+
]
170+
},
171+
{
172+
"cell_type": "markdown",
173+
"metadata": {},
174+
"source": [
175+
"<h5 style=\"color: #ADD8E6;\">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?</h5>\n",
176+
"\n",
177+
"Sea $X_0 = 1$, queremos calcular:\n",
178+
"\n",
179+
"$$\n",
180+
"P(X_2 = 5 \\mid X_0 = 1)\n",
181+
"$$\n",
182+
"\n",
183+
"Esta probabilidad corresponde a una transición de dos pasos, por lo tanto:\n",
184+
"\n",
185+
"$$\n",
186+
"P(X_2 = 5 \\mid X_0 = 1) = [P^2]_{1,5}\n",
187+
"$$"
188+
]
20189
},
21190
{
22191
"cell_type": "code",
23192
"execution_count": null,
24-
"metadata": {
25-
"vscode": {
26-
"languageId": "plaintext"
193+
"metadata": {},
194+
"outputs": [
195+
{
196+
"name": "stdout",
197+
"output_type": "stream",
198+
"text": [
199+
"La probabilidad de que el inventario sea igual a 5 en 2 días, dado que hoy hay 1 unidad, es: 0.0108\n"
200+
]
27201
}
28-
},
29-
"outputs": [],
30-
"source": []
202+
],
203+
"source": [
204+
"# Probabilidad de estar en estado 5 en dos días si hoy estamos en estado 1\n",
205+
"alpha = np.zeros(len(estados))\n",
206+
"alpha[1] = 1 # Estado inicial = 1\n",
207+
"\n",
208+
"dist_paso2 = mc.transient_probabilities(n=2, alpha=alpha)\n",
209+
"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}\")\n"
210+
]
211+
},
212+
{
213+
"cell_type": "markdown",
214+
"metadata": {},
215+
"source": [
216+
"<h5 style=\"color: #ADD8E6;\">Literal b. ¿Qué pasaría si el umbral de reposición cambiara de 2 a 3? ¿Se reduciría el riesgo de quiebre?</h5>\n",
217+
"\n",
218+
"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.\n",
219+
"\n",
220+
"Quedaría como:\n",
221+
"$$\n",
222+
"p_{ij} = \n",
223+
"\\left\\{\n",
224+
"\\begin{array}{ll}\n",
225+
"P(D = i - j), & i \\geq 3 , j \\neq 0 \\\\\n",
226+
"P(D \\geq i), & i \\geq 3 , j = 0 \\\\\n",
227+
"P(D = 5 - j), & 1 \\leq i \\lt 3 , j \\neq 0 \\\\\n",
228+
"P(D \\geq 5), & 1 \\leq i \\lt 3 , j = 0 \\\\\n",
229+
"P(Y = j), & i = 0 , j \\in \\{0,1,\\dots,5\\} \\\\\n",
230+
"0, & \\text{d.l.c.}\n",
231+
"\\end{array}\n",
232+
"\\right.\n",
233+
"$$\n",
234+
"\n",
235+
"\n",
236+
"Luego de modificar la política, podemos comparar la probabilidad de que el sistema llegue al estado 0 (sin inventario) en el largo plazo.\n",
237+
"\n",
238+
"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.\n",
239+
"\n",
240+
"En este contexto, la **probabilidad de que el sistema esté en quiebre (inventario = 0)** en el largo plazo corresponde al valor $\\pi[0]$. \n",
241+
"\n",
242+
"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.\n"
243+
]
244+
},
245+
{
246+
"cell_type": "code",
247+
"execution_count": null,
248+
"metadata": {},
249+
"outputs": [
250+
{
251+
"name": "stdout",
252+
"output_type": "stream",
253+
"text": [
254+
"π[0] con umbral <2 = 0.5648\n",
255+
"π[0] con umbral <3 = 0.4759\n"
256+
]
257+
}
258+
],
259+
"source": [
260+
"# Matriz de transición ajustando el umbral de reposición de <2 a <3\n",
261+
"P_alt = np.zeros((len(estados), len(estados)))\n",
262+
"\n",
263+
"for i in estados:\n",
264+
" for j in estados:\n",
265+
" if i>=3 and j!=0:\n",
266+
" P_alt[i,j] = poisson.pmf(i - j, demanda_lambda)\n",
267+
" elif i>=3 and j==0:\n",
268+
" P_alt[i,j] = poisson.sf(i-1, demanda_lambda)\n",
269+
" elif 1<=i and i<3 and j!=0:\n",
270+
" P_alt[i,j] = poisson.pmf(inventario_max - j, demanda_lambda)\n",
271+
" elif 1<=i and i<3 and j==0:\n",
272+
" P_alt[i,j] = poisson.sf(inventario_max-1, demanda_lambda)\n",
273+
" elif i==0:\n",
274+
" P_alt[i,j] = binom.pmf(j, hospitales_n, prob_hospitales)\n",
275+
"\n",
276+
"\n",
277+
"# Comparar probabilidad de inventario agotado a largo plazo\n",
278+
"\n",
279+
"mc_old = dtmc(P) # Política con umbral <2\n",
280+
"mc_new = dtmc(P_alt) # Política con umbral <3\n",
281+
"\n",
282+
"pi_old = mc_old.steady_state()\n",
283+
"pi_new = mc_new.steady_state()\n",
284+
"\n",
285+
"print(f\"π[0] con umbral <2 = {pi_old[0]:.4f}\")\n",
286+
"print(f\"π[0] con umbral <3 = {pi_new[0]:.4f}\")"
287+
]
288+
},
289+
{
290+
"cell_type": "markdown",
291+
"metadata": {},
292+
"source": [
293+
"<h5 style=\"color: #ADD8E6;\">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?</h5>\n",
294+
"\n",
295+
"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).\n"
296+
]
297+
},
298+
{
299+
"cell_type": "code",
300+
"execution_count": null,
301+
"metadata": {},
302+
"outputs": [
303+
{
304+
"name": "stdout",
305+
"output_type": "stream",
306+
"text": [
307+
"Comparación:\n",
308+
"Recuperación con p = 0.07: 0.3043\n",
309+
"Recuperación con p = 0.20: 0.6723\n"
310+
]
311+
}
312+
],
313+
"source": [
314+
"# Aumentar la probabilidad de colaboración a 0.2\n",
315+
"nueva_prob_hospitales = 0.2\n",
316+
"\n",
317+
"# Recalcular fila i=0 de la matriz\n",
318+
"P_colab = P.copy()\n",
319+
"for j in estados:\n",
320+
" P_colab[0,j] = binom.pmf(j, hospitales_n, nueva_prob_hospitales)\n",
321+
"\n",
322+
"# Crear cadena con nueva p de colaboración\n",
323+
"mc_alt = dtmc(P_colab)\n",
324+
"\n",
325+
"# Vector inicial desde estado 0\n",
326+
"alpha_cuarentena = np.zeros(len(estados))\n",
327+
"alpha_cuarentena[0] = 1\n",
328+
"\n",
329+
"# Distribución después de 1 paso (reapertura)\n",
330+
"dist_reapertura = mc_alt.transient_probabilities(n=1, alpha=alpha_cuarentena)\n",
331+
"\n",
332+
"# Probabilidad de tener inventario > 0\n",
333+
"prob_recuperacion = 1 - dist_reapertura[0]\n",
334+
"\n",
335+
"# Comparación con la política original\n",
336+
"dist_reapertura_orig = mc.transient_probabilities(n=1, alpha=alpha_cuarentena)\n",
337+
"prob_recuperacion_orig = 1 - dist_reapertura_orig[0]\n",
338+
"\n",
339+
"print(f\"Comparación:\")\n",
340+
"print(f\"Recuperación con p = 0.07: {prob_recuperacion_orig:.4f}\")\n",
341+
"print(f\"Recuperación con p = 0.20: {prob_recuperacion:.4f}\")"
342+
]
31343
},
32344
{
33345
"cell_type": "markdown",

0 commit comments

Comments
 (0)