Symulacja Monte Carlo w R.

Próbuję rozwiązać następujące ćwiczenie, ale tak naprawdę nie mam pojęcia, jak zacząć to robić. W mojej książce znalazłem kod, który wygląda tak, ale jest to zupełnie inne ćwiczenie i nie wiem, jak je ze sobą powiązać. Jak rozpocząć symulowanie przylotów i skąd mam wiedzieć, kiedy są one ukończone? Wiem, jak je przechowywać i obliczyć odpowiednio a, b, c, d. Ale nie wiem, jak mam symulować symulację Monte Carlo. Czy ktoś mógłby mi pomóc zacząć? Wiem, że nie jest to miejsce, w którym otrzymujesz odpowiedzi na twoje pytania, ale tylko rozwiązane. Problem w tym, że nie wiem jak zacząć.

Dział pomocy technicznej IT reprezentuje system kolejkowania z pięcioma asystentami odbierającymi połączenia od klientów. Połączenia odbywają się zgodnie z procesem Poissona ze średnią szybkością jednego połączenia co 45 sekund. Czasy obsługi dla 1., 2., 3., 4. i 5. asystenta to wszystkie wykładnicze zmienne losowe o parametrach λ1 = 0,1, λ2 = 0,2, λ3 = 0,3, λ4 = 0,4 i λ5 = 0,5 min − 1, odpowiednio ( jth asystent pomocy technicznej ma λk = k / 10 min − 1). Oprócz klientów objętych pomocą można wstrzymać do dziesięciu innych klientów. W chwili osiągnięcia tej pojemności nowi dzwoniący otrzymują sygnał zajętości. Użyj metod Monte Carlo, aby oszacować następujące parametry wydajności,

(a) część klientów, którzy otrzymują sygnał zajętości;

(b) oczekiwany czas reakcji;

(c) średni czas oczekiwania;

(d) część klientów obsługiwanych przez każdego asystenta pomocy technicznej;

EDYCJA: to, co do tej pory mam (niewiele):

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++;
    }
}

Jest to jeden z najbardziej pouczających i zabawnych rodzajów symulacji do wykonania: tworzysz niezależnych agentów na komputerze, pozwalasz im na interakcję, śledzisz, co robią, i analizujesz, co się dzieje. Jest to cudowny sposób na poznanie złożonych systemów, zwłaszcza (ale nie tylko) tych, których nie można zrozumieć za pomocą analizy czysto matematycznej.

Najlepszym sposobem na zbudowanie takich symulacji jest projektowanie od góry w dół.

Na najwyższym poziomie kod powinien wyglądać mniej więcej tak

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

(Ten i wszystkie kolejne przykłady to kod wykonywalny R , a nie tylko pseudo-kod.) Pętla jest symulacją sterowaną zdarzeniami : get.next.event()znajduje dowolne „zdarzenie” będące przedmiotem zainteresowania i przekazuje jego opis process, który coś z nim robi (w tym rejestrowanie dowolnego informacje na ten temat). Powraca TRUEtak długo, jak wszystko działa dobrze; po zidentyfikowaniu błędu lub zakończeniu symulacji wraca FALSE, kończąc pętlę.

Jeśli wyobrażamy sobie fizyczną implementację tej kolejki, na przykład osoby czekające na prawo jazdy w Nowym Jorku lub prawo jazdy lub bilet kolejowy prawie wszędzie, myślimy o dwóch rodzajach agentów: klientów i „asystentów” (lub serwerów) . Klienci ogłaszają się, pokazując się; asystenci ogłaszają swoją dostępność, włączając światło, znak lub otwierając okno. Są to dwa rodzaje zdarzeń do przetworzenia.

