Distribución Estacionaria de una Cadena de Markov en Tiempo Continuo

En la Teoría de Probabilidades, una Cadena de Markov en Tiempo Continuo (CTMCContinuous Time Markov Chain) corresponde a un modelo matemático donde una variable aleatoria toma valores en un espacio finito y donde el tiempo en que la variable se encuentra en un determinado estado adopta valores reales no negativos que siguen una Distribución Exponencial. Dicha característica denominada Propiedad de No Memoria determina que el comportamiento a futuro del proceso estocástico depende sólo del estado actual y no del comportamiento histórico de dicha variable.

En este contexto a continuación presentamos un ejemplo del cálculo de una Distribución Estacionaria o de Largo Plazo para una Cadena de Markov en Tiempo Continuo. Nuestro objeto de análisis en términos intuitivos es similar al caso de la Distribución Limite de una Cadena de Markov en Tiempo Discreto, es decir, estimar en el largo plazo (e independiente de la distribución inicial) cuál es la probabilidad que el proceso estocástico se encuentre en uno de los N estados posibles.

Distribución Estacionaria Cadena de Markov en Tiempo Continuo

Aves arriban a uno de los 4 alimentadores (comederos) ubicados en el patio de una casa de acuerdo a un Proceso de Poisson con tasa promedio de un ave por minuto. Si todos los alimentadores (4 en total) están ocupados, el ave se va sin esperar que un alimentador se desocupe y en caso contrario el ave come de un alimentador por un tiempo determinado que sigue una Distribución Exponencial con media de un minuto. Se desea modelar la cantidad de alimentadores o comederos ocupados X(t) a través de una Cadena de Markov en Tiempo Continuo.

Sea X(t) la variable aleatoria que representa el número de alimentadores ocupados. En este contexto X(t) representa el número de entidades presentes en los alimentadores en el instante t, cuyo tamaño aumenta con la llegada (nacimiento) de entidades (aves) o disminuye con la salida (muerte) de entidades (aves).

Al existir 4 alimentadores, la cantidad de éstos que se pueden encontrar ocupados en un instante del tiempo son 0, 1, 2, 3 o 4, (X(t)\in\begin{Bmatrix}0,1,2,3,4\end{Bmatrix}t\geq 0) según se representa en el diagrama de transición a continuación:

proceso nacimiento y muerte ctmc

Los valores en las flechas que unen los estados 0 con el 1, 1 con el 2, 2 con el 3 y 3 con el 4, corresponden a la tasa de llegada λ (tasas de nacimiento) que en este caso corresponde a \lambda_{i}=1[\frac{ave}{minuto}] para todo i.

Por otro lado las tasas de servicio (salida o muerte) µ corresponderán a \mu_{i}=i\mu, de modo que de esta forma se obtienen las tasas indicadas en la parte inferior del diagrama en aquella transiciones que consisten en disminuir el número de alimentadores ocupados. Notar que en este caso \mu=1[\frac{ave}{minuto}] .

La cadena propuesta es claramente irreducible (es decir, todos los estados se comunican entre sí y existe por tanto una única clase de estados), de modo que podemos estimar una distribución estacionaria o de largo plazo. En este sentido la ecuación correspondiente es:

ecuaciones largo plazo ctmc

Donde por cierto \sum\pi _{i}=1. La Matriz General G, que tiene como característica que la sumatoria de los valores en sus respectivas filas corresponde a cero, a su vez queda definida por:

matriz g ctmc

De este modo se obtiene el siguiente sistema de ecuaciones:

sistema ecuaciones largo plazo ctmc

La solución del sistema corresponde a:

soluciones ecuaciones ctmc

Que representa para cada estado la proporción del tiempo (en el largo plazo) que cada uno de ellos se encontrará ocupado, por ejemplo, se espera que en el 36,92% del tiempo no se encuentren alimentadores (comederos) ocupados.

Los resultados anteriores han sido corroborados haciendo uso de una planilla de Simulación de Sistemas de Espera disponible en el Libro Investigación de Operaciones de H. Taha. Notar que se ha establecido un límite de 4 entidades (en este caso aves) para el sistema:

simulación ctmc

El número de alimentadores ocupados en el largo plazo son:

alimentadores ocupados largo plazo

Se concluye que en el largo plazo se encuentren ocupados 0,9845 comederos. Cabe destacar que según el resultado de la simulación anterior L_{s}=0,98462 y por cierto la diferencia respecto a nuestro resultado obedece exclusivamente a la cantidad de decimales utilizados para la estimación.

