Introduction

Two mechanics are changing oil filters for the arriving customers. The service time has an Exponential distribution with mean 12 minutes for the first mechanic, and mean 3 minutes for the second mechanic. When you arrive to have your oil filter changed, your probability of being served by the faster mechanic is 0.8.

Part 1

Use simulation to generate 10000 service times and estimate the mean service time for yourself.

N=10000 #numbers of times it is simulated 
P1 = .8 #probability of the faster mechanic
P2 = .2 #probability of the slower mechanic
L1= 1/3 #lambda of the faster mechanic
L2= 1/12 #lambda of the slower mechanic

WhichMechanic = runif(N,0,1) #generates random numbers between 0 and 1.
Fast = WhichMechanic <.8 #creates a true or false for each customer.

time = numeric(N) #creates an empty array for the servise time. 

time[Fast] = rexp(sum(Fast), rate=L1) #Generate service times for customers assigned to the fast mechanic
time[!Fast] = rexp(sum(!Fast), rate=L2)#Generate service times for customers assigned to the slow mechanic

#displays the first 50 simulated service times. 
cat("First 50 simulated servise times:\n") 
## First 50 simulated servise times:
print(time[1:50])
##  [1]  0.85859220  0.90761921 12.55341779  1.44721072  0.22573712  1.30795865
##  [7]  0.62598079  0.63009230  4.45570934  6.01891606  0.06112106  0.36974863
## [13] 11.31099934  1.24067140  1.65961690  6.56261279  3.57342681  3.89351293
## [19]  2.27496175  1.08086301  0.19556428  1.89288624  1.73361727  1.91220699
## [25]  1.09145449  0.56349608  6.95654125  4.15937605  7.15359124  1.05610787
## [31]  5.41892007  0.69928533  6.30605540  2.90053747  0.75469683  7.32000309
## [37]  0.09578361  2.40654220  1.05195212  4.38616004 29.73185816  3.30616774
## [43]  4.26617135  0.75794470  0.41136227  5.95932823  1.38848575  6.32480077
## [49]  1.69446210  0.92593280
#simulated mean
simulatedMean = mean(time)
simulatedMean
## [1] 4.839786

Part 2

Summarize your data with R code hist(x, probability = TRUE) and mean(x), where x is an array representing your simulated data. Add to the histogram two vertical lines indicating the simulated mean and the theoretical mean.

N=10000 #numbers of times it is simulated 
P1 = .8 #probability of the faster mechanic
P2 = .2 #probability of the slower mechanic
L1= 1/3 #lambda of the faster mechanic
L2= 1/12 #lambda of the slower mechanic

WhichMechanic = runif(N,0,1) #generates random numbers between 0 and 1.
Fast = WhichMechanic <.8 #creates a true or false for each customer.

time = numeric(N) #creates an empty array for the service time. 

time[Fast] = rexp(sum(Fast), rate=L1) #Generate service times for customers assigned to the fast mechanic
time[!Fast] = rexp(sum(!Fast), rate=L2)#Generate service times for customers assigned to the slow mechanic

#displays the first 50 simulated service times. 
cat("First 50 simulated servise times:\n") 
## First 50 simulated servise times:
print(time[1:50])
##  [1] 10.99542830  1.32834616 10.20666958  4.58684065  3.37289266  0.02673063
##  [7]  5.69820067  0.69753669  2.20617374 23.47622421  2.91827035  1.81913930
## [13]  4.31444674  0.75058305  0.31938386  0.39875327  2.85578850  4.02180209
## [19]  8.78236330  1.07186677  5.53693737  1.77773942  0.15106583  1.45910650
## [25]  0.02998215 18.29433777  1.25695268  1.52591579  0.30237418  2.02429699
## [31]  0.80723802  0.17193155  5.94624625 17.98043091  0.77181718  5.19671150
## [37] 10.09532724  5.37026654  2.35908112  0.97618704  1.16722392  1.99027846
## [43]  2.52551315  0.24929546 13.88764568 42.52652610  5.18505309  1.18047357
## [49] 13.28918129  2.10044878
#density Histogram.  
hist(time, probability = TRUE,xlim =c(0,100)) #the probability = TRUE produces the density histogram. 

#simulated mean
simulatedMean = mean(time)
abline(v = simulatedMean, col = "blue", lwd = 2, lty = 2) #makes a vertical line at the simulated mean

#Theoretical mean
TheorMean = (P1*3)+(P2*12)
abline(v = TheorMean, col = "hotpink", lwd = 2) #makes a vertical line at the theoretical mean

legend('topright',
       legend = c("Simulated Mean", "Theoretical Mean"),
       col = c("blue","hotpink"), 
       lwd = 2, 
       lty = c(1,1,1), 
       bty = "n")

Part 3

Summarize your data with R code hist(x, probability = TRUE) and mean(x), where x is an array representing your simulated data. Add to the histogram two vertical lines indicating the simulated mean and the theoretical mean.

N=10000 #numbers of times it is simulated 
P1 = .8 #probability of the faster mechanic
P2 = .2 #probability of the slower mechanic
L1= 1/3 #lambda of the faster mechanic
L2= 1/12 #lambda of the slower mechanic

WhichMechanic = runif(N,0,1) #generates random numbers between 0 and 1.
Fast = WhichMechanic <.8 #creates a true or false for each customer.

time = numeric(N) #creates an empty array for the service time. 

time[Fast] = rexp(sum(Fast), rate=L1) #Generate service times for customers assigned to the fast mechanic
time[!Fast] = rexp(sum(!Fast), rate=L2)#Generate service times for customers assigned to the slow mechanic

#displays the first 50 simulated service times. 
cat("First 50 simulated servise times:\n") 
## First 50 simulated servise times:
print(time[1:50])
##  [1]  1.3459926  1.7877715  4.8344126  9.5325502  6.9318455  6.8382860
##  [7]  0.6887722  1.3506187  1.7293530 48.5932641  3.3534825  0.3768873
## [13]  1.3103898  0.9295194  9.7218978  1.6816068 10.3974814  7.8078668
## [19]  3.5375630  1.2111768  1.3736448  2.0881822  0.1230781  3.7039315
## [25]  0.4571957  6.5306451  1.1876801  8.7519079  2.4828681 11.1870258
## [31] 12.1628783  1.0324912  9.1934515  3.3819874  0.3788565  0.1141851
## [37]  2.6890452  1.5311458  2.7351164  1.2375901 14.2063466  0.8393695
## [43]  2.4400093  1.6661802  2.6147572  2.4644552  1.5463716  1.1081537
## [49]  2.3000062  1.3424675
#density Histogram.  
hist(time, probability = TRUE,xlim =c(0,100)) #the probability = TRUE produces the density histogram. 

#simulated mean
simulatedMean = mean(time)
abline(v = simulatedMean, col = "blue", lwd = 2, lty = 2) #makes a vertical line at the simulated mean

#Theoretical mean
TheorMean = (P1*3)+(P2*12)
abline(v = TheorMean, col = "hotpink", lwd = 2) #makes a vertical line at the theoretical mean

#Mixture of two exponential formula 
curve(P1*L1*exp(-L1*x)+P2*L2*exp(-L2*x), from =0, to = 100, col = "cyan", lwd = 2, add = TRUE)

legend('topright',
       legend = c("Simulated Mean", "Theoretical Mean", "Mixture"),
       col = c("blue","hotpink", "cyan"), 
       lwd = 2, 
       lty = c(1,1,1), 
       bty = "n")

Discussion

The simulated and theoretical results for the mean are near identical. This supports the accuracy of the theoretical formula for the mean of a mixture of exponential distributions.