Simulación de Monte Carlo en R

Estoy tratando de resolver el siguiente ejercicio, pero en realidad no tengo idea de cómo comenzar a hacerlo. Encontré un código en mi libro que parece, pero es un ejercicio completamente diferente y no sé cómo relacionarlos entre sí. ¿Cómo puedo comenzar a simular llegadas y cómo sé cuándo están terminadas? Sé cómo almacenarlos y calcular a, b, c, d de acuerdo con eso. Pero no sé cómo realmente necesito simular la simulación de Monte Carlo. ¿Podría alguien ayudarme a comenzar? Sé que este no es un lugar donde sus preguntas son respondidas por usted, sino que solo se resuelven. Pero el problema es que no sé cómo empezar.

Una mesa de ayuda de soporte de TI representa un sistema de colas con cinco asistentes que reciben llamadas de los clientes. Las llamadas se realizan de acuerdo con un proceso de Poisson con la tasa promedio de una llamada cada 45 segundos. Los tiempos de servicio para los asistentes primero, segundo, tercero, cuarto y quinto son variables aleatorias exponenciales con parámetros λ1 = 0.1, λ2 = 0.2, λ3 = 0.3, λ4 = 0.4 y λ5 = 0.5 min − 1, respectivamente (el jth help desk assistant tiene λk = k / 10 min − 1). Además de los clientes que reciben asistencia, se puede poner en espera a otros diez clientes. En los momentos en que se alcanza esta capacidad, las nuevas personas que llaman reciben una señal de ocupado. Utilice los métodos de Monte Carlo para estimar las siguientes características de rendimiento:

(a) la fracción de clientes que reciben una señal de ocupado;

(b) el tiempo de respuesta esperado;

(c) el tiempo de espera promedio;

(d) la porción de clientes atendidos por cada asistente de mesa de ayuda;

EDITAR: lo que tengo hasta ahora es (no mucho):

pa = 1/45sec-1

jobs = rep(1,5); onHold = rep(1,10);

jobsIndex = 0;

onHoldIndex = 0;

u = runif(1)
for (i in 1:1000) {

    if(u  <= pa){ # new arrival

        if(jobsIndex < 5) # assistant is free, #give job to assistant

            jobsIndex++;

        else #add to onHold array

            onHoldIndex++;
    }
}

Este es uno de los tipos de simulación más instructivos y divertidos para realizar: crea agentes independientes en la computadora, los deja interactuar, realiza un seguimiento de lo que hacen y estudia lo que sucede. Es una forma maravillosa de aprender sobre sistemas complejos, especialmente (pero no limitado a) aquellos que no se pueden entender con un análisis puramente matemático.

La mejor manera de construir tales simulaciones es con un diseño de arriba hacia abajo.

En el nivel más alto, el código debería ser similar a

initialize(...)
while (process(get.next.event())) {}

(Este y todos los ejemplos posteriores son código ejecutable R , no solo pseudocódigo). El bucle es una simulación dirigida por eventos : get.next.event()encuentra cualquier "evento" de interés y le pasa una descripción process, que hace algo con él (incluido el registro de cualquier información al respecto). Regresa TRUEmientras las cosas funcionen bien; Al identificar un error o el final de la simulación, regresa FALSEy finaliza el ciclo.

Si imaginamos una implementación física de esta cola, como personas que esperan una licencia de matrimonio en la ciudad de Nueva York o una licencia de conducir o un boleto de tren en casi cualquier lugar, pensamos en dos tipos de agentes: clientes y "asistentes" (o servidores) . Los clientes se anuncian apareciendo; los asistentes anuncian su disponibilidad encendiendo una luz o señal o abriendo una ventana. Estos son los dos tipos de eventos para procesar.

El entorno ideal para tal simulación es una verdadera orientación a objetos en la que los objetos son mutables : pueden cambiar de estado para responder de forma independiente a las cosas que los rodean. Res absolutamente terrible para esto (¡incluso Fortran sería mejor!). Sin embargo, aún podemos usarlo si tenemos cuidado. El truco consiste en mantener toda la información en un conjunto común de estructuras de datos a las que se puede acceder (y modificar) mediante muchos procedimientos separados e interactivos. Adoptaré la convención de usar nombres de variables EN TODAS LAS MAYÚSCULAS para dichos datos.

