Examples: Geometric distribution as applied to the 100-year flood. https://tonyladson.wordpress.com/2015/07/07/100-year-flood-geometric-distribution-number-of-years-to-the-first-flood/

OneHundredYearGeom.R

# Geometric Distribution of 1% floods


# probability of 4 or fewer safe years
pgeom(4, 0.01)
# [1] 0.04900995

# probability of exactly 4 safe years and then a flood in the 5th
dgeom(4, 0.01)
#[1] 0.00960596

# Probability of at least one flood in 5 years
# 1 minus the probability of zero floods in 5 years
1 - pbinom(0, 5, 0.01)

# Graphs

# Probability of being flooded as a function of the number of years living in the camp

y <- 0:20
plot(y+1, pgeom(y, 0.01), 
     pch = 21,
     bg = grey(0.8),
     type = 'p',
     las = 1,
     xlab = 'Number of years living in the flood prone camp',
     ylab = '',
     yaxt = 'n'
)
mtext(text = 'Probability of being flooded',side = 2, line = 4)
axis(side = 2, at = axTicks(2), label = paste(100*axTicks(2),'%', sep = ''), las = 1)
grid()

# Probability of being safe for a certain number of years and then being flooded in the next year
op <- par(oma = c(1,2,0,0))
y <- 0:20
plot(y, dgeom(y, 0.01), 
     pch = 21,
     bg = grey(0.8),
     type = 'p',
     las = 1,
     xlab = 'Number of years to wait before a 100-year flood',
     ylab = '',
     yaxt = 'n'
     )
mtext(text = 'Probability of being flooded in the next year',side = 2, line = 4)
axis(side = 2, at = axTicks(2), label = paste(100*axTicks(2),'%', sep = ''), las = 1)
grid()



# probability mass of the geometric distirubiton

p <- 0.01
my.mean = (1-p)/p
my.median = -1/(log(1 - p, 2)) - 1

op = par(oma = c(1, 2, 0, 0))
plot(dgeom(0:200, p),
     xlab = '',
     ylab = '',
     yaxt = 'n',
     pch = 19,
     cex = 0.5,
     bg = grey(0.5))
title(xlab = 'Years')
axis(side = 2, at = axTicks(2), label = paste0(100*axTicks(2), '%'), las = 1)
abline(v = my.mean, lty = 2) 
abline(v = my.median, lty = 3) 

my.median.label <- round(my.median, 1)
text(labels = bquote(paste('median = ',.(my.median.label))), x = 67 , y = 0.009, adj = 1)
text(labels = bquote(paste('mean = ',.(my.mean))), x = 100 , y = 0.008, adj = 0)
text(labels = 'mode = 0', x = 5, y = 0.01, adj = 0)
mtext(side = 2, line = 4, text = 'Probability')

par(op)

# simulation of the number of years till a flood occurs
library(ggplot2)

set.seed(2015)
df <- expand.grid(x = 1:100, y = 1:100)
df$z <- rbinom(100*100,1,0.01)


ggplot(data=df, aes(x,y)) + geom_tile(aes(fill=factor(z))) + 
  scale_x_continuous(breaks = NULL) +
  scale_y_continuous(breaks = NULL) +
  scale_fill_manual(values=c('light blue', 'dark blue')) +
  theme(legend.position="none") +
  xlab('') +
  ylab('') 
  
  

# Length of time between floods
z <- rle(df$z)
gap <- z$lengths[z$values == 0]
mean(gap)
# [1] 104.2632
median(gap)
# [1] 71
range(gap)
# [1]   1 415
hist(gap, breaks = 'FD')