Idealnym środowiskiem do takiej symulacji jest środowisko zorientowane obiektowo, w którym obiekty są zmienne : mogą zmieniać stan, aby reagować niezależnie na otaczające je rzeczy. Rjest do tego absolutnie okropny (nawet Fortran byłby lepszy!). Możemy jednak nadal z niego korzystać, jeśli zachowamy ostrożność. Sztuką jest utrzymanie wszystkich informacji we wspólnym zestawie struktur danych, do których można uzyskać dostęp (i je zmodyfikować) za pomocą wielu oddzielnych, oddziałujących na siebie procedur. Przyjmę konwencję używania nazw zmiennych WSZYSTKIMI KAPITAMI dla takich danych.

Kolejnym poziomem odgórnego projektu jest kodowanie process. Odpowiada na pojedynczy deskryptor zdarzenia 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)
}

Musi reagować na zdarzenie zerowe, gdy get.next.eventnie ma żadnych zdarzeń do zgłoszenia. W przeciwnym razie processimplementuje „reguły biznesowe” systemu. Praktycznie pisze się z opisu w pytaniu. To, jak to działa, powinno wymagać niewielkiego komentarza, z wyjątkiem zaznaczenia, że ​​w końcu będziemy musieli kodować podprogramy put.on.holdi release.hold(wdrażanie kolejki obsługującej klientów) i serve(wdrażanie interakcji klient-asystent).

Co to jest „wydarzenie”? Musi zawierać informacje o tym, kto działa, jaki rodzaj działania podejmuje i kiedy ma miejsce. Mój kod wykorzystuje zatem listę zawierającą te trzy rodzaje informacji. Jednakże, get.next.eventtylko musi sprawdzać czasy. Odpowiada tylko za utrzymanie kolejki zdarzeń, w której

1. Każde zdarzenie może zostać umieszczone w kolejce po jego odebraniu i

2. Najwcześniejsze zdarzenie w kolejce można łatwo wyodrębnić i przekazać dzwoniącemu.

Najlepszą implementacją tej kolejki priorytetowej byłaby kupa, ale to jest zbyt wybredne R. Zgodnie z sugestią w The Art of R Programming Normana Matloffa (który oferuje bardziej elastyczny, abstrakcyjny, ale ograniczony symulator kolejek), użyłem ramki danych do przechowywania zdarzeń i po prostu przeszukiwania go w celu znalezienia minimalnego czasu wśród jego rekordów.

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

Można to zakodować na wiele sposobów. Ostateczna wersja pokazana tutaj odzwierciedla wybór, którego dokonałem przy kodowaniu, jak processreaguje na zdarzenie „Asystent” i jak new.customerdziała: get.next.eventpo prostu usuwa klienta z kolejki wstrzymania, a następnie siada i czeka na kolejne zdarzenie. Czasami konieczne będzie poszukiwanie nowego klienta na dwa sposoby: po pierwsze, aby sprawdzić, czy ktoś czeka (przy drzwiach), a po drugie, czy ktoś wszedł, gdy nie patrzyliśmy.

Oczywiście, new.customeri next.customer.timeto ważne rutyny , więc niech się nimi zająć w przyszłym.

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
}

CUSTOMERSto tablica 2D z danymi dla każdego klienta w kolumnach. Ma cztery wiersze (działające jak pola), które opisują klientów i rejestrują ich doświadczenia podczas symulacji : „Przybył”, „Obsługiwany”, „Czas trwania” i „Asystent” (dodatni identyfikator numeryczny asystenta, jeśli taki był, który służył je, a poza tym -1dla sygnałów zajętości). W bardzo elastycznej symulacji kolumny te byłyby generowane dynamicznie, ale ze względu na to, jak Rlubi pracować, wygodnie jest wygenerować wszystkich klientów na początku, w jednej dużej matrycy, a czasy ich przybycia są już generowane losowo. next.customer.timemoże zajrzeć do następnej kolumny tej macierzy, aby zobaczyć, kto będzie następny. Zmienna globalnaCUSTOMER.COUNTwskazuje ostatniego klienta, który przybył. Klientami zarządza się bardzo prosto za pomocą tego wskaźnika, przesuwając go w celu pozyskania nowego klienta i spoglądając poza niego (bez przechodzenia), aby zerknąć na następnego klienta.