El siguiente nivel del diseño de arriba hacia abajo es codificar process. Responde a un solo descriptor de eventos e:

process <- function(e) {
  if (is.null(e)) return(FALSE)
  if (e$type == "Customer") {
    i <- find.assistant(e$time)
    if (is.null(i)) put.on.hold(e$x, e$time) else serve(i, e$x, e$time)
  } else {
    release.hold(e$time)
  }
  return(TRUE)
}

Tiene que responder a un evento nulo cuando get.next.eventno tiene eventos para informar. De lo contrario, processimplementa las "reglas de negocio" del sistema. Prácticamente se escribe a partir de la descripción en la pregunta. Su funcionamiento debería requerir pocos comentarios, excepto para señalar que eventualmente tendremos que codificar subrutinas put.on.holde release.hold(implementar una cola de retención de clientes) e serve(implementar las interacciones cliente-asistente).

¿Qué es un "evento"? Debe contener información sobre quién está actuando, qué tipo de acción están tomando y cuándo está ocurriendo. Por lo tanto, mi código usa una lista que contiene estos tres tipos de información. Sin embargo, get.next.eventsolo necesita inspeccionar los tiempos. Es responsable solo de mantener una cola de eventos en los que

1. Cualquier evento se puede poner en la cola cuando se recibe y

2. El primer evento en la cola se puede extraer y pasar fácilmente a la persona que llama.

La mejor implementación de esta cola prioritaria sería un montón, pero eso es demasiado exigente R. Siguiendo una sugerencia en The Art of R Programming de Norman Matloff (que ofrece un simulador de cola más flexible, abstracto pero limitado), he usado un marco de datos para contener los eventos y simplemente buscarlo por el tiempo mínimo entre sus registros.

get.next.event <- function() {
  if (length(EVENTS$time) <= 0) new.customer()               # Wait for a customer$
  if (length(EVENTS$time) <= 0) return(NULL)                 # Nothing's going on!$
  if (min(EVENTS$time) > next.customer.time()) new.customer()# See text
  i <- which.min(EVENTS$time)
  e <- EVENTS[i, ]; EVENTS <<- EVENTS[-i, ]
  return (e)
}

Hay muchas formas en que esto podría haberse codificado. La versión final que se muestra aquí refleja una elección que hice al codificar cómo processreacciona a un evento "Asistente" y cómo new.customerfunciona: get.next.eventsimplemente saca a un cliente de la cola de espera, luego se sienta y espera otro evento. A veces será necesario buscar un nuevo cliente de dos maneras: primero, para ver si uno está esperando en la puerta (por así decirlo) y segundo, si uno ha entrado cuando no estábamos buscando.

Claramente, new.customery next.customer.timeson rutinas importantes , así que cuidémoslas a continuación.

new.customer <- function() {  
  if (CUSTOMER.COUNT < dim(CUSTOMERS)[2]) {
    CUSTOMER.COUNT <<- CUSTOMER.COUNT + 1
    insert.event(CUSTOMER.COUNT, "Customer", 
                 CUSTOMERS["Arrived", CUSTOMER.COUNT])
  }
  return(CUSTOMER.COUNT)
}
next.customer.time <- function() {
  if (CUSTOMER.COUNT < dim(CUSTOMERS)[2]) {
    x <- CUSTOMERS["Arrived", CUSTOMER.COUNT]
  } else {x <- Inf}
  return(x) # Time when the next customer will arrive
}

CUSTOMERSes una matriz 2D, con datos para cada cliente en columnas. Tiene cuatro filas (que actúan como campos) que describen a los clientes y registran sus experiencias durante la simulación : "Llegado", "Servido", "Duración" y "Asistente" (un identificador numérico positivo del asistente, si lo hay, que sirvió ellos, y de lo contrario -1para señales de ocupado). En una simulación altamente flexible, estas columnas se generarían dinámicamente, pero debido a cómo Rle gusta trabajar, es conveniente generar todos los clientes desde el principio, en una única matriz grande, con sus tiempos de llegada ya generados al azar. next.customer.timepuede echar un vistazo a la siguiente columna de esta matriz para ver quién vendrá después. La variable globalCUSTOMER.COUNTindica el último cliente en llegar. Los clientes se gestionan de manera muy simple mediante este puntero, avanzando para obtener un nuevo cliente y mirando más allá (sin avanzar) para echar un vistazo al próximo cliente.

