Simulation (Project-1)
Author: Simeon Paynter
Last Updated: July 13, 2021
Problem 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.
Use simulation to generate 10000 service times and estimate the mean service time for yourself.
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.
Add the curve of the probability density function (which is a mixture of two exponential pdfs) of the service time. For a mixture of two exponential distributions, refer to formula (15) of the paper.
What do we know?
From the problem description, we know the service times are exponentially distributed. This means when we run our simulation, we can use Algorithm 5.2 to generate continuous variables.
Algorithm 5.2 says to obtain a standard uniform random variable from a number generator. Then we compute \(F(X) = U for X\).
In order to estimate the mean service time, we’ll use Monte Carlo simulation to run the algorithm (5.2) over and over again: 10,000 times. The algorithm will use rexp(n, rate) to generate a random variable with exponential distribution, where \(rate = \lambda\). The rate will be in jobs per minute, either 1/12 or 1/3 (from \(3,12 = E(x) = 1 /\lambda\)) depending which mechanic is assigned.
First 50 values of the simulated data
The first 50 simulated service times have a mean of 4.9 (with a seed set to 100). After 50 simulations, we see the mean approaching the theoretical mean, as shown below.
## U Mechanic Jobs_per_minute simulated_service_time
## 1 0.30776611 2 - faster 0.33333333 2.1715116
## 2 0.55232243 2 - faster 0.33333333 9.2920870
## 3 0.81240262 2 - faster 0.33333333 3.5232880
## 4 0.54655860 2 - faster 0.33333333 5.2451723
## 5 0.62499648 2 - faster 0.33333333 0.5829794
## 6 0.76255108 2 - faster 0.33333333 1.0141303
## 7 0.20461216 2 - faster 0.33333333 3.3697398
## 8 0.35947511 2 - faster 0.33333333 1.1417432
## 9 0.53581115 2 - faster 0.33333333 1.2648231
## 10 0.53834870 2 - faster 0.33333333 1.4938334
## 11 0.42010145 2 - faster 0.33333333 5.2729682
## 12 0.77030161 2 - faster 0.33333333 0.5775103
## 13 0.48830599 2 - faster 0.33333333 0.7250846
## 14 0.69527414 2 - faster 0.33333333 0.3751463
## 15 0.98956414 2 - faster 0.33333333 4.2858160
## 16 0.33066053 2 - faster 0.33333333 1.6169414
## 17 0.60332436 2 - faster 0.33333333 2.5113704
## 18 0.30708590 2 - faster 0.33333333 3.0457997
## 19 0.19867907 1 - slower 0.08333333 18.9219753
## 20 0.25339065 2 - faster 0.33333333 6.6779364
## 21 0.46370118 2 - faster 0.33333333 0.8826072
## 22 0.96057309 2 - faster 0.33333333 1.0583890
## 23 0.44514802 2 - faster 0.33333333 3.3727269
## 24 0.45573146 2 - faster 0.33333333 2.5890973
## 25 0.41223704 2 - faster 0.33333333 3.0121520
## 26 0.57256477 2 - faster 0.33333333 1.2990224
## 27 0.77477889 2 - faster 0.33333333 2.0041626
## 28 0.09151028 1 - slower 0.08333333 13.3034199
## 29 0.98282408 2 - faster 0.33333333 8.9468138
## 30 0.57793740 2 - faster 0.33333333 1.3998850
## 31 0.24874240 2 - faster 0.33333333 2.6882798
## 32 0.73346670 2 - faster 0.33333333 0.4363017
## 33 0.44829914 2 - faster 0.33333333 0.8098162
## 34 0.12523909 1 - slower 0.08333333 44.3121929
## 35 0.38947869 2 - faster 0.33333333 9.2588306
## 36 0.36139663 2 - faster 0.33333333 0.4258685
## 37 0.68488024 2 - faster 0.33333333 0.0240082
## 38 0.83657146 2 - faster 0.33333333 1.8412005
## 39 0.08093264 1 - slower 0.08333333 16.9488819
## 40 0.91640643 2 - faster 0.33333333 1.3573595
## 41 0.20077621 2 - faster 0.33333333 2.0412045
## 42 0.39666235 2 - faster 0.33333333 3.7918687
## 43 0.47255689 2 - faster 0.33333333 0.5019415
## 44 0.35238717 2 - faster 0.33333333 12.3883945
## 45 0.55165118 2 - faster 0.33333333 8.1321146
## 46 0.23791515 2 - faster 0.33333333 1.0341318
## 47 0.57992109 2 - faster 0.33333333 10.5044177
## 48 0.47212988 2 - faster 0.33333333 9.3871801
## 49 0.46315442 2 - faster 0.33333333 0.7799686
## 50 0.67324929 2 - faster 0.33333333 7.4159888Theoretical mean service time
The theoretical mean service time is 3 min * (80%) + 12 min * (20%).
## [1] 4.8Summary graphs/tables
After all 10,000 service times have been simulated, we get a simulated mean within 0.2 of the theoretical mean. This isn’t because of luck, this is because Monte Carlo simulations help demonstrate long-term probabilities.
Simulated mean service time:
## [1] 4.781547Of the 10,000 services, there were not exactly 8,000 and 2,000 services performed by mechanic 2 (fast) and mechanic 1 (slow).
## 1 - slower 2 - faster
## 2016 7984These were probabilities, but in experiment we will not always get the exact probability. We can see the mean of U was close to 0.5.
## U Mechanic Jobs_per_minute
## Min. :0.0003098 Length:10000 Min. :0.08333
## 1st Qu.:0.2458366 Class :character 1st Qu.:0.33333
## Median :0.4972072 Mode :character Median :0.33333
## Mean :0.4979808 Mean :0.28293
## 3rd Qu.:0.7516270 3rd Qu.:0.33333
## Max. :0.9999495 Max. :0.33333
## simulated_service_time
## Min. : 0.00003
## 1st Qu.: 1.02360
## Median : 2.58136
## Mean : 4.78155
## 3rd Qu.: 5.50228
## Max. :89.19821Using a histogram, we can graphically see the distribution of the mean service times, which are exponentially distributed.
The blue line is the curve for the probability density function, which is a mixture of two pdfs from the link in part 3 in the project description.
Discussion
The results show that after 50 simulations, the mean of our simulated service times is close to the theoretical mean. After 10,000 simulations, the simulated mean service is very close to the theoretical mean.
Code
set.seed(100)
# initialize variables and data frames
n <- 10000 # simulates 10,000 times
df <- c() # empty df to store results
# monte carlo simulation: 10,000 service times
for( i in 1:n ){
# generates a uniform random variable between 0 and 1.
U <- runif(1)
# 20% of time, you get the slow mechanic whose rate is 12 minutes between each job
if(U < 0.2){
mech = "1 - slower"
rate = 1/12
simulated_service_time <- rexp(1, rate = rate) # rexp() produces random variable with exponential distribution
# 80% of the time, you get the fast mechanic whose rate is 3 minutes between each job
}else{
mech = "2 - faster"
rate = 1/3
simulated_service_time <- rexp(1, rate = rate) # rexp() produces random variable with exponential distribution
}
# Write results to a data frame
temp_df <- data.frame(U, Mechanic = mech, Jobs_per_minute = rate, simulated_service_time)
if (i == 1){
df <- temp_df
}
else {
df <- rbind(df, temp_df)
}
}
first_50 <- df[1:50, ]
# Average of 3 min, 80% of time + Average of 12 min, 20% of time
theoretical_mean <- 3 * 0.8 + 12 * 0.2
simulated_mean <- mean(df$simulated_service_time)
summary(as.factor(df$Mechanic))
summary(df)
hist(df$simulated_service_time, probability = TRUE)
abline(v = theoretical_mean, col="black", lwd=2)
abline(v = simulated_mean, col="red")
curve(expr = ((0.2)*(1/12)*exp(-(1/12)*x) + (0.8)*(1/3)*exp(-(1/3)*x)), add=TRUE, col="blue")
legend(40, 0.1, c("True mean", "Simulated mean"), text.col = c("black", "red"), bty = "n")