Simulação de Monte Carlo em R

Estou tentando resolver o exercício a seguir, mas na verdade não tenho idéia de como começar a fazer isso. Encontrei um código no meu livro que parece com isso, mas é um exercício completamente diferente e não sei como relacioná-lo um com o outro. Como posso começar a simular chegadas e como sei quando elas terminam? Eu sei como armazená-los e calcular a, b, c, d de acordo com isso. Mas não sei como preciso simular a simulação de monte carlo. Alguém poderia me ajudar a começar? Sei que este não é um lugar onde suas perguntas sejam respondidas, mas apenas resolvidas. Mas o problema é que não sei por onde começar.

Um help desk de suporte de TI representa um sistema de filas com cinco assistentes atendendo chamadas de clientes. As chamadas ocorrem de acordo com um processo de Poisson com a taxa média de uma chamada a cada 45 segundos. Os tempos de serviço para o 1º, 2º, 3º, 4º e 5º assistentes são todas variáveis ​​aleatórias exponenciais com os parâmetros λ1 = 0,1, λ2 = 0,2, λ3 = 0,3, λ4 = 0,4 e λ5 = 0,5 min-1, respectivamente (o O enésimo assistente de help desk possui λk = k / 10 min − 1). Além dos clientes que estão sendo assistidos, até dez outros clientes podem ser colocados em espera. Às vezes, quando essa capacidade é atingida, os novos chamadores recebem um sinal de ocupado. Use os métodos de Monte Carlo para estimar as seguintes características de desempenho,

(a) a fração de clientes que recebem um sinal de ocupado;

(b) o tempo de resposta esperado;

(c) o tempo médio de espera;

(d) a parcela de clientes atendidos por cada assistente de help desk;

EDIT: o que tenho até agora é (não muito):

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 é um dos tipos de simulação mais instrutivos e divertidos de executar: você cria agentes independentes no computador, permite que eles interajam, acompanhem o que fazem e estudam o que acontece. É uma maneira maravilhosa de aprender sobre sistemas complexos, especialmente (mas não limitado a) aqueles que não podem ser entendidos com análise puramente matemática.

A melhor maneira de construir essas simulações é com o design de cima para baixo.

No nível mais alto, o código deve parecer algo como

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

(Este e todos os exemplos subseqüentes são código executável R , não apenas pseudo-código.) O loop é uma simulação orientada a eventos : get.next.event()encontra qualquer "evento" de interesse e passa uma descrição para ele process, o que faz algo com ele (incluindo registrar qualquer informações sobre ele). Ele retorna TRUEenquanto as coisas estiverem funcionando bem; ao identificar um erro ou no final da simulação, ele retorna FALSE, encerrando o loop.

Se imaginarmos uma implementação física dessa fila, como pessoas aguardando uma licença de casamento na cidade de Nova York ou uma carteira de motorista ou passagem de trem em quase qualquer lugar, pensamos em dois tipos de agentes: clientes e "assistentes" (ou servidores) . Os clientes se anunciam aparecendo; os assistentes anunciam sua disponibilidade acendendo uma luz ou sinal ou abrindo uma janela. Esses são os dois tipos de eventos a serem processados.

O ambiente ideal para tal simulação é um verdadeiro ambiente orientado a objetos no qual os objetos são mutáveis : eles podem mudar de estado para responder independentemente às coisas ao seu redor. Ré absolutamente terrível para isso (até Fortran seria melhor!). No entanto, ainda podemos usá-lo se tomarmos alguns cuidados. O truque é manter todas as informações em um conjunto comum de estruturas de dados que podem ser acessadas (e modificadas) por muitos procedimentos separados e em interação. Adotarei a convenção de usar nomes de variáveis ​​IN ALL CAPS para esses dados.

O próximo nível do design de cima para baixo é codificar process. Ele responde a um único descritor 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)
}

Ele deve responder a um evento nulo quando get.next.eventnão houver eventos a serem relatados. Caso contrário, processimplementa as "regras de negócios" do sistema. Ele praticamente se escreve a partir da descrição na pergunta. Como funciona deve exigir pouco comentário, exceto para salientar que, eventualmente, precisaremos codificar sub-rotinas put.on.holde release.hold(implementar uma fila de espera do cliente) e serve(implementar as interações cliente-assistente).

O que é um "evento"? Ele deve conter informações sobre quem está agindo, que tipo de ação está sendo executada e quando está ocorrendo. Meu código, portanto, usa uma lista contendo esses três tipos de informações. No entanto, get.next.eventsó precisa inspecionar os horários. É responsável apenas por manter uma fila de eventos em que

1. Qualquer evento pode ser colocado na fila quando é recebido e

2. O evento mais antigo da fila pode ser facilmente extraído e passado para o chamador.

A melhor implementação dessa fila de prioridade seria uma pilha, mas isso é muito complicado R. Seguindo uma sugestão de The Art of R Programming, de Norman Matloff (que oferece um simulador de filas mais flexível, abstrato, mas limitado), usei um quadro de dados para armazenar os eventos e simplesmente pesquisá-lo pelo tempo mínimo entre seus 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)
}