serve implementa las reglas de negocio en la simulación.

serve <- function(i, x, time.now) {
  #
  # Serve customer `x` with assistant `i`.
  #
  a <- ASSISTANTS[i, ]
  r <- rexp(1, a$rate)                       # Simulate the duration of service
  r <- round(r, 2)                           # (Make simple numbers)
  ASSISTANTS[i, ]$available <<- time.now + r # Update availability
  #
  # Log this successful service event for later analysis.
  #
  CUSTOMERS["Assistant", x] <<- i
  CUSTOMERS["Served", x] <<- time.now
  CUSTOMERS["Duration", x] <<- r
  #
  # Queue the moment the assistant becomes free, so they can check for
  # any customers on hold.
  #
  insert.event(i, "Assistant", time.now + r)
  if (VERBOSE) cat(time.now, ": Assistant", i, "is now serving customer", 
                   x, "until", time.now + r, "\n")
  return (TRUE)
}

Esto es sencillo. ASSISTANTSes un marco de datos con dos campos: capabilities(indicando su tasa de servicio) y available, que marca la próxima vez que el asistente estará libre. Se atiende a un cliente generando una duración de servicio aleatoria de acuerdo con las capacidades del asistente, actualizando el tiempo cuando el asistente esté disponible y registrando el intervalo de servicio en la CUSTOMERSestructura de datos. El VERBOSEindicador es útil para probar y depurar: cuando es verdadero, emite una secuencia de oraciones en inglés que describe los puntos clave de procesamiento.

Cómo se asignan los asistentes a los clientes es importante e interesante. Uno puede imaginar varios procedimientos: asignación al azar, por algún orden fijo, o de acuerdo con quién ha estado libre el tiempo más largo (o más corto). Muchos de estos se ilustran en código comentado:

find.assistant <- function(time.now) {
  j <- which(ASSISTANTS$available <= time.now)
  #if (length(j) > 0) {
  #  i <- j[ceiling(runif(1) * length(j))]
  #} else i <- NULL                                    # Random selection
  #if (length(j) > 0) i <- j[1] else i <- NULL         # Pick first assistant
  #if (length(j) > 0) i <- j[length(j)] else i <- NULL # Pick last assistant
  if (length(j) > 0) {
    i <- j[which.min(ASSISTANTS[j, ]$available)]
  } else i <- NULL                                     # Pick most-rested assistant
  return (i)
}

El resto de la simulación es realmente solo un ejercicio de rutina para persuadir Ra implementar estructuras de datos estándar, principalmente un búfer circular para la cola en espera. Debido a que no quieres correr sin control con los globales, puse todo esto en un solo procedimiento sim. Sus argumentos describen el problema: el número de clientes a simular ( n.events), la tasa de llegada de clientes, las capacidades de los asistentes y el tamaño de la cola de espera (que puede establecerse en cero para eliminar la cola por completo).

r <- sim(n.events=250, arrival.rate=60/45, capabilities=1:5/10, hold.queue.size=10)

Devuelve una lista de las estructuras de datos mantenidas durante la simulación; el de mayor interés es la CUSTOMERSmatriz. Rhace que sea bastante fácil trazar la información esencial en esta matriz de una manera interesante. Aquí hay una salida que muestra los últimos clientes en una simulación más larga de clientes.$50$$250$

La experiencia de cada cliente se traza como una línea de tiempo horizontal, con un símbolo circular al momento de la llegada, una línea negra sólida para cualquier espera en espera y una línea de color durante la interacción con un asistente (el color y el tipo de línea diferenciar entre los asistentes). Debajo de esta trama de Clientes hay una que muestra las experiencias de los asistentes, marcando los momentos en que estaban y no estaban comprometidos con un cliente. Los puntos finales de cada intervalo de actividad están delimitados por barras verticales.

Cuando se ejecuta con verbose=TRUE, la salida de texto de la simulación se ve así:

...
160.71 : Customer 211 put on hold at position 1 
161.88 : Customer 212 put on hold at position 2 
161.91 : Assistant 3 is now serving customer 213 until 163.24 
161.91 : Customer 211 put on hold at position 2 
162.68 : Assistant 4 is now serving customer 212 until 164.79 
162.71 : Assistant 5 is now serving customer 211 until 162.9 
163.51 : Assistant 5 is now serving customer 214 until 164.05 
...