Simulación de una Línea de Espera M/M/1 (Teoría de Colas) en Excel

Un sistema de espera M/M/1 es aquel que considera un servidor, con tiempos exponenciales de servicio y entre llegadas de clientes. La implicancia que los tiempos de servicio se distribuyan exponencial es que existe una preponderancia de tiempos de servicio menores al promedio combinados con algunos pocos tiempos extensos. Un ejemplo de ello es lo que sucede en las cajas de los bancos donde la mayoría de las transacciones requieren poco tiempo de proceso por parte del cajero, no obstante algunas transacciones más complejas consumen bastante tiempo. Por otra parte afirmar que los tiempos entre llegadas se distribuyen exponencial implica una preponderancia de tiempos entre llegadas menores que el promedio en combinación con algunos tiempos más extensos. Lo anterior tiene relación con la aleatoriedad del proceso de llegada de clientes que permite establecer la Propiedad de Falta de Memoria o Amnesia de la Distribución Exponencial y con los conceptos presentados en el artículo Qué son las Líneas de Espera (Teoría de Colas), donde queda en evidencia que la formación de las colas o filas esta asociada a la variabilidad del sistema.

En este contexto consideremos la siguiente notación, donde valores usuales para A y B son M (distribución exponencial) y G (distribución general).

notacion-de-kendall

El siguiente ejemplo disponible en el artículo Qué es la Ley de Little y su aplicación en el análisis de Líneas de Espera, nos permitirá ilustrar la simulación en Excel del comportamiento de un sistema de espera M/M/1.

Simulación de una Línea de Espera M/M/1

Ejemplo: Un pequeño banco está considerando abrir un servicio para que los clientes paguen desde su automóvil. Se estima que los clientes llegarán a una tasa promedio de 15 por hora (λ=15). El cajero que trabajará en la ventanilla puede atender a los clientes a un ritmo promedio de uno cada tres minutos (es decir, la capacidad promedio del servidor es de µ=20). Suponga que el patrón de llegadas es Poisson y el patrón de servicios es Exponencial.

Al hacer uso de la Planilla de Fórmulas de Sistema de Espera, se alcanzan los resultados que se resumen en la imagen a continuación.

salida-planilla-linea-de-es

Con esto la utilización promedio del servidor es de un 75%, el número esperado de clientes en el sistema Ls es 3, el número esperado de clientes en la fila Lq son 2,5, el tiempo promedio que un cliente permanece en el sistema Ws (espera más atención) son 0,2 horas (0 12 minutos) y el tiempo promedio que un cliente permanece en la fila Wq (esperando su atención) es de 0,15 horas (o 9 minutos).

Otra alternativa es simular el comportamiento de la línea de espera con configuración M/M/1 haciendo uso de Excel. Para ello ingresamos en la planilla Queueing_Simulator el número de servidores (1), la distribución para el tiempo entre llegadas (exponencial con media de 4 minutos, esto es, si llegan en promedio 15 clientes por hora, entonces en promedio llega un cliente cada 1/15 de hora o equivalentemente cada 4 minutos) y una distribución para el tiempo de servicio también exponencial con media de 3 minutos. Finalmente ingresamos el número de llegadas que se desea simular (arbitrariamente se ha seleccionado 100.000 llegadas para evaluar un comportamiento del sistema en el largo plazo) y luego Run Simulation.

simulacion-mm1-excel

Se puede apreciar que los resultados obtenidos en la columna F son (aproximadamente) similares a los obtenidos utilizando los resultados que establece la Ley de Little. Por ejemplo, el número esperado de clientes en el sistema L es 3,0157; el número esperado de clientes en la fila Lq es 2,2665; el tiempo esperado que un cliente permanece en el sistema W son 12,0612 minutos y así sucesivamente.

Importante: Los resultados mostrados anteriormente corresponden a aquellos obtenidos con una simulación tipo. Si una vez alcanzados los resultados presionamos nuevamente Run Simulation obtendremos cambios en los resultados los cuales de todos modos deberían aproximar los resultados de la Ley de Little bajo el supuesto de considerar un número significativo de llegadas a simular.

¿Quieres tener el archivo Excel con la Simulación de una Línea de Espera M/M/1 utilizada en este ejemplo?

[sociallocker]

