project Description

Two mechanics are changing oil filters with different service times: - Fast mechanic: Exponential distribution with mean 3 minutes - Slow mechanic: Exponential distribution with mean 12 minutes
Customers have an 80% probability of being served by the faster mechanic.
We will simulate 10,000 service times and analyze the results.

N <- 10000 # Number of simulations

Step 1] Randomly assign to one of the mechanics (1 = fast, 0 = slow)

assignments <- rbinom(N, size = 1, prob = 0.8)

Step 2 Generate service times based on who was assigned 1 or 0

service_times <- ifelse( assignments == 1, rexp(N, rate = 1/3), # Faster mechanic (rate = 1/3) rexp(N, rate = 1/12) # Slower mechanic (rate = 1/12) )

Step 3 compute mean service time

print(mean(service_times))

##Theoretical Verification: “The expected service time can be calculated as: E[Service Time]= (0.8×3)+(0.2×12)=2.4 + 2.4= 4.8 minutes. Our simulation result (~4.76 minutes) is very close to the expected value.” ### Freq Histogram N <- 10000 assignments <- rbinom(N, 1, 0.8) x <- ifelse(assignments == 1, rexp(N, 1/3), rexp(N, 1/12))

Calculate means

simulated_mean <- mean(x) theoretical_mean <- 0.83 + 0.212 # 4.8 minutes #histogram (default) hist(x, main = “Frequency Histogram of Service Times”, xlab = “Service Time (minutes)”, ylab = “Frequency”, col = “plum”, border = “white”) abline(v = simulated_mean, col = “red”, lwd = 2, lty = 2) abline(v = theoretical_mean, col = “blue”, lwd = 2, lty = 2) legend(“topright”, legend = c(paste(“Simulated mean:”, round(simulated_mean, 2)), paste(“Theoretical mean:”, theoretical_mean)), col = c(“red”, “blue”), lwd = 2, lty = 2) # Density histogram (probability = TRUE) hist(x, probability = TRUE, main = “Density Histogram of Service Times”, xlab = “Service Time (minutes)”, ylab = “Density”, col = “lightgreen”, border = “white”) lines(density(x), col = “darkgreen”, lwd = 2) # Add density curve abline(v = simulated_mean, col = “red”, lwd = 2, lty = 2) abline(v = theoretical_mean, col = “blue”, lwd = 2, lty = 2) legend(“topright”, legend = c(“Density curve”, paste(“Simulated mean:”, round(simulated_mean, 2)), paste(“Theoretical mean:”, theoretical_mean)), col = c(“darkgreen”, “red”, “blue”), lwd = 2, lty = c(1, 2, 2)) # Create density histogram hist(service_times, probability = TRUE, main = “Service Time Distribution with Theoretical PDF”, xlab = “Service Time (minutes)”, ylab = “Density”, col = “lightgreen”, border = “white”, ylim = c(0, 0.12))

Add empirical density curve

lines(density(service_times), col = “deeppink”, lwd = 2)

Add theoretical PDF curve (mixture of two exponentials)

curve(0.8(1/3)exp(-(1/3)x) + 0.2(1/12)exp(-(1/12)x), col = “blue”, lwd = 2, add = TRUE)

Add vertical lines for means

abline(v = mean(service_times), col = “red”, lty = 2, lwd = 2) # Simulated mean abline(v = 4.8, col = “blue”, lty = 2, lwd = 2) # Theoretical mean

Add precisely positioned legend (using x,y coordinates)

legend(25, 0.10, # Positioned at x=25 minutes, y=0.10 density legend = c(“Empirical Density”, “Theoretical PDF”, paste(“Simulated mean:”, round(mean(service_times), 2)), “Theoretical mean: 4.8”), col = c(“red”, “blue”, “red”, “blue”), lty = c(1, 1, 2, 2), lwd = 2, bty = “n”) # Remove legend box