## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

## ----setup--------------------------------------------------------------------
library(distionary)
library(famish)

## ----data---------------------------------------------------------------------
x <- c(4.0, 2.7, 3.5, 3.2, 7.1, 3.1, 2.5, 5.0, 2.3, 4.5, 3.0, 3.8)

## ----fit-gev------------------------------------------------------------------
gev <- fit_dst_gev(x)
gev

## ----fit-lp3------------------------------------------------------------------
lp3 <- fit_dst("lp3", x, method = "lmom-log")
lp3

## ----density-plot, fig.width=6, fig.height=4----------------------------------
hist(x, freq = FALSE, ylim = c(0, 0.5), main = NULL, xlab = "Flow (cms)")
plot(gev, "density", add = TRUE, n = 400, lty = 2, col = "blue4")
plot(lp3, "density", add = TRUE, n = 400, lty = 3, col = "orange4")
legend(
  "topright",
  legend = c("Data histogram", "GEV density", "LP3 density"),
  lty = c(NA, 2, 3),
  pch = c(15, NA, NA),
  col = c("gray70", "blue4", "orange4"),
  pt.cex = c(1.5, NA, NA),
  bty = "n"
)

## ----return-levels------------------------------------------------------------
quantiles <- enframe_return(
  gev, lp3,
  at = c(2, 5, 10, 20, 50, 100, 200),
  arg_name = "return_period",
  fn_prefix = "flow"
)
quantiles

## ----freq-mag, fig.width=6, fig.height=4--------------------------------------
x_return_periods <- rpscore(x, pos = "Weibull")

plot(
  sort(x_return_periods), sort(x),
  xlab = "Return period (years)",
  ylab = "Flow (cms)",
  pch = 16, col = "black"
)
lines(quantiles$return_period, quantiles$flow_gev, col = "blue4", lty = 2, lwd = 2)
lines(quantiles$return_period, quantiles$flow_lp3, col = "orange4", lty = 3, lwd = 2)
legend(
  "topleft",
  legend = c("Empirical", "GEV", "LP3"),
  col = c("black", "blue4", "orange4"),
  pch = c(16, NA, NA),
  lty = c(NA, 2, 3),
  lwd = 2,
  bty = "n"
)

## ----quantile-score-----------------------------------------------------------
gev_100y <- quantiles$flow_gev[quantiles$return_period == 100]
lp3_100y <- quantiles$flow_lp3[quantiles$return_period == 100]

qs_gev <- quantile_score(x, gev_100y, tau = 1 - 1 / 100)
qs_lp3 <- quantile_score(x, lp3_100y, tau = 1 - 1 / 100)

c(GEV = mean(qs_gev), LP3 = mean(qs_lp3))