(Los números a la izquierda son las horas en que se emitió cada mensaje). Puede hacer coincidir estas descripciones con las partes del diagrama de Clientes que se encuentran entre los tiempos y .$160$$165$

Podemos estudiar la experiencia de los clientes en espera trazando las duraciones en espera por identificador de cliente, usando un símbolo especial (rojo) para mostrar a los clientes que reciben una señal de ocupado.

(¡No todas estas parcelas serían un maravilloso tablero en tiempo real para cualquiera que administre esta cola de servicio!)

Es fascinante comparar los gráficos y las estadísticas que obtienes al variar los parámetros que se pasan sim. ¿Qué sucede cuando los clientes llegan demasiado rápido para ser procesados? ¿Qué sucede cuando la cola de espera se hace más pequeña o se elimina? ¿Qué cambia cuando los asistentes son seleccionados de diferentes maneras? ¿Cómo influyen los números y las capacidades de los asistentes en la experiencia del cliente? ¿Cuáles son los puntos críticos en los que algunos clientes comienzan a ser rechazados o quedan en espera por mucho tiempo?

Normalmente, para preguntas evidentes de autoaprendizaje como esta, nos detendríamos aquí y dejaríamos los detalles restantes como un ejercicio. Sin embargo, no quiero decepcionar a los lectores que pueden haber llegado tan lejos y están interesados ​​en probar esto por sí mismos (y tal vez modificarlo y desarrollarlo para otros fines), así que a continuación se adjunta el código de trabajo completo.

(El procesamiento de en este sitio estropeará la sangría en cualquier línea que contenga un símbolo , pero la sangría legible debe restaurarse cuando el código se pega en un archivo de texto).$\TeX$$\$$

