set.seed(123)
# Number of simulations
n <- 10000
# Probabilities and means
p_faster <- 0.8
mean_faster <- 3
mean_slower <- 12
# Simulate the service times
service_times <- ifelse(runif(n) < p_faster, rexp(n, rate = 1/mean_faster), rexp(n, rate = 1/mean_slower))
# Display the first 50 simulated values
first_50_service_times <- service_times[1:50]
print(first_50_service_times)
## [1] 2.80654185 2.97368223 0.26138837 12.15831266 6.70129238 11.57743569
## [7] 0.78742121 6.08725224 1.30346734 5.61774300 6.92074467 4.04490518
## [13] 3.00851350 4.12936643 3.46665858 1.77933332 1.50879619 3.02913671
## [19] 4.21603684 2.84565162 0.72422371 1.12975997 4.92870661 4.88760015
## [25] 5.89203647 0.23580571 0.16918321 5.68575078 10.98087693 2.38747224
## [31] 41.21853179 8.63788196 1.75157292 0.29489719 2.25793441 1.82756745
## [37] 6.42691461 7.07578241 2.41906886 2.94403882 0.23447013 4.59718510
## [43] 0.04211809 6.55128600 2.25664708 1.50229495 3.19274140 0.81392848
## [49] 6.22694222 6.64737532
# Calculate the theoretical mean
theoretical_mean <- p_faster * mean_faster + (1 - p_faster) * mean_slower
# Calculate the simulated mean
simulated_mean <- mean(service_times)
# Plot histogram with density
hist(service_times, probability = TRUE, main = "Density Histogram of Service Times",
xlab = "Service Time", ylim = c(0, 0.2), breaks = 50)
# Add vertical lines for the means
abline(v = simulated_mean, col = "blue", lwd = 2, lty = 2)
abline(v = theoretical_mean, col = "red", lwd = 2, lty = 2)
# Add legend
legend("topright", legend = c("Simulated Mean", "Theoretical Mean"),
col = c("blue", "red"), lwd = 2, lty = 2)
# Add the mixture of exponential PDFs
curve(p_faster * dexp(x, rate = 1/mean_faster) + (1 - p_faster) * dexp(x, rate = 1/mean_slower),
add = TRUE, col = "darkgreen", lwd = 2)
# Print the means
cat("Simulated Mean:", simulated_mean, "\n")
## Simulated Mean: 4.79265
cat("Theoretical Mean:", theoretical_mean, "\n")
## Theoretical Mean: 4.8