sim_users <- function(users, periods, prob) {
  active <- rbinom(periods, size = users, prob = prob) 
  avg_active <- mean(active)
  return(avg_active)
}

users <- 500  # Number of users
periods <- 10  # Time periods
prob <- 0.2  # Probability of a user being active

avg_active_users <- sim_users(users, periods, prob)
print(paste0("Avg active users: ", format(avg_active_users, digits = 2)))

[1] "Avg active users: 99"

sim_users <- function(users, prob, periods) {
  active <- numeric(periods)
  
  for (i in 1:periods) {
    activity <- rbinom(users, 1, prob)
    active[i] <- sum(activity)
  }
  
  return(active)
}

sim_resource <- function(users, prob, periods) {
  blocked <- numeric(periods)
  
  for (i in 1:periods) {
    numbers <- sample(1:10, users, replace = TRUE)
    active_users <- sample(c(TRUE, FALSE), users, replace = TRUE, prob = c(prob, 1 - prob))
    
    indices <- which(active_users)
    active_numbers <- numbers[indices]
    
    if (length(indices) > 2) {
      lowest_two <- sort(active_numbers)[1:2]
      blocked_indices <- indices[!(active_numbers %in% lowest_two)]
      blocked[i] <- length(blocked_indices)
    } else {
      blocked[i] <- 0 # No users blocked if less than 2 active users
    }
  }
  
  return(blocked)
}

N_vals <- c(4, 8)
p_vals <- c(0.25, 0.5, 0.75)
periods <- 1000

results <- matrix(NA, nrow = length(N_vals), ncol = length(p_vals))
rownames(results) <- N_vals
colnames(results) <- p_vals

for (i in seq_along(N_vals)) {
  for (j in seq_along(p_vals)) {
    blocked_users <- sim_resource(N_vals[i], p_vals[j], periods)
    avg_blocked <- mean(blocked_users)
    results[i, j] <- avg_blocked
  }
}

print(results)

   0.25   0.5  0.75
4 0.056 0.327 0.932
8 0.416 1.832 3.730