MUCHAS GRACIAS!. DESCARGA AQUÍ EL ARCHIVO

[/sociallocker]

Cálculo de la Probabilidad de un Número de Llegadas en un Tiempo Determinado utilizando la Distribución de Poisson

Cuando los clientes llegan a un servicio de forma totalmente aleatoria (es decir, no hay forma de pronosticar cuándo va a llegar alguien) la función de densidad de probabilidad para describir la cantidad de llegadas durante un tiempo determinado se representa por la Distribución de Poisson y automáticamente la distribución del tiempo entre llegadas sigue una Distribución Exponencial según lo expuesto en el artículo Propiedad de Falta de Memoria o Amnesia de la Distribución Exponencial.

En este contexto la fórmula que permite calcular la probabilidad exacta de n llegadas dentro de un período T es la siguiente:

probabilidad-poisson

Consideremos por ejemplo un taller que se dedica a labores de reparación y que la llegada de éstos diariamente se comporta de forma aleatoria con una tasa de 10 trabajos diarios. ¿Cuál es la probabilidad de que no lleguen trabajos durante una hora cualquiera bajo el supuesto que el taller opera 8 horas al día?.

probabilida-cero-llegadas-p

Notar que \lambda =\frac{10}{8}=1,25[\frac{trabajos}{hora}]. Es decir, la probabilidad de no recibir trabajos durante una hora cualquiera es aproximadamente a un 28,65%.

Asumamos ahora una nueva situación. Un proceso que tiene una tasa promedio de llegada de 6 clientes por hora (\lambda =6[clientes/hora]) y se desea evaluar cuál es la probabilidad de que lleguen exactamente 0, 1, 2,…,n clientes en un intervalo de tiempo de 0,5 horas (30 minutos). El siguiente vídeo proporciona una simulación de dicho escenario:

En el gráfico, el área amarilla, por ejemplo, significa exactamente la probabilidad que 3 personas lleguen en las 0,5 horas. El área amarilla más el área roja, por ejemplo, significa la probabilidad de que lleguen 2 o 3 personas en los 30 minutos.

Adicionalmente haciendo uso del software Geogebra y su herramienta cálculos de probabilidad, se puede representar la Distribución de Poisson para los parámetros descritos anteriormente de forma de obtener rápidamente los resultados para distintos números de llegadas (notar que la Distribución de Poisson es discreta).

distribucion-poisson-geogeb

Cálculo de la Probabilidad de Aceptación de Lotes Productivos en la Evaluación de Proveedores

Los criterios para la evaluación de proveedores son múltiples y en particular los procedimientos que se utilizan para discriminar si un lote productivo se acepta o rechaza son críticos para garantizar la calidad de los insumos sobre los cuales se desarrollará un proceso productivo. En el siguiente artículo se presenta un ejemplo que consiste en la comparación de 2 proveedores en términos de la probabilidad de aceptación de lotes para distintos planes de muestreo, en conjunto con una estimación de la calidad promedio a la salida AOQ (Average Outgoing Quality) luego de la rectificación de las unidades no conformantes en cada caso.

Cálculo de la Probabilidad de Aceptación de Lotes Productivos

Usted está considerando la evaluación de 2 proveedores que le ofrecen uno de los principales insumos para su proceso productivo. El proveedor A, entrega lotes de 2.000 unidades, con un número promedio histórico de defectuosos por lote de 50 unidades. El proveedor B, entrega lotes de 4.000 unidades, teniendo un número histórico de defectuosos de 120 unidades por lote.

Al primer proveedor se le aplica un plan de muestreo simple de n=30 y c=2, es decir, se toma una muestra aleatoria de 30 unidades y se decide aceptar el lote si el número de unidades defectuosas detectadas es menor o igual a 2.

Al segundo proveedor se le aplica un plan de muestreo según norma MIL-STD-105, nivel riguroso y con una calidad media aceptable, AQL de 1%.

Se sabe que el costo de inspección es de $100 por unidad, independiente del proveedor que se tenga, mientras que el costo por artículo fallado que pasa la inspección es de $1.500, el cual está compuesto principalmente por garantías efectivas, las cuales no implican la reposición del producto sino su reparación.

Determine la probabilidad de aceptación del lote para ambos planes de muestreo, dados los tamaños de muestra elegidos, además de establecer las diferencias generadas en términos de calidad promedio de salida del lote.

