Simulation de Monte Carlo en R

11

J'essaie de résoudre l'exercice suivant, mais je n'ai en fait aucune idée de la façon de commencer à le faire. J'ai trouvé un code dans mon livre qui lui ressemble mais c'est un exercice complètement différent et je ne sais pas comment les relier les uns aux autres. Comment puis-je commencer à simuler des arrivées et comment savoir quand elles sont terminées? Je sais comment les stocker et calculer a, b, c, d en fonction de cela. Mais je ne sais pas comment j'ai réellement besoin de simuler la simulation de monte carlo. Quelqu'un pourrait-il m'aider à démarrer? Je sais que ce n'est pas un endroit où vos questions obtiennent une réponse pour vous, mais seulement résolues à la place. Mais le problème est que je ne sais pas par où commencer.

Un service d'assistance informatique représente un système de mise en file d'attente avec cinq assistants prenant les appels des clients. Les appels se produisent selon un processus de Poisson avec le taux moyen d'un appel toutes les 45 secondes. Les durées de service des 1er, 2e, 3e, 4e et 5e assistants sont toutes des variables aléatoires exponentielles avec les paramètres λ1 = 0,1, λ2 = 0,2, λ3 = 0,3, λ4 = 0,4 et λ5 = 0,5 min − 1, respectivement (le jth help desk assistant a λk = k / 10 min − 1). Outre les clients assistés, jusqu'à dix autres clients peuvent être mis en attente. Lorsque cette capacité est atteinte, les nouveaux appelants reçoivent un signal occupé. Utilisez les méthodes de Monte Carlo pour estimer les caractéristiques de performance suivantes,

a) la fraction de clients qui reçoivent un signal occupé;

(b) le temps de réponse attendu;

c) le temps d'attente moyen;

d) la proportion de clients servis par chaque assistant du service d'assistance;

EDIT: ce que j'ai jusqu'à présent est (pas beaucoup):

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++;
    }
}
user3485470
la source
Il ne s'agit pas exactement de "comment faire MC", mais connaissez-vous ce paquet: r-bloggers.com/… ? Il semble parfaitement adapté au type de problèmes que vous décrivez (bien qu'en utilisant un modèle différent).
Tim
J'essaie actuellement de résoudre ce problème sans bibliothèques externes, mais si je ne peux pas le faire, j'utiliserai le vôtre à coup sûr :)
user3485470
Montrez ce que vous avez fait jusqu'à présent. Vous ne pouvez pas simplement venir ici et demander une solution de travail à domicile.
Aksakal

Réponses:

22

C'est l'un des types de simulation les plus instructifs et amusants à effectuer: vous créez des agents indépendants dans l'ordinateur, les laissez interagir, gardez une trace de ce qu'ils font et étudiez ce qui se passe. C'est une merveilleuse façon d'en apprendre davantage sur les systèmes complexes, en particulier (mais sans s'y limiter) ceux qui ne peuvent pas être compris avec une analyse purement mathématique.

La meilleure façon de construire de telles simulations est la conception descendante.

Au plus haut niveau, le code devrait ressembler à quelque chose

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

(Cet exemple et tous les exemples suivants sont du code exécutable R , pas seulement du pseudo-code.) La boucle est une simulation événementielle : get.next.event()trouve tout "événement" d'intérêt et en transmet une description process, ce qui en fait quelque chose (y compris la journalisation de tout informations à ce sujet). Il revient TRUEtant que les choses fonctionnent bien; lors de l'identification d'une erreur ou de la fin de la simulation, il revient FALSE, terminant la boucle.

Si nous imaginons une implémentation physique de cette file d'attente, comme des personnes attendant un permis de mariage à New York ou un permis de conduire ou un billet de train presque n'importe où, nous pensons à deux types d'agents: les clients et les "assistants" (ou serveurs) . Les clients s'annoncent en se présentant; les assistants annoncent leur disponibilité en allumant une lumière ou une pancarte ou en ouvrant une fenêtre. Ce sont les deux types d'événements à traiter.

L'environnement idéal pour une telle simulation est un véritable environnement orienté objet dans lequel les objets sont mutables : ils peuvent changer d'état pour répondre indépendamment aux choses qui les entourent. Rest absolument terrible pour cela (même Fortran serait mieux!). Cependant, nous pouvons toujours l'utiliser si nous y prenons garde. L'astuce consiste à conserver toutes les informations dans un ensemble commun de structures de données accessibles (et modifiées) par de nombreuses procédures distinctes et interactives. J'adopterai la convention d'utilisation des noms de variables EN TOUTES MAJUSCULES pour ces données.

Le niveau suivant de la conception descendante consiste à coder process. Il répond à un seul descripteur d'événement 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)
}