sim <- function(n.events, verbose=FALSE, ...) {
  #
  # Simulate service for `n.events` customers.
  #
  # Variables global to this simulation (but local to the function):
  #
  VERBOSE <- verbose         # When TRUE, issues informative message
  ASSISTANTS <- list()       # List of assistant data structures
  CUSTOMERS <- numeric(0)    # Array of customers that arrived
  CUSTOMER.COUNT <- 0        # Number of customers processed
  EVENTS <- list()           # Dynamic event queue   
  HOLD <- list()             # Customer on-hold queue
  #............................................................................#
  #
  # Start.
  #
  initialize <- function(arrival.rate, capabilities, hold.queue.size) {
    #
    # Create common data structures.
    #
    ASSISTANTS <<- data.frame(rate=capabilities,     # Service rate
                              available=0            # Next available time
    )
    CUSTOMERS <<- matrix(NA, nrow=4, ncol=n.events, 
                         dimnames=list(c("Arrived",  # Time arrived
                                         "Served",   # Time served
                                         "Duration", # Duration of service
                                         "Assistant" # Assistant id
                         )))
    EVENTS <<- data.frame(x=integer(0),              # Assistant or customer id
                          type=character(0),         # Assistant or customer
                          time=numeric(0)            # Start of event
    )
    HOLD <<- list(first=1,                           # Index of first in queue
                  last=1,                            # Next available slot
                  customers=rep(NA, hold.queue.size+1))
    #
    # Generate all customer arrival times in advance.
    #
    CUSTOMERS["Arrived", ] <<- cumsum(round(rexp(n.events, arrival.rate), 2))
    CUSTOMER.COUNT <<- 0
    if (VERBOSE) cat("Started.\n")
    return(TRUE)
  }
  #............................................................................#
  #
  # Dispatching.
  #
  # Argument `e` represents an event, consisting of an assistant/customer 
  # identifier `x`, an event type `type`, and its time of occurrence `time`.
  #
  # Depending on the event, a customer is either served or an attempt is made
  # to put them on hold.
  #
  # Returns TRUE until no more events occur.
  #
  process <- function(e) {
    if (is.null(e)) return(FALSE)
    if (e$type == "Customer") {
      i <- find.assistant(e$time)
      if (is.null(i)) put.on.hold(e$x, e$time) else serve(i, e$x, e$time)
    } else {
      release.hold(e$time)
    }
    return(TRUE)
  }#$
  #............................................................................#
  #
  # Event queuing.
  #
  get.next.event <- function() {
    if (length(EVENTS$time) <= 0) new.customer()
    if (length(EVENTS$time) <= 0) return(NULL)
    if (min(EVENTS$time) > next.customer.time()) new.customer()
    i <- which.min(EVENTS$time)
    e <- EVENTS[i, ]; EVENTS <<- EVENTS[-i, ]
    return (e)
  }
  insert.event <- function(x, type, time.occurs) {
    EVENTS <<- rbind(EVENTS, data.frame(x=x, type=type, time=time.occurs))
    return (NULL)
  }
  # 
  # Customer arrivals (called by `get.next.event`).
  #
  # Updates the customers pointer `CUSTOMER.COUNT` and returns the customer
  # it newly points to.
  #
  new.customer <- function() {  
    if (CUSTOMER.COUNT < dim(CUSTOMERS)[2]) {
      CUSTOMER.COUNT <<- CUSTOMER.COUNT + 1
      insert.event(CUSTOMER.COUNT, "Customer", 
                   CUSTOMERS["Arrived", CUSTOMER.COUNT])
    }
    return(CUSTOMER.COUNT)
  }
  next.customer.time <- function() {
    if (CUSTOMER.COUNT < dim(CUSTOMERS)[2]) {
      x <- CUSTOMERS["Arrived", CUSTOMER.COUNT]
    } else {x <- Inf}
    return(x) # Time when the next customer will arrive
  }
  #............................................................................#
  #
  # Service.
  #
  find.assistant <- function(time.now) {
    #
    # Select among available assistants.
    #
    j <- which(ASSISTANTS$available <= time.now) 
    #if (length(j) > 0) {
    #  i <- j[ceiling(runif(1) * length(j))]
    #} else i <- NULL                                    # Random selection
    #if (length(j) > 0) i <- j[1] else i <- NULL         # Pick first assistant
    #if (length(j) > 0) i <- j[length(j)] else i <- NULL # Pick last assistant
    if (length(j) > 0) {
      i <- j[which.min(ASSISTANTS[j, ]$available)]
    } else i <- NULL # Pick most-rested assistant
    return (i)
  }#$
  serve <- function(i, x, time.now) {
    #
    # Serve customer `x` with assistant `i`.
    #
    a <- ASSISTANTS[i, ]
    r <- rexp(1, a$rate)                       # Simulate the duration of service
    r <- round(r, 2)                           # (Make simple numbers)
    ASSISTANTS[i, ]$available <<- time.now + r # Update availability
    #
    # Log this successful service event for later analysis.
    #
    CUSTOMERS["Assistant", x] <<- i
    CUSTOMERS["Served", x] <<- time.now
    CUSTOMERS["Duration", x] <<- r
    #
    # Queue the moment the assistant becomes free, so they can check for
    # any customers on hold.
    #
    insert.event(i, "Assistant", time.now + r)
    if (VERBOSE) cat(time.now, ": Assistant", i, "is now serving customer", 
                     x, "until", time.now + r, "\n")
    return (TRUE)
  }
  #............................................................................#
  #
  # The on-hold queue.
  #
  # This is a cicular buffer implemented by an array and two pointers,
  # one to its head and the other to the next available slot.
  #
  put.on.hold <- function(x, time.now) {
    #
    # Try to put customer `x` on hold.
    #
    if (length(HOLD$customers) < 1 || 
          (HOLD$first - HOLD$last %% length(HOLD$customers) == 1)) {
      # Hold queue is full, alas.  Log this occurrence for later analysis.
      CUSTOMERS["Assistant", x] <<- -1 # Busy signal
      if (VERBOSE) cat(time.now, ": Customer", x, "got a busy signal.\n")
      return(FALSE)
    }
    #
    # Add the customer to the hold queue.
    #
    HOLD$customers[HOLD$last] <<- x
    HOLD$last <<- HOLD$last %% length(HOLD$customers) + 1
    if (VERBOSE) cat(time.now, ": Customer", x, "put on hold at position", 
                 (HOLD$last - HOLD$first - 1) %% length(HOLD$customers) + 1, "\n")
    return (TRUE)
  }
  release.hold <- function(time.now) {
    #
    # Pick up the next customer from the hold queue and place them into
    # the event queue.
    #
    if (HOLD$first != HOLD$last) {
      x <- HOLD$customers[HOLD$first]   # Take the first customer
      HOLD$customers[HOLD$first] <<- NA # Update the hold queue
      HOLD$first <<- HOLD$first %% length(HOLD$customers) + 1
      insert.event(x, "Customer", time.now)
    }
  }$
  #............................................................................#
  #
  # Summaries.
  #
  # The CUSTOMERS array contains full information about the customer experiences:
  # when they arrived, when they were served, how long the service took, and
  # which assistant served them.
  #
  summarize <- function() return (list(c=CUSTOMERS, a=ASSISTANTS, e=EVENTS,
                                       h=HOLD))
  #............................................................................#
  #
  # The main event loop.
  #
  initialize(...)
  while (process(get.next.event())) {}
  #
  # Return the results.
  #
  return (summarize())
}
#------------------------------------------------------------------------------#
#
# Specify and run a simulation.
#
set.seed(17)
n.skip <- 200  # Number of initial events to skip in subsequent summaries
system.time({
  r <- sim(n.events=50+n.skip, verbose=TRUE, 
           arrival.rate=60/45, capabilities=1:5/10, hold.queue.size=10)
})
#------------------------------------------------------------------------------#
#
# Post processing.
#
# Skip the initial phase before equilibrium.
#
results <- r$c
ids <- (n.skip+1):(dim(results)[2])
arrived <- results["Arrived", ]
served <- results["Served", ]
duration <- results["Duration", ]
assistant <- results["Assistant", ]
assistant[is.na(assistant)] <- 0   # Was on hold forever
ended <- served + duration
#
# A detailed plot of customer experiences.
#
n.events <- length(ids)
n.assistants <- max(assistant, na.rm=TRUE) 
colors <- rainbow(n.assistants + 2)
assistant.color <- colors[assistant + 2]
x.max <- max(results["Served", ids] + results["Duration", ids], na.rm=TRUE)
x.min <- max(min(results["Arrived", ids], na.rm=TRUE) - 2, 0)
#
# Lay out the graphics.
#
layout(matrix(c(1,1,2,2), 2, 2, byrow=TRUE), heights=c(2,1))
#
# Set up the customers plot.
#
plot(c(x.min, x.max), range(ids), type="n",
     xlab="Time", ylab="Customer Id", main="Customers")
