Skip to content

Instantly share code, notes, and snippets.

@TonyLadson
Last active June 2, 2020 11:10
Show Gist options
  • Star 1 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save TonyLadson/9c9b908e395ca5bdd5fa9e6a4921d099 to your computer and use it in GitHub Desktop.
Save TonyLadson/9c9b908e395ca5bdd5fa9e6a4921d099 to your computer and use it in GitHub Desktop.
library(tidyverse)
#pbinom(q, size, prob, lower.tail = TRUE, log.p = FALSE)
# Probability of 1, 1% flood in 100 years
dbinom(1, size = 100, prob = 0.01, log = FALSE)
choose(100,1) *0.01^1*0.99^99
# Probability of 5 or more
1-pbinom(4, size = 100, 0.01)
# Create data frame for plotting
sampling_dist <- expand_grid(num_events = 0:10, years = 100, event_prob = 0.01)
sampling_dist <- sampling_dist %>%
rowwise() %>%
mutate(prob_num = dbinom(x = num_events, size = years, prob = event_prob))
p <- sampling_dist %>%
ggplot(aes(x = num_events, y = prob_num)) +
geom_col(fill = 'blue', colour = 'black') +
scale_x_continuous(name = 'Number of events', breaks = 0:10) +
scale_y_continuous(name = 'Probability of the number of events', limits = c(0, 0.4)) +
theme_gray(base_size = 7)
ggsave(file.path('../figures','Sampling_dist_1pc_100y.png'), p, width = 4, height = 3)
# Expected number of floods
# Sum x * Pr(x)
sampling_dist %>%
mutate(mean_component = num_events * prob_num) %>%
ungroup() %>%
dplyr::summarize(sum(mean_component), n = n())
# or
100*0.01
##################
# AEP = 2%
sampling_dist <- expand_grid(num_events = 0:10, years = 100, event_prob = 0.02)
sampling_dist <- sampling_dist %>%
rowwise() %>%
mutate(prob_num = dbinom(x = num_events, size = years, prob = event_prob))
sampling_dist %>%
ggplot(aes(x = num_events, y = prob_num)) +
geom_col(fill = 'blue', colour = 'black') +
scale_x_continuous(name = 'Number of events', breaks = 0:10) +
scale_y_continuous(name = 'Probability of the number of events', limits = c(0, 0.4))
p <- sampling_dist %>%
ggplot(aes(x = num_events, y = prob_num)) +
geom_col(fill = 'blue', colour = 'black') +
scale_x_continuous(name = 'Number of events', breaks = 0:10) +
scale_y_continuous(name = 'Probability of the number of events', limits = c(0, 0.4)) +
theme_gray(base_size = 7)
ggsave(file.path('../figures','Sampling_dist_2pc_100y.png'), p, width = 4, height = 3)
# -------------------------------------------------------------------------
#Change in Risk bewteen a 1% and 2% flood
# Probability of 1 or more floods
# Comparing 2% AEP to 1% AEP
1-pbinom(0, size = 100, prob = 0.02, log = FALSE) #0.867
1-pbinom(0, size = 100, prob = 0.01, log = FALSE) #0.633
(1-pbinom(0, size = 100, prob = 0.02, log = FALSE))/(1-pbinom(0, size = 100, prob = 0.01, log = FALSE))
# Probability of 5 or more floods
# Comparing 2% AEP to 1% AEP
(1-pbinom(4, size = 100, prob = 0.02, log = FALSE))
(1-pbinom(4, size = 100, prob = 0.02, log = FALSE))/(1-pbinom(4, size = 100, prob = 0.01, log = FALSE))
#14.8
# Compare 1% with 0.05%
1-pbinom(0, size = 100, prob = 0.005, log = FALSE) #0.39
1-pbinom(0, size = 100, prob = 0.01, log = FALSE) #0.633
(1-pbinom(4, size = 100, prob = 0.005, log = FALSE))
(1-pbinom(4, size = 100, prob = 0.005, log = FALSE))/(1-pbinom(4, size = 100, prob = 0.01, log = FALSE))
1/((1-pbinom(4, size = 100, prob = 0.005, log = FALSE))/(1-pbinom(4, size = 100, prob = 0.01, log = FALSE)))
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment