Simulating Bus Failures in R

[swaheera]

I am working with the R programming language.

• Suppose there are 100 buses (bus_1, bus_2,... bus_100).
• The bus company sends the first 5 buses (bus_1, bus_2... bus_5) on the first day.
• On the first day, each bus has a 0.5 probability of breaking down (once a bus breaks down, it is out of service permanently).
• If the bus does not break down on the first day, its sent out on the second day but now it has a probability of breaking down of 0.49. That is each day a bus doesn't break down, the probability of breaking down reduces by 0.01.
• When a bus breaks down, it is replaced with the next bus (e..g if bus1, bus2, bus3,bus4, bus5 are there bus2 and bus4 break on the same day, the next day we will have bus1, bus6, bus3, bus7 , bus5).
• At any given day, there can only be a max of 5 buses out on the road ... towards the end, there might be 4 buses, 3 buses ... until they all break down.

My Question: I am trying to write a simulation which simulates this situation until all buses break down? The final result should be a data frame with 11 columns: day_number (1,2,3...n), x1, x2, x3, x4, x5 (these represent the bus number in that position), p1,p2,p3,p4,p5 (these represent the probabilities of each bus breaking down on that day's row).

Here is what I tried so far using basic loops:

   bus_count <- 100
bus_prob <- rep(0.5, bus_count)
bus_active <- 1:5
results <- data.frame(matrix(ncol = 11, nrow = 0))
colnames(results) <- c("day_number", "x1", "x2", "x3", "x4", "x5", "p1", "p2", "p3", "p4", "p5")

day_number <- 0
while(length(bus_active) > 0) {
    day_number <- day_number + 1
    breakdown <- runif(length(bus_active)) < bus_prob[bus_active]
    bus_prob[bus_active] <- pmax(bus_prob[bus_active] - 0.01, 0)
    next_bus <- (max(bus_active)+1):(max(bus_active)+sum(breakdown))
    next_bus <- next_bus[next_bus <= bus_count]
    bus_prob[next_bus] <- 0.5
    bus_active <- c(bus_active[!breakdown], next_bus)
    results <- rbind(results, c(day_number, bus_active, bus_prob[bus_active]))
}

# rename bus values to actual names
results[,2:6] <- apply(results[,2:6], 2, function(x) paste0("bus_", x))

print(results)

This is not correct. I am noticing the following errors (column names did not come out correct):

Problem 1: The bus ordering is incorrect. On the first day, the buses should be bus1, bus2, bus3, bus4, bus5.

   X1      X3       X5       X6       X7       X8  X0.49 X0.49.1   X0.5 X0.5.1 X0.5.2
1   1   bus_3    bus_5    bus_6    bus_7    bus_8   0.49    0.49   0.50   0.50    0.5
2   2   bus_3    bus_6    bus_9   bus_10   bus_11   0.48    0.49   0.50   0.50    0.5

Problem 2: Buses are "coming back from the dead". E.g. Day 15 vs Day 47, bus47 comes back

15 15  bus_38   bus_45   bus_46   bus_47   bus_48   0.46    0.49   0.50   0.50    0.5
47 47  bus_47   bus_47   bus_47   bus_47   bus_47  47.00   47.00  47.00  47.00   47.0

Problem 3: The same bus appears multiple times in the same row:

47 47  bus_47   bus_47   bus_47   bus_47   bus_47  47.00   47.00  47.00  47.00   47.0

Problem 4: The probabilities are not in the correct ranges (e.g. can only be between 0.5 and 0)

37 37  bus_98   bus_99   bus_98  bus_0.5  bus_0.5   0.50   37.00  98.00  99.00   98.0
38 38  bus_99   bus_98  bus_100 bus_0.49 bus_0.49   0.50   38.00  99.00  98.00  100.0

Can someone please help me fix these problems and write this code correctly?

Thanks!

[vedoa]

Hi @swaheera

i would do it like this

• Since you decrement by 0.01 it can happen (very special occasion) that you have an never ending rotation since the probability drops to 0 (never breaks)
• all is done with base R
• you can change the number of busses per day and overall busses
• set.seed is to replicate the results - for me it stops on day 42 as everything is broken with bus 98 being the last standing one

Hope it gives some pointers for you.

set.seed(1)
bus_count <- 100
prob_start <- 0.5
prob_decrement <- 0.01
bus_per_shift <- 10
df <- data.frame(busId = 1:bus_count, 
                 daysWithoutBreak = 0, 
                 broken = FALSE, 
                 prob_current = prob_start)
day <- 1
shift <- data.frame(row.names = paste0("BusNumber_",1:bus_per_shift), Day1 = 1:bus_per_shift)

# Day 1
brokenDay1 <- rbinom(bus_per_shift, 1, prob_start)
idxBroken <- which(brokenDay1 == 1)
idxNotBroken <- which(brokenDay1 == 0)
if(!(length(idxBroken)==0)){
  df$broken[idxBroken] <- TRUE
}
if(!(length(idxNotBroken)==0)){
  df$daysWithoutBreak[idxNotBroken] <- df$daysWithoutBreak[idxNotBroken] + 1
  tempDay1 <- df$prob_current[idxNotBroken] - prob_decrement
  df$prob_current[idxNotBroken] <- ifelse(tempDay1 >= 0, tempDay1, 0)
}
# all other days
while(sum(!df$broken) > 0 & sum(df$prob_current) > 0) {
  day <- day + 1
  idx <- which(df$broken == FALSE)[1:bus_per_shift]
  shift[[paste0("Day",day)]] <- idx
  if(!(sum(df$prob_current[idx[!is.na(idx)]]) > 0)){
    cat("this shift goes on forever")
    break 
  }
  brokenDayN <- unlist(Map(rbinom, 1, 1, df$prob_current[idx]))
  idxBroken <- idx[which(brokenDayN == 1)]
  idxNotBroken <- idx[which(brokenDayN == 0)]
  if(!(length(idxBroken)==0)){
    df$broken[idxBroken] <- TRUE
  }
  if(!(length(idxNotBroken)==0)){
    df$daysWithoutBreak[idxNotBroken] <- df$daysWithoutBreak[idxNotBroken] + 1
    tempDayN <- df$prob_current[idxNotBroken] - prob_decrement
    df$prob_current[idxNotBroken] <- ifelse(tempDayN >= 0, tempDayN, 0)
  }    
}