Há muitas maneiras pelas quais isso poderia ter sido codificado. A versão final mostrada aqui reflete uma escolha que eu fiz na codificação de como processreage a um evento "Assistant" e como new.customerfunciona: get.next.eventapenas retira um cliente da fila de espera, depois senta e espera por outro evento. Às vezes, será necessário procurar um novo cliente de duas maneiras: primeiro, para ver se alguém está esperando na porta (por assim dizer) e, segundo, se alguém entrou quando não estávamos olhando.

Claramente, new.customere next.customer.timesão rotinas importantes , então vamos cuidar delas em seguida.

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
}

CUSTOMERSé uma matriz 2D, com dados para cada cliente em colunas. Possui quatro linhas (atuando como campos) que descrevem os clientes e registram suas experiências durante a simulação : "Chegou", "Servido", "Duração" e "Assistente" (um identificador numérico positivo do assistente, se houver, que serviu -los, e de outro modo -1para sinais de ocupado). Em uma simulação altamente flexível, essas colunas seriam geradas dinamicamente, mas, devido à maneira como as Rpessoas gostam de trabalhar, é conveniente gerar todos os clientes desde o início, em uma única matriz grande, com seus tempos de chegada já gerados aleatoriamente. next.customer.timepode dar uma olhada na próxima coluna desta matriz para ver quem vem em seguida. A variável globalCUSTOMER.COUNTindica o último cliente a chegar. Os clientes são gerenciados com muita simplicidade por meio desse ponteiro, avançando para obter um novo cliente e olhando além (sem avançar) para espiar o próximo cliente.

serve implementa as regras de negócios na simulação.

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

Isso é direto. ASSISTANTSé um quadro de dados com dois campos: capabilities(fornecendo a taxa de serviço) e available, que sinaliza na próxima vez em que o assistente estará livre. Um cliente é atendido gerando uma duração de serviço aleatória de acordo com os recursos do assistente, atualizando o horário em que o assistente se torna disponível e registrando o intervalo de serviço na CUSTOMERSestrutura de dados. O VERBOSEsinalizador é útil para teste e depuração: quando verdadeiro, emite um fluxo de frases em inglês que descreve os principais pontos de processamento.

Como os assistentes são atribuídos aos clientes é importante e interessante. Pode-se imaginar vários procedimentos: atribuição aleatória, por alguma ordem fixa ou de acordo com quem está livre há mais tempo (ou mais curto). Muitos destes são ilustrados no 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)
}

O restante da simulação é realmente apenas um exercício de rotina para persuadir Ra implementar estruturas de dados padrão, principalmente um buffer circular para a fila de espera. Como você não quer ficar louco com os globais, coloquei tudo isso em um único procedimento sim. Seus argumentos descrevem o problema: o número de clientes a serem simulados ( n.events), a taxa de chegada de clientes, as capacidades dos assistentes e o tamanho da fila de espera (que pode ser definida como zero para eliminar a fila completamente).

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

Retorna uma lista das estruturas de dados mantidas durante a simulação; o de maior interesse é a CUSTOMERSmatriz. Rtorna bastante fácil plotar as informações essenciais nessa matriz de uma maneira interessante. Aqui está uma saída que mostra os últimos clientes em uma simulação mais longa de clientes.$50$$250$

A experiência de cada cliente é plotada como uma linha do tempo horizontal, com um símbolo circular no momento da chegada, uma linha preta sólida para qualquer espera em espera e uma linha colorida pela duração da interação com um assistente (o tipo de cor e linha diferenciar entre os assistentes). Abaixo desse enredo de Clientes, há uma que mostra as experiências dos assistentes, marcando os horários em que eles estavam e não estavam envolvidos com um cliente. Os pontos finais de cada intervalo de atividade são delimitados por barras verticais.

Quando executada verbose=TRUE, a saída de texto da simulação se parece com:

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

(Os números à esquerda são os horários em que cada mensagem foi emitida.) Você pode corresponder essas descrições às partes da plotagem Customers entre os tempos e .$160$$165$

Podemos estudar a experiência dos clientes em espera, plotando as durações em espera pelo identificador do cliente, usando um símbolo especial (vermelho) para mostrar aos clientes que estão recebendo um sinal de ocupado.

(Todos esses gráficos não seriam um painel maravilhoso em tempo real para qualquer pessoa que gerencia essa fila de serviço!)

É fascinante comparar os gráficos e estatísticas que você obtém variando os parâmetros passados sim. O que acontece quando os clientes chegam muito rapidamente para serem processados? O que acontece quando a fila de espera é reduzida ou eliminada? O que muda quando os assistentes são selecionados de maneiras diferentes? Como os números e as capacidades dos assistentes influenciam a experiência do cliente? Quais são os pontos críticos em que alguns clientes começam a se afastar ou a ficar em espera por muito tempo?

Normalmente, para questões óbvias de auto-estudo como essa, paramos aqui e deixamos os detalhes restantes como um exercício. No entanto, não quero decepcionar os leitores que podem ter chegado tão longe e estão interessados ​​em tentar fazer isso sozinhos (e talvez modificá-lo e desenvolvê-lo para outros propósitos); portanto, em anexo abaixo está o código completo de funcionamento.

(O processamento neste site atrapalha a indentação em todas as linhas que contêm um símbolo , mas a indentação legível deve ser restaurada quando o código for colado em um arquivo 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)