Il doit répondre à un événement nul lorsqu'il get.next.eventn'a aucun événement à signaler. Sinon, processimplémente les "règles métier" du système. Il s'écrit pratiquement à partir de la description de la question. Comment cela fonctionne devrait nécessiter peu de commentaires, sauf pour souligner que nous devrons éventuellement coder des sous put.on.hold- routines et release.hold(implémenter une file d'attente client) et serve(implémenter les interactions client-assistant).

Qu'est-ce qu'un "événement"? Il doit contenir des informations sur qui agit, sur le type d’action qu’il entreprend et sur le moment où il se produit. Mon code utilise donc une liste contenant ces trois types d'informations. Cependant, il get.next.eventsuffit d'inspecter les temps. Il est uniquement responsable du maintien d'une file d'attente d'événements dans lesquels

  1. Tout événement peut être placé dans la file d'attente lors de sa réception et

  2. Le premier événement de la file d'attente peut facilement être extrait et transmis à l'appelant.

La meilleure implémentation de cette file d'attente prioritaire serait un tas, mais c'est trop difficile R. Suite à une suggestion dans The Art of R Programming de Norman Matloff (qui propose un simulateur de file d'attente plus flexible, abstrait mais limité), j'ai utilisé un bloc de données pour contenir les événements et simplement le rechercher le temps minimum parmi ses enregistrements.

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)
}

Il y a plusieurs façons de coder cela. La version finale présentée ici reflète un choix que j'ai fait en codant comment processréagit à un événement "Assistant" et comment new.customerfonctionne: il get.next.eventsuffit de retirer un client de la file d'attente, puis de s'asseoir et d'attendre un autre événement. Il sera parfois nécessaire de rechercher un nouveau client de deux manières: premièrement, pour voir si l'on attend à la porte (pour ainsi dire) et deuxièmement, si on est entré alors que nous ne regardions pas.

De toute évidence, new.customeret next.customer.timesont des routines importantes , alors prenons-en soin ensuite.

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
}

CUSTOMERSest un tableau 2D, avec des données pour chaque client en colonnes. Il comporte quatre lignes (faisant office de champs) qui décrivent les clients et enregistrent leurs expériences au cours de la simulation : "Arrivé", "Servi", "Durée" et "Assistant" (un identificateur numérique positif de l'assistant, le cas échéant, qui a servi et sinon -1pour les signaux occupés). Dans une simulation très flexible, ces colonnes seraient générées dynamiquement, mais en raison de la façon dont Raime travailler, il est pratique de générer tous les clients au départ, dans une seule grande matrice, avec leurs heures d'arrivée déjà générées au hasard. next.customer.timepeut jeter un œil à la colonne suivante de cette matrice pour voir qui vient ensuite. La variable globaleCUSTOMER.COUNTindique le dernier client arrivé. Les clients sont gérés très simplement au moyen de ce pointeur, l'avançant pour obtenir un nouveau client et regardant au-delà (sans avancer) pour jeter un coup d'œil au client suivant.

serve implémente les règles métier dans la simulation.

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)
}

C'est simple. ASSISTANTSest une trame de données avec deux champs: capabilities(donnant leur taux de service) et available, qui marque la prochaine fois que l'assistant sera libre. Un client est servi en générant une durée de service aléatoire en fonction des capacités de l'assistant, en mettant à jour l'heure à laquelle l'assistant devient disponible et en enregistrant l'intervalle de service dans la CUSTOMERSstructure de données. Le VERBOSEdrapeau est pratique pour les tests et le débogage: lorsqu'il est vrai, il émet un flux de phrases en anglais décrivant les points de traitement clés.

La façon dont les assistants sont affectés aux clients est importante et intéressante. On peut imaginer plusieurs procédures: assignation au hasard, par une commande fixe, ou selon qui a été libre le plus longtemps (ou le plus court). Beaucoup d'entre eux sont illustrés dans du code commenté:

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)
}