serve implementuje reguły biznesowe w symulacji.

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

To jest proste. ASSISTANTSto ramka danych z dwoma polami: capabilities(podająca stawkę usługi) i available, która oznacza następny czas, w którym asystent będzie wolny. Klient jest obsługiwany przez generowanie losowego czasu trwania usługi zgodnie z możliwościami asystenta, aktualizowanie czasu, kiedy asystent będzie następny dostępny, i rejestrowanie interwału usługi w CUSTOMERSstrukturze danych. VERBOSEFlaga jest przydatny do testowania i debugowania: gdy prawdziwe, to emituje strumień angielskich zdań opisujących kluczowych punktów przetwarzania.

Sposób przydzielania asystentów do klientów jest ważny i interesujący. Można sobie wyobrazić kilka procedur: losowe przydzielanie, pewne stałe porządkowanie lub według tego, kto był wolny najdłużej (lub najkrócej). Wiele z nich zostało zilustrowanych w skomentowanym kodzie:

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

Reszta symulacji to tak naprawdę zwykłe ćwiczenie polegające na przekonaniu Rdo wdrożenia standardowych struktur danych, głównie okrągłego bufora dla kolejki wstrzymania. Ponieważ nie chcesz uruchamiać amoku z globalsami, umieściłem je wszystkie w jednej procedurze sim. Argumenty opisują problem: liczbę klientów do symulacji ( n.events), wskaźnik przybycia klientów, możliwości asystentów i wielkość kolejki wstrzymania (którą można ustawić na zero, aby całkowicie wyeliminować kolejkę).

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

CUSTOMERSR$50$$250$

Doświadczenie każdego klienta jest wykreślane jako pozioma linia czasu, z okrągłym symbolem w momencie przybycia, ciągła czarna linia dla każdego oczekującego oczekiwania i kolorowa linia na czas interakcji z asystentem (kolor i rodzaj linii rozróżniać asystentów). Pod spiskiem Klientów znajduje się fabuła pokazująca doświadczenia asystentów, oznaczająca czasy, w których byli i nie byli zaangażowani z klientem. Punkty końcowe każdego przedziału aktywności są ograniczone pionowymi paskami.

Po uruchomieniu z verbose=TRUEtekstem symulacji wygląda następująco:

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

$160$$165$

Możemy badać doświadczenie klientów zawieszonych, wykreślając czasy oczekiwania według identyfikatora klienta, używając specjalnego (czerwonego) symbolu, aby pokazać klientom odbierającym sygnał zajętości.

(Czy wszystkie te wykresy nie byłyby wspaniałym pulpitem nawigacyjnym w czasie rzeczywistym dla każdego, kto zarządza tą kolejką usług!)

Porównywanie wykresów i statystyk, które uzyskujesz po zmianie parametrów, jest fascynujące sim. Co się stanie, gdy klienci pojawią się zbyt szybko, aby je przetworzyć? Co się stanie, gdy kolejka wstrzymania zostanie zmniejszona lub wyeliminowana? Jakie zmiany, gdy asystenci są wybierani na różne sposoby? Jak liczba i możliwości asystentów wpływają na doświadczenie klienta? Jakie są krytyczne punkty, w których niektórzy klienci zaczynają się odwracać lub zaczynają być zawieszani na długi czas?

Zwykle w przypadku oczywistych pytań do samodzielnego studiowania, takich jak to, zatrzymalibyśmy się tutaj i pozostawiliśmy pozostałe szczegóły jako ćwiczenie. Nie chcę jednak zawieść czytelników, którzy mogli zaatakować tak daleko i są zainteresowani wypróbowaniem tego na własną rękę (i być może modyfikacją i rozbudowaniem go do innych celów), dlatego poniżej znajduje się pełny działający kod.

$\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)