An Original Description of the Problem

The original problem was two mechanics service oil filters with different mean times. Customers have a probability of being served by the faster mechanic. The goal is to simulate service times and compare them with theoretical expectations.

First 50 Values of the Simulated Data

2.75540704 13.77091563 0.78768420 4.00014043 1.05879673 4.76747013 4.69852369 3.54052731 2.67325557 0.65408601 0.57184704 18.35970060 0.06636194 2.58389111 1.37804827 0.67451215 2.34797157 3.62703941 3.67298820 30.07988439 11.29619527 4.66627120 1.01910955 6.61869901 2.97329201 3.45420786 0.20129083 0.99970615 1.52206145 1.63153887 16.14671453 4.36669256 5.31876880 0.27517701 3.59040610 1.25194722 10.79713133 2.98249845 8.79241671 5.20597381 1.73970501 3.06750545 3.43362047 1.78744529 2.29410437 0.59999070 7.26246459 5.38674280 0.53892809 10.28736350

Summary Graphs

Findings Discussion

The simulated mean closely approximates the theoretical mean derived from the mixture of exponential distributions. The histogram and PDF curve illustrate the distribution of service times, which aligns with theoretical expectations. The comparison between simulated and theoretical means provides insights into real-world service scenarios.

Code

num_simulations <- 10000 p_faster <- 0.8 lambda_faster <- 1/3 lambda_slower <- 1/12

service_times <- numeric(num_simulations) for (i in 1:num_simulations) { if (runif(1) < p_faster) { service_times[i] <- rexp(1, rate = lambda_faster) } else { service_times[i] <- rexp(1, rate = lambda_slower) } }

hist(service_times, probability = TRUE, main = “Simulation of Service Times”, xlab = “Service Time (minutes)”, ylab = “Density”) simulated_mean <- mean(service_times) theoretical_mean <- p_faster * (1 / lambda_faster) + (1 - p_faster) * (1 / lambda_slower) abline(v = simulated_mean, col = “red”, lty = 2, lwd = 2) abline(v = theoretical_mean, col = “blue”, lty = 2, lwd = 2)

pdf_mixture <- function(x) { p_faster * lambda_faster * exp(-lambda_faster * x) + (1 - p_faster) * lambda_slower * exp(-lambda_slower * x) }

x_values <- seq(0, max(service_times), length.out = 100)

lines(x_values, pdf_mixture(x_values), col = “green”, lty = 2, lwd = 2)

legend_x <- “topright” legend_y <- max(hist(service_times, plot = FALSE)$density) * 1.1 legend(legend_x, legend_y, legend = c(“Simulated Mean”, “Theoretical Mean”, “Mixture PDF”), col = c(“red”, “blue”, “green”), lty = 2, lwd = 2)

service_times[1:50] mean(service_times)