En relación a la información anterior se tiene que para el Proveedor A el porcentaje promedio de defectuosos según el registro histórico es p=50/2.000=2,5% y para el Proveedor B es p=120/4.000=3,0%.

proveedores-inspeccion

Para determinar la probabilidad de aceptación de un lote proveniente del Proveedor A utilizamos la Distribución de Poisson dado que se cumplen los requisitos para su aplicación, a saber: n>15 ; p<10% ; N>10*n. Luego podemos hacer uso de Excel para obtener la probabilidad de aceptación del lote dada por la fórmula =POISSON(2;0,75;VERDADERO), donde c=2 corresponde al número de aceptación y 0,75 a n*p (30*0,025).

formula-poisson-excel

En el caso del Proveedor B se utiliza Military Standard 105D (conocido también por MIL-STD-105) donde el plan de muestreo esta dado por:

muestreo-mil-std-105

Notar que se considera nivel riguroso, con un tamaño de lote en el intervalo entre [3.201,10.000] unidades y AQL de 1%. El código asociado a este plan de muestreo es M que determina un tamaño de muestral de 315 unidades y un número de aceptación c=7. Esto se puede corroborar con algunas herramientas online para el cálculo del tamaño del lote utilizando la norma MIL-STD-105 como se observa a continuación:

tamaño-muestra-mil-std

En resumen la probabilidad de aceptación para ambos proveedores es:

plan-de-muestreo

Además:

tabla-calculo-aoq

Ambos valores calculados como Pa*p en vista de que N>>n. Por lo tanto, la inspección realizada al Proveedor B tendrá una mejor calidad promedio a la salida luego de la rectificación. Lo anterior se debe principalmente a la baja probabilidad de aceptación de estos lotes debido a un muestreo más riguroso, lo cual aumenta la calidad de los lotes efectivamente recibidos.

Cómo construir una Curva Característica de Operación (CO) con Excel

La Curva Característica de Operación (o Curva Característica Operativa) consiste en una representación gráfica que muestra para un plan de muestreo específico (n,c) la probabilidad de aceptación del lote, para varios valores de calidad del lote a la entrada p (% de unidades defectuosas). Una Curva Característica de Operación tendrá entre sus puntos uno definido por el NCA y 1-α y otro punto definido por PTDL y β.

Para determinar la probabilidad de aceptación de un lote, se pueden utilizar las distribuciones Binomial o Poisson. Cuando el tamaño de la muestra es al menos 15 unidades (n>15), el tamaño del lote es al menos 10 veces el tamaño de la muestra (N>10*n) y el porcentaje de unidades defectuosas históricamente es menor a un 10% (p<10%), es preferible la Distribución de Poisson debido a la facilidad de los cálculos.

Consideremos el siguiente ejemplo: Un fabricante de pistones para motocicletas comenzará a vender diariamente 1.200 unidades para un nuevo cliente. Este último determina condiciones contractuales para la inspección del lote diario, especificando que tomará muestras de 100 unidades (n=100) y que sólo aceptará los pedidos con 4 o menos defectos (c=4). El fabricante menciona en el contrato, que históricamente ha obtenido un porcentaje defectivo del 2% (p=2%). Determinar la probabilidad de aceptación del lote por parte del cliente.

Para determinar la probabilidad de aceptación del lote podemos utilizar la siguiente fórmula haciendo uso de una planilla de cálculo Excel: =POISSON(4;100*2%;VERDADERO). En consecuencia la probabilidad de aceptación del lote por parte del cliente es de un 94,73% (aproximado).

poisson-probabilidad-acepta

Una alternativa a Excel para estos efectos es hacer uso de la herramienta de cálculo de probabilidades del software Geogebra donde se debe ingresar los parámetros de la distribución µ (equivalente a n*p=100*2%=2) y el valor correspondiente al número de aceptación c (en la imagen c=4).

poisson-geogebra

De esta forma se puede extender el procedimiento calculando la probabilidad de aceptación del lote Pa para distintos valores de calidades a la entrada p. Un extracto de ello se presenta en la siguiente tabla:

extracto-tabla-calculo-prob

A continuación se construye un gráfico con la Curva Característica de Operación. Se ha destacado con una etiqueta color amarillo el dato que hemos calculado previamente.

curva-caracteristica-operat

En un próximo artículo discutiremos el impacto que tiene en el plan de muestreo y en particular en la Curva Característica de Operación un cambio en el tamaño de la muestra o un cambio en el número de aceptación.