### 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 ### ### ###






### Problem 1: A frequency histogram of total travel times

# Create the plot and the aesthetic

  
  # Create a histogram with orange bins outlined in black, each of width 1
  
  
  # Define the labels as title, x-axis and y-axis




  
  # Apply a minimalist theme (good for academic and professional purposes)
 






### Problem 2: Create a line graph showing the distributions of the lengths of 
### delays due to congestion and due to incidents

# For a line graph we need x and y coordinates
# This means we have to calculate the frequencies of each possible delay value
# The easiest way to do this is to create a dataframe with each Dc or Dix value
# and a corresponding description
# The table() function is useful for creating frequency tables
# This uses geom_line and geom_point

# Add the 10000 Dix data to the end of the 10000 Dc data



# Put in 10000 labels for Dc, then another 10000 labels for Dix




# Build these into a dataframe




# Remove NA delays (Dix is NA when no incident occurred)




# Build frequency tables for each source




# Build a dataframe for the Dc frequencies





# Build a dataframe for the Dix frequencies





# Bind the two datasets together




# (Optional but helpful) Order by delay for each line




# Create the plot and the aesthetic




  # Create the lines




  # Create the points




  # Create the labels





  # Apply a minimalist theme







### Problem 3: Create a scatter plot showing the travel time verses incidence
### occurrence. This can be done using geom_point

# If you're doing something binary like this, it's best to have Incident as 
# a factor

# Convert incident to factor





# Set up dataset and aesthetic



  
  # Create points



  
  # Create labels




# Apply a minimalist theme




### Extra note: In cases like this where you expect many points to be on top
### of each other you can use geom_jitter to spread them out vertically

# Set up dataset and aesthetic


  # Create the points with added jitter
 

  
  # Create labels



  
  # Apply a minimalist theme