12_extremalni_rozdeleni.qmd

# Rozdělení extrémních hodnot {#sec-rozdeleniextremnichhodnota}

## Vymezení extrémní hodnoty

### Bloková maxima

### Překročení limitu
Z dat jsou užita pouze pozorování překračující stanovenou mez, například 95 \% kvantil.

::: callout-tip
## Úloha

1.  a) Nahrajte do prostředí data z klementiské řady měření
    b) Metodou blokových maxim stanovte maxima 

:::


```{r, echo=FALSE, eval=TRUE}
klementinum <- readRDS("data/klementinum.rds")
klementinum$date <- as.POSIXct(paste(klementinum$yr, 
                             sprintf(fmt = "%02d", klementinum$mon),
                             sprintf(fmt = "%02d", klementinum$day),
                       sep = "-"), tz = "CET", format = "%Y-%m-%d")

threshold <- quantile(x = klementinum$tavg, 
                      probs = c(0.99, 0.05))
```


```{r, echo=FALSE, eval=FALSE}
plot(main = "Klementinum 1776 - 2023", 
     sub = "Identifikace extrémů s pomocí prahových hodnot",  
     x = klementinum$date, 
     y = klementinum$tavg, 
     type = "p", 
     cex = 0.1, 
     col = "#00000050", 
     ylab = expression(T[avg]), 
     xlab = "")
rect(xleft = min(klementinum$date), 
     xright = max(klementinum$date), 
     ytop = threshold[1], 
     ybottom = threshold[2], 
     col = "#ffffffee", 
     border = "#ffffff")
abline(h = threshold, 
       col = "#660033", 
       lwd = 2)
```

![](images/chapter_9_fig_1.png)


```{r}
klementinum_short <- klementinum |> 
  dplyr::filter(yr %in% 2000:2020) 
window <- 2000:2020


i <- seq(from = as.Date("2000-01-01"), to = as.Date("2020-12-31"), by = "1 year")
j <- ifelse(as.numeric(i) %% 2 == 0, "black", "#00000090")

plot(main = "Klementinum 1776 - 2023", 
     sub = "Identifikace extrémů s pomocí blokových maxim",  
     x = klementinum_short$date,
     y = klementinum_short$tavg, 
     type = "p", 
     cex = 0.1, 
     col = "#00000050", 
     ylab = expression(T[avg]), 
     xlab = "")
## draw rectangles with bottom left (100, 300)+i
## and top right (150, 380)+i
rect(0+i, 365+i, -10, 20, col = j, 
     border = j)
abline(h = threshold, 
       col = "#660033", 
       lwd = 2)
```


## Gumbellovo rozdělení

Hustota funkce Gumbellova rozdělení je dána předpisem

$$
f(x) = \dfrac{1}{d}\exp\left(-\dfrac{x-c}{d}\right)\cdot\exp\left[-exp\left(-\dfrac{x-c}{d}\right)\right]
$$
kde $c$ a $d$ jsou parametry.
Kvantilová funkce je pak
$$
x_T = c - d\cdot\ln\left[-\ln\left(1-\dfrac{1}{T}\right)\right]
$$

kde $T$ je doba opakování a parametry odhadnuté metodou mometů jsou:

$$
d = \dfrac{\sqrt{6}}{\pi}\sigma, \quad c = \mu-0.5772\cdot d
$$

a $\mu$ $\sigma$ jsou momenty celkového souboru.


```{r, fig.align='center'}
fitgumbel <- function(x, ...) {
  mX <- mean(x)
  sdX <- sd(x)
  dX <- sqrt(6)/pi*sdX
  cX <- mX - 0.5772*dX
  curve(cX - d*log(-log(1-1/x)), add = TRUE, ...)
}

plot(x = c(1,100), 
     y = c(0,100), 
     type = "n", 
     xlab = "Doba opakování (roky)", 
     ylab = "Maximální denní průtok v roce (mm/d)")
# fitgumbel()
```