Distribution of Initial and continuing losses. Gist to support a blog https://tonyladson.wordpress.com/2019/07/23/the-distribution-of-losses/

Loss.R

# Generating initial loss values from the standardised initial losses from ARR
# Book 5
# Chapter 3.6.1
# Table 5.3.13

# The original source of this is Nathan et al. (2003)
# Nathan, R.J., Weinmann, P.E. and Hill, P.I. (2003), Use of a Monte-Carlo Simulation to
# estimate the Expected Probability of large to extreme floods, Proceedings of the 28th
# International Hydrology and Water Resources Symposium, pp: 1105-1112, 10-14 November,
# Wollongong.

# Also
# Sampling from an empirical distribution
# Chapter 7.2 of 
# Nathan, R. and Weinmann, E. (2013)
# Monte Carlo Simulation Techniques
# ARR Discussion Paper
#  http://arr.ga.gov.au/__data/assets/pdf_file/0017/40562/ARR_DiscussionPaper_MonteCarlo.pdf

# And
# Hill, P., Graszkiewicz, Z, Taylor, M., Nathan, R.  (2014)
# Project 6: Loss models for catchment simulation: Phase 4 analysis of rural catchments. 
# Australian Rainfall and Runoff Revision Projects.  Engineers Australia (link). 
# http://arr.ga.gov.au/__data/assets/pdf_file/0005/40496/ARR_Project6_Phase-4_Report.pdf

library(tidyverse)


loss_std <- structure(list(Percentile = c(100, 90, 80, 70, 60, 50, 40, 30, 
                              20, 10, 0), IL = c(0.14, 0.39, 0.53, 0.68, 0.85, 1, 1.2, 1.4, 
                                                 1.71, 2.26, 3.19), CL = c(0.15, 0.35, 0.48, 0.61, 0.79, 1, 1.24, 
                                                                           1.5, 1.88, 2.48, 3.85), prob = c(0, 0.1, 0.2, 0.3, 0.4, 0.5, 
                                                                                                            0.6, 0.7, 0.8, 0.9, 1)), class = c("tbl_df", "tbl", "data.frame"
                                                                                                            ), row.names = c(NA, -11L))
loss_std <- loss_std %>% 
  mutate(prob = (100-Percentile)/100) %>% 
  arrange(prob)

loss_std


# IL

p <- loss_std %>% 
  ggplot(aes(x = IL, y = prob)) +
  geom_point() +
  geom_line() +
  scale_x_continuous(name = 'Standardised Initial Loss', minor_breaks = seq(0,3.2,0.2)) +
  scale_y_continuous(name = 'Percentile') +
  geom_hline(yintercept = 1.00, linetype = 'dashed') +
  theme_grey(base_size = 7)

p

ggsave(here('figures', 'IL_percentile.png'), p, width = 4, height = 3)

# CL
p <- loss_std %>% 
  ggplot(aes(x = CL, y = prob)) +
  geom_point() +
  geom_line() +
  scale_x_continuous(name = 'Standardised Initial Loss', minor_breaks = seq(0,3.2,0.2)) +
  scale_y_continuous(name = 'Percentile') +
  geom_hline(yintercept = 1.00, linetype = 'dashed') +
  theme_grey(base_size = 7)

p



# Generate IL -------------------------------------------------------------



# 1. Generate a random number between 0 and 1
# 2. Linearly interpolate in the table

IL_f <- approxfun(x = 100-loss_std$Percentile, y= loss_std$IL)


# Random initial loss
IL_f(100*runif(10))

ILs <- tibble(IL = IL_f(100*runif(10000)))
mean(ILs$IL)
median(ILs$IL)

mean(ILs$IL>1.5)

# plot
hist(ILs$IL, breaks = 'FD')
# Check the FD breaks
hist(ILs$IL, breaks = 'FD', plot = FALSE)


# Create graph for blog
p <- ILs %>% 
  ggplot(aes(IL)) +
  geom_histogram(binwidth = 0.1, fill = 'blue', colour = 'black') +
  scale_x_continuous(name = 'Standardised Initial Loss', minor_breaks = seq(0,3.2,0.2)) +
  scale_y_continuous(name = 'Count') +
  theme_grey(base_size = 7)
  
p  
ggsave(here('figures', 'IL_hist.png'), p, width = 4, height = 3)



# Generate CL -------------------------------------------------------------

CL_f <- approxfun(x = 100-loss_std$Percentile, y= loss_std$CL)


# Random initial loss
CL_f(100*runif(10))

CLs <- tibble(CL = CL_f(100*runif(10000)))
mean(CLs$CL)
median(CLs$CL)

mean(CLs$CL>1.5)

# plot
hist(CLs$CL, breaks = 'FD')
# Check the FD breaks
hist(CLs$CL, breaks = 'FD', plot = FALSE)


# Create graph for blog
p <- CLs %>% 
  ggplot(aes(CL)) +
  geom_histogram(binwidth = 0.1, fill = 'blue', colour = 'black') +
  scale_x_continuous(name = 'Standardised Continuing loss Loss', minor_breaks = seq(0,4,0.2)) +
  scale_y_continuous(name = 'Count') +
  theme_grey(base_size = 7)

p
ggsave(here('figures', 'CL_hist.png'), p, width = 4, height = 3)