### 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
ggplot(simdata, aes(x = T)) +
  
  # Create a histogram with orange bins outlined in black, each of width 1
  geom_histogram(binwidth = 1, fill="orange", colour="black") +
  
  # Define the labels as title, x-axis and y-axis
  labs(title="Distribution of Commute Times",
       x = "Journey time",
       y = "Frequency") +
  
  # Apply a minimalist theme (good for academic and professional purposes)
  theme_minimal()






### 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
delay = c(simdata$Dc, simdata$Dix)

# Put in 10000 labels for Dc, then another 10000 labels for Dix
source = c(rep("Dc (congestion)", nrow(simdata)),
           rep("Dix (incident)",   nrow(simdata)))

# Build these into a dataframe
plotdata <- data.frame(delay = delay, source = source)

# Remove NA delays (Dix is NA when no incident occurred)
plotdata <- plotdata[!is.na(plotdata$delay), ]

# Build frequency tables for each source
tab_dc <- table(plotdata$delay[plotdata$source == "Dc (congestion)"])
tab_dix <- table(plotdata$delay[plotdata$source == "Dix (incident)"])

# Build a dataframe for the Dc frequencies
freq_dc <- data.frame(
  delay = as.integer(names(tab_dc)),
  freq = as.integer(tab_dc),
  source = "Dc (congestion)"
)

# Build a dataframe for the Dix frequencies
freq_dix <- data.frame(
  delay = as.integer(names(tab_dix)),
  freq = as.integer(tab_dix),
  source = "Dix (incident)"
)

# Bind the two datasets together
freqdata <- rbind(freq_dc, freq_dix)

# (Optional but helpful) Order by delay for each line
freqdata <- freqdata[order(freqdata$source, freqdata$delay), ]

# Create the plot and the aesthetic
ggplot(freqdata, aes(x = delay, y = freq, colour = source, group = source)) +

  # Create the lines
  geom_line(linewidth = 1) +

  # Create the points
  geom_point(size = 1.8) +

  # Create the labels
  labs(title  = "Frequency Distributions of Commute Delays",
       x      = "Delay (minutes)",
       y      = "Frequency",
       colour = "Delay source") +

  # Apply a minimalist theme
  theme_minimal()






### 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
simdata$Incident <- factor(
  simdata$I,
  levels = c(0, 1),
  labels = c("No incident", "Incident")
)

# Set up dataset and aesthetic
ggplot(simdata, aes(x = T, y = Incident)) +
  
  # Create points
  geom_point(alpha = 0.4) +
  
  # Create labels
  labs(
    title = "Total Travel Time vs Incident Occurrence",
    x = "Total travel time (minutes)",
    y = "Incident occurred"
  ) +
  
  
  theme_minimal()

### 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
ggplot(simdata, aes(x = T, y = Incident)) +

  # Create the points with added jitter
  geom_jitter(height = 0.1, alpha = 0.4) +
  
  # Create labels
  labs(
    title = "Total Travel Time vs Incident Occurrence",
    x = "Total travel time (minutes)",
    y = "Incident occurred"
  ) +
  
  # Apply a minimalist theme
  theme_minimal()