Original Description

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.

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

  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. A histogram can be made based on frequencies. A density histogram can also be made. To get a density histogram, the height of each bar is divided by the total number of values, and then by the length of the bar. Such a histogram has a vertical scale comparable with the probability density function of the distribution from which data are randomly drawn.

  3. Add the curve of the probability density function (which is a mixture of two exponential pdfs) of the service time to the density histogram you got in question (2). For a mixture of two exponential distributions, refer to formula (15) of the paper http://www.columbia.edu/~ww2040/OnApproximationsIII84.pdf.

Code

#Part 1
N=10000 #number of service times generate
x=NULL #keeps track of service times
mech1=12 #mean time for 1st mechanic
mech2=3 #mean time for 2nd mechanic
for (k in 1:N){
  U=runif(1)
  if(U<=0.8){ 
    x[k]=rexp(1, rate=1/3) #faster mechanic (mech2)
  }
  else{ 
    x[k]=rexp(1, rate=1/12) #slower mechanic (mech1)
  }
}
mean(x) #simulated mean
## [1] 4.811821
#Part 2
tm=3*0.8+12*0.2 #theoreticalmean
hist(x, probability=TRUE, main="Histogram")
abline(v=mean(x), col=2, lty=2, lwd=2) #simulated mean vertical line
abline(v=tm, col=3, lty=3, lwd=2) #theoretical mean vertical line
legend(60, 0.05, c("Simulated Mean", "Theoretical Mean"), col=2:3, text.col=2:3, lty=2:3, lwd=2, bty="N")
#Part 3
f=function(x){
  (0.8*(1/12)*exp(-0.8*x)+0.2*(1/3)*exp(-0.2*x))
}
curve(f,add=TRUE, col="green")

The First 50 Data Values

x[1:50]
##  [1]  2.94415929  2.32334825  4.94899451 49.83670903  1.19360083  5.72420284
##  [7]  1.03408533  3.68973448  0.09474122  3.36302351  1.65074115  0.70274181
## [13]  0.49945835  2.95049683  3.25949707  0.22707443 17.93324150  2.26274207
## [19] 16.15341886  4.71768522 25.19880511  5.92117645  2.08367251 14.77324313
## [25]  2.78472058 19.08816614  2.48782360  4.07146369  0.91802677  6.61332032
## [31] 28.40459532  1.61297306  2.14543670  6.73162362  3.44906266  8.99383795
## [37]  6.79880357 14.27310791  5.18692500  0.61112773  2.93643784  1.30199096
## [43]  3.91259967 10.16910414  1.66875284 11.12873679 11.14177184  3.42569406
## [49] 14.07864705  0.73467506

Discussion

When the R Code is ran, the simulated mean varies due to the random generation. Although, the simulated mean continuously resulted to differ only very slightly from the theoretical mean, which is 4.8. Thus, the simulation seems to have worked effectively since the theoretical mean is quite close to the various simulated means after running this code.This is shown on the histogram with the simulated mean (red line) being nearly plotted on the theoretical mean (green line).