#
# Place points at customer arrival times.
#
points(arrived[ids], ids, pch=21, bg=assistant.color[ids], col="#00000070")
#
# Show wait times on hold.
#
invisible(sapply(ids, function(i) {
  if (!is.na(served[i])) lines(x=c(arrived[i], served[i]), y=c(i,i))
}))
#
# More clearly show customers getting a busy signal.
#
ids.not.served <- ids[is.na(served[ids])]
ids.served <- ids[!is.na(served[ids])]
points(arrived[ids.not.served], ids.not.served, pch=4, cex=1.2)
#
# Show times of service, colored by assistant id.
#
invisible(sapply(ids.served, function(i) {
  lines(x=c(served[i], ended[i]), y=c(i,i), col=assistant.color[i], lty=assistant[i])
}))
#
# Plot the histories of the assistants.
#
plot(c(x.min, x.max), c(1, n.assistants)+c(-1,1)/2, type="n", bty="n",
     xlab="", ylab="Assistant Id", main="Assistants")
abline(h=1:n.assistants, col="#808080", lwd=1)
invisible(sapply(1:(dim(results)[2]), function(i) {
  a <- assistant[i]
  if (a > 0) {
    lines(x=c(served[i], ended[i]), y=c(a, a), lwd=3, col=colors[a+2])
    points(x=c(served[i], ended[i]), y=c(a, a), pch="|", col=colors[a+2])
  }
}))
#
# Plot the customer waiting statistics.
#
par(mfrow=c(1,1))
i <- is.na(served)
plot(served - arrived, xlab="Customer Id", ylab="Minutes",
     main="Service Wait Durations")
lines(served - arrived, col="Gray")
points(which(i), rep(0, sum(i)), pch=16, col="Red")
#
# Summary statistics.
#
mean(!is.na(served)) # Proportion of customers served
table(assistant)