### 19-6 Exercise 2 / Homework
### Visualising the output of a commute time Monte Carlo simulator

# Include the ggplot2 library

library(ggplot2)

# set the pseudorandom seed for data reproducability

set.seed(137)

### Build the simulation

### The total timme taken for a commute can be calculated from the following:
###
### T = T0 + Dc + (I * Di)
###
### where
###
### T = total commute time
### T0 = minimum commute time with no problems (constant)
### Dc = Delay due to congestion (randomised variable)
### I = Probability that an incident occurs and causes a delay (randomised variable)
### Di = Delay due to incident (randomised variable)

### Don't worry about understanding the maths in this section - it's just to build
### data of the right shape

# Set number of trials n

n <- 10000

# Set minimum commute time T

T0 <- 30

# Set random delay due to congestion
# This uses a normal distribution with a mean delay of 8 minutes and a standard
# deviation of 4 minutes

Dc <- round(rnorm(n,8,4),0)

# Negative check - if delay is less than 0 make it 0

Dc <- ifelse(Dc<0,0,Dc)

# Set probability of delay due to an incident pI

pI <- 0.03   # 3% chance of incident

# For each trial, did an incident occur?

I <- runif(n,0,1)

# Logic check - if incidence is less than pI make it 0, else make it 1

I <- ifelse(I<=pI,1,0)

# Set random delay due to incident
# This uses a lognormal distribution with a mean delay of 10 minutes and a standard
# deviation of 5 minutes. NB - these values are converted to meanlog and sdlog.
# Lognormal functions are normal functions with a strong skewness to the right.

Di <- round(rlnorm(n, 
                   meanlog = log(10^2 / sqrt(5^2 + 10^2)), 
                   sdlog = sqrt(log(1 + (5^2 / 10^2)))),
            0)

# Negative check - if delay is less than 0 make it 0

Di <- ifelse(Di<0,0,Di)

### All elements of the simulation equation have now been calculated
### The final output value, T, can now be calculated

T <- T0 + Dc + (I * Di)

### ggplot2 expects all data to be in a dataframe, so bind everything into one

simdata <- data.frame(T = T, T0 = T0, Dc = Dc, I = I, Di = Di)

# It might be useful to have a column that only shows the actual delays due to
# incidents, not the potential delay, so lets add one called Dix

simdata$Dix <- ifelse(simdata$I>0,simdata$Di,NA)

### ### ### End of simulation model ### ### ###