Le reste de la simulation n'est en fait qu'un exercice de routine pour persuader Rde mettre en œuvre des structures de données standard, principalement un tampon circulaire pour la file d'attente en attente. Parce que vous ne voulez pas vous fâcher avec les globaux, j'ai mis tout cela en une seule procédure sim. Ses arguments décrivent le problème: le nombre de clients à simuler ( n.events), le taux d'arrivée des clients, les capacités des assistants et la taille de la file d'attente (qui peut être définie sur zéro pour éliminer complètement la file d'attente).

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

Il renvoie une liste des structures de données maintenues pendant la simulation; le plus intéressant est le CUSTOMERStableau. Rpermet de tracer assez facilement les informations essentielles dans ce tableau de manière intéressante. Voici une sortie montrant les derniers clients dans une simulation plus longue de clients.25050250

Figure 1

L'expérience de chaque client est tracée sous la forme d'une ligne de temps horizontale, avec un symbole circulaire au moment de l'arrivée, une ligne noire continue pour toute attente, et une ligne colorée pour la durée de leur interaction avec un assistant (la couleur et le type de ligne différencier les assistants). Sous ce graphique des clients se trouve celui montrant les expériences des assistants, marquant les moments où ils étaient et n'étaient pas engagés avec un client. Les points limites de chaque intervalle d'activité sont délimités par des barres verticales.

Lorsqu'il est exécuté avec verbose=TRUE, la sortie texte de la simulation ressemble à ceci:

...
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 
...

(Les nombres à gauche sont les heures d'émission de chaque message.) Vous pouvez faire correspondre ces descriptions aux parties du tracé Clients situées entre les heures et .165160165

Nous pouvons étudier l'expérience des clients en attente en traçant les durées d'attente par identifiant client, en utilisant un symbole spécial (rouge) pour montrer aux clients qui reçoivent un signal occupé.

Figure 2

(Tous ces graphiques ne feraient-ils pas un merveilleux tableau de bord en temps réel pour quiconque gère cette file d'attente de services!)

Il est fascinant de comparer les graphiques et les statistiques que vous obtenez en variant les paramètres transmis sim. Que se passe-t-il lorsque les clients arrivent trop rapidement pour être traités? Que se passe-t-il lorsque la file d'attente est réduite ou supprimée? Qu'est-ce qui change lorsque les assistants sont sélectionnés de différentes manières? Comment le nombre et les capacités des assistants influencent-ils l'expérience client? Quels sont les points critiques où certains clients commencent à se faire refuser ou à se mettre en attente pendant longtemps?


Normalement, pour des questions d'autoformation évidentes comme celle-ci, nous nous arrêterions ici et laisserions les détails restants comme un exercice. Cependant, je ne veux pas décevoir les lecteurs qui sont peut-être allés si loin et qui sont intéressés à essayer cela par eux-mêmes (et peut-être à le modifier et à le développer à d'autres fins), donc ci-dessous est le code de travail complet.

(Le traitement sur ce site gâchera l'indentation sur toutes les lignes contenant un symbole , mais l'indentation lisible doit être restaurée lorsque le code est collé dans un fichier texte.)$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)
whuber
la source
2
+1 Incroyable! Pourriez-vous répondre à toutes les questions avec ce niveau d'exhaustivité et d'attention aux détails? Des rêves, juste des rêves ...
Aleksandr Blekh
+1 Que puis-je dire? Aujourd'hui, j'ai appris tant de choses intéressantes! Voulez-vous ajouter un livre pour une lecture plus approfondie, s'il vous plaît?
mugen
1
@mugen J'ai mentionné le livre de Matloff dans le texte. Cela pourrait être approprié pour ceux Rqui débutent et qui souhaitent une autre perspective (mais assez similaire) sur les simulations de files d'attente. En écrivant ce petit simulateur, je me suis beaucoup réfléchi à ce que j'avais appris en étudiant le code dans (la première édition de) le texte d' Andrew Tanenbaum Operating Systems / Design and Implementation. J'ai également appris des structures de données pratiques, telles que des tas, grâce aux articles de Jon Bentley dans CACM et sa série de livres Programming Pearls . Tanenbaum et Bentley sont de grands auteurs que tout le monde devrait lire.
whuber
1
@ Mugen, il y a un manuel en ligne gratuit sur la théorie des files d'attente par Moshe ici . Le cours sur les processus stochastoc discrets du professeur Gallager couvre également ce sujet sur le MIT OCW . Les conférences vidéo sont vraiment bonnes.
Aksakal
@whuber, une excellente réponse. Bien que je ne pense pas que vous puissiez faire des enfants ces jours-ci pour lire Tanenbaum et Bentley :)
Aksakal