Commit 313e7813 authored by Felix Fauer's avatar Felix Fauer
Browse files

updated readme.rmd and manual files to new eta2_handling

parent c0d02a9c
......@@ -30,7 +30,7 @@
#' * The dependence of scale and location parameter on duration, \eqn{\sigma(d)} and \eqn{\mu(d)}, can be extended by multiscaling
#' and flattening, if requested. Multiscaling introduces a second duration exponent \eqn{\eta_2}, enabling the model to change slope
#' linearly with return period. Flattening adds a parameter \eqn{\tau}, that flattens the IDF curve for long durations:
#' \deqn{\sigma(x)=\sigma_0(d+\theta)^{-\eta_2}+\tau }
#' \deqn{\sigma(x)=\sigma_0(d+\theta)^{-(\eta+\eta_2)}+\tau }
#' \deqn{\mu(x)=\tilde{\mu}(\sigma_0(d+\theta)^{-\eta_1}+\tau)}
#' * A useful introduction to __Maximum Likelihood Estimation__ for fitting for the
#' generalized extreme value distribution (GEV) is provided by Coles (2001). It should be noted, however, that this method uses
......
......@@ -71,7 +71,7 @@ plot(ann.max$ds,ann.max$xdat,log='xy',xlab = 'Duration [h]',ylab='Intensity [mm/
```{r fit,warning=FALSE}
fit <- gev.d.fit(xdat = ann.max$xdat,ds = ann.max$ds,sigma0link = make.link('log'))
# checking the fit
gev.d.diag(fit,pch=1,)
gev.d.diag(fit,pch=1,cl=TRUE)
# parameter estimates
params <- gev.d.params(fit)
print(params)
......@@ -137,7 +137,7 @@ set.seed(42)
# durations
ds <- 1/60*2^(seq(0,13,1))
# random data for each duration
xdat <- sapply(ds,rgev.d,n = 20,mut = 2,sigma0 =3,xi = 0.2,theta = 0.1,eta = 0.6,tau = 0.1,eta2 = 0.8)
xdat <- sapply(ds,rgev.d,n = 20,mut = 2,sigma0 =3,xi = 0.2,theta = 0.1,eta = 0.6,tau = 0.1,eta2 = 0.2)
# transform to data.frame
example <- data.frame(xdat=as.numeric(xdat),ds=rep(ds,each=dim(xdat)[1]))
......@@ -171,6 +171,7 @@ for(i.fit in 1:length(all.fits)){
}
for(i.p in 1:length(idf.probs)){
# plotting IDF curves for each model (different colors) and probability (different lty)
print(gev.d.params(all.fits[[i.fit]]))
IDF.plot(1/60*2^(seq(0,13,0.5)),gev.d.params(all.fits[[i.fit]]),probs = idf.probs[i.p]
,add = TRUE,legend = FALSE,lty = i.p,cols = fit.cols[i.fit])
}
......
<!-- README.md is generated from README.Rmd. Please edit that file -->
# IDF
IDF
===
<!-- badges: start -->
<!-- badges: end -->
Intensity-duration-frequency (IDF) curves are a widely used analysis-tool in hydrology to assess the characteristics of extreme precipitation. The package 'IDF' functions to estimate IDF relations for given precipitation time series on the basis of a duration-dependent generalized extreme value (GEV) distribution. The central function is , which uses the method of maximum-likelihood estimation for the d-GEV parameters, whereby it is possible to include generalized linear modeling for each parameter. For more detailed information on the methods and the application of the package for estimating IDF curves with spatial covariates, see Ulrich et. al (2020, <https://doi.org/10.3390/w12113119>).
Intensity-duration-frequency (IDF) curves are a widely used
analysis-tool in hydrology to assess the characteristics of extreme
precipitation. The package ‘IDF’ functions to estimate IDF relations for
given precipitation time series on the basis of a duration-dependent
generalized extreme value (GEV) distribution. The central function is ,
which uses the method of maximum-likelihood estimation for the d-GEV
parameters, whereby it is possible to include generalized linear
modeling for each parameter. For more detailed information on the
methods and the application of the package for estimating IDF curves
with spatial covariates, see Ulrich et. al (2020,
<https://doi.org/10.3390/w12113119>).
## Installation
Installation
------------
You can install the released version of IDF from
[CRAN](https://CRAN.R-project.org) with:
You can install the released version of IDF from [CRAN](https://CRAN.R-project.org) with:
``` r
install.packages("IDF")
......@@ -34,14 +22,12 @@ or from gitlab using:
devtools::install_git("https://gitlab.met.fu-berlin.de/Rpackages/idf_package")
```
## Example
Example
-------
Here are a few examples to illustrate the order in which the functions
are intended to be used.
Here are a few examples to illustrate the order in which the functions are intended to be used.
- Step 0: sample 20 years of example hourly ‘precipitation’ data
<!-- end list -->
- Step 0: sample 20 years of example hourly 'precipitation' data
``` r
set.seed(999)
......@@ -50,9 +36,7 @@ sample.precip <- rgamma(n = length(dates), shape = 0.05, rate = 0.4)
precip.df <- data.frame(date=dates,RR=sample.precip)
```
- Step 1: get annual maxima
<!-- end list -->
- Step 1: get annual maxima
``` r
library(IDF)
......@@ -65,9 +49,7 @@ plot(ann.max$ds,ann.max$xdat,log='xy',xlab = 'Duration [h]',ylab='Intensity [mm/
<img src="man/figures/README-annualmax-1.png" width="100%" />
- Step 2: fit d-GEV to annual maxima
<!-- end list -->
- Step 2: fit d-GEV to annual maxima
``` r
fit <- gev.d.fit(xdat = ann.max$xdat,ds = ann.max$ds,sigma0link = make.link('log'))
......@@ -83,7 +65,7 @@ fit <- gev.d.fit(xdat = ann.max$xdat,ds = ann.max$ds,sigma0link = make.link('log
#> $se
#> [1] 4.018974e-01 6.677404e-02 4.931317e-02 2.000063e-06 1.141522e-02
# checking the fit
gev.d.diag(fit,pch=1,)
gev.d.diag(fit,pch=1,cl=TRUE)
```
<img src="man/figures/README-fit-1.png" width="100%" />
......@@ -92,8 +74,8 @@ gev.d.diag(fit,pch=1,)
# parameter estimates
params <- gev.d.params(fit)
print(params)
#> mut sigma0 xi theta eta eta2 tau
#> 1 6.478887 1.4648 -0.01833254 2.843666e-09 0.792231 0.792231 0
#> mut sigma0 xi theta eta eta2 tau
#> 1 6.478887 1.4648 -0.01833254 2.843666e-09 0.792231 0 0
# plotting the probability density for a single duration
q.min <- floor(min(ann.max$xdat[ann.max$ds%in%1:2]))
......@@ -112,9 +94,7 @@ legend('topright',col=1:2,lwd=1,legend = paste('d=',1:2,'h'),title = 'Duration')
<img src="man/figures/README-fit-2.png" width="100%" />
- Step 3: adding the IDF-curves to the data
<!-- end list -->
- Step 3: adding the IDF-curves to the data
``` r
plot(ann.max$ds,ann.max$xdat,log='xy',xlab = 'Duration [h]',ylab='Intensity [mm/h]')
......@@ -123,63 +103,41 @@ IDF.plot(durations,params,add=TRUE)
<img src="man/figures/README-idf-1.png" width="100%" />
## IDF Features
This Example depicts the different features that can be used to model
the IDF curves (submitted for publication: Fauer et al. An Extended
Model in Estimating Consistent Quantiles for
Intensity-Duration-Frequency Curves, 2020). Here we assume, that the
block maxima of each duration can be modeled with the GEV distribution
(![\\xi\\neq0](https://latex.codecogs.com/png.latex?%5Cxi%5Cneq0
"\\xi\\neq0")):
![G(z;\\mu,\\sigma,\\xi)=\\exp \\left\\lbrace -\\left\[ &#10;1+\\xi
\\left( \\frac{z-\\mu}{\\sigma} \\right)&#10;\\right\]^{-1/\\xi}
\\right\\rbrace,](https://latex.codecogs.com/png.latex?G%28z%3B%5Cmu%2C%5Csigma%2C%5Cxi%29%3D%5Cexp%20%5Cleft%5Clbrace%20-%5Cleft%5B%20%0A1%2B%5Cxi%20%5Cleft%28%20%5Cfrac%7Bz-%5Cmu%7D%7B%5Csigma%7D%20%5Cright%29%0A%5Cright%5D%5E%7B-1%2F%5Cxi%7D%20%5Cright%5Crbrace%2C
"G(z;\\mu,\\sigma,\\xi)=\\exp \\left\\lbrace -\\left[
IDF Features
------------
This Example depicts the different features that can be used to model the IDF curves (submitted for publication: Fauer et al. An Extended Model in Estimating Consistent Quantiles for Intensity-Duration-Frequency Curves, 2020). Here we assume, that the block maxima of each duration can be modeled with the GEV distribution (![\\xi\\neq0](https://latex.codecogs.com/png.latex?%5Cxi%5Cneq0 "\xi\neq0")):
![G(z;\\mu,\\sigma,\\xi)=\\exp \\left\\lbrace -\\left\[
1+\\xi \\left( \\frac{z-\\mu}{\\sigma} \\right)
\\right]^{-1/\\xi} \\right\\rbrace,")
\\right\]^{-1/\\xi} \\right\\rbrace,](https://latex.codecogs.com/png.latex?G%28z%3B%5Cmu%2C%5Csigma%2C%5Cxi%29%3D%5Cexp%20%5Cleft%5Clbrace%20-%5Cleft%5B%20%0A1%2B%5Cxi%20%5Cleft%28%20%5Cfrac%7Bz-%5Cmu%7D%7B%5Csigma%7D%20%5Cright%29%0A%5Cright%5D%5E%7B-1%2F%5Cxi%7D%20%5Cright%5Crbrace%2C "G(z;\mu,\sigma,\xi)=\exp \left\lbrace -\left[
1+\xi \left( \frac{z-\mu}{\sigma} \right)
\right]^{-1/\xi} \right\rbrace,")
where the GEV parameters depend on duration according to:
![&#10;\\sigma(d)=\\sigma\_0(d+\\theta)^{-\\eta\_2}+\\tau,
\\\\&#10;\\mu(d) =
\\tilde{\\mu}\\cdot\\sigma\_0(d+\\theta)^{-\\eta}+\\tau,
\\\\&#10;\\xi(d) = \\text{const.}
&#10;](https://latex.codecogs.com/png.latex?%0A%5Csigma%28d%29%3D%5Csigma_0%28d%2B%5Ctheta%29%5E%7B-%5Ceta_2%7D%2B%5Ctau%2C%20%5C%5C%0A%5Cmu%28d%29%20%3D%20%5Ctilde%7B%5Cmu%7D%5Ccdot%5Csigma_0%28d%2B%5Ctheta%29%5E%7B-%5Ceta%7D%2B%5Ctau%2C%20%20%5C%5C%0A%5Cxi%28d%29%20%3D%20%5Ctext%7Bconst.%7D%20%0A
"
\\sigma(d)=\\sigma_0(d+\\theta)^{-\\eta_2}+\\tau, \\\\
\\mu(d) = \\tilde{\\mu}\\cdot\\sigma_0(d+\\theta)^{-\\eta}+\\tau, \\\\
![
\\sigma(d)=\\sigma\_0(d+\\theta)^{-\\eta\_2}+\\tau, \\\\
\\mu(d) = \\tilde{\\mu}\\cdot\\sigma\_0(d+\\theta)^{-\\eta}+\\tau, \\\\
\\xi(d) = \\text{const.}
")
](https://latex.codecogs.com/png.latex?%0A%5Csigma%28d%29%3D%5Csigma_0%28d%2B%5Ctheta%29%5E%7B-%5Ceta_2%7D%2B%5Ctau%2C%20%5C%5C%0A%5Cmu%28d%29%20%3D%20%5Ctilde%7B%5Cmu%7D%5Ccdot%5Csigma_0%28d%2B%5Ctheta%29%5E%7B-%5Ceta%7D%2B%5Ctau%2C%20%20%5C%5C%0A%5Cxi%28d%29%20%3D%20%5Ctext%7Bconst.%7D%20%0A "
\sigma(d)=\sigma_0(d+\theta)^{-\eta_2}+\tau, \\
\mu(d) = \tilde{\mu}\cdot\sigma_0(d+\theta)^{-\eta}+\tau, \\
\xi(d) = \text{const.}
")
The function `gev.d.fit` provides the options:
- `theta_zero = TRUE` ![\\theta
= 0](https://latex.codecogs.com/png.latex?%5Ctheta%20%3D%200
"\\theta = 0")
- `eta2_zero = TRUE` (default) ![\\eta\_2 =
\\eta](https://latex.codecogs.com/png.latex?%5Ceta_2%20%3D%20%5Ceta
"\\eta_2 = \\eta")
- `tau_zero = TRUE` (default) ![\\tau
= 0](https://latex.codecogs.com/png.latex?%5Ctau%20%3D%200
"\\tau = 0")
- `theta_zero = TRUE` ![\\theta = 0](https://latex.codecogs.com/png.latex?%5Ctheta%20%3D%200 "\theta = 0")
- `eta2_zero = TRUE` (default) ![\\eta\_2 = \\eta](https://latex.codecogs.com/png.latex?%5Ceta_2%20%3D%20%5Ceta "\eta_2 = \eta")
- `tau_zero = TRUE` (default) ![\\tau = 0](https://latex.codecogs.com/png.latex?%5Ctau%20%3D%200 "\tau = 0")
resulting in the following features for IDF-curves:
- simple scaling: using only parameters ![\\tilde{\\mu}, \\sigma\_0,
\\xi,
\\eta](https://latex.codecogs.com/png.latex?%5Ctilde%7B%5Cmu%7D%2C%20%5Csigma_0%2C%20%5Cxi%2C%20%5Ceta
"\\tilde{\\mu}, \\sigma_0, \\xi, \\eta")
- curvature for small durations: allowing ![\\theta
\\neq 0](https://latex.codecogs.com/png.latex?%5Ctheta%20%5Cneq%200
"\\theta \\neq 0") (default)
- multi-scaling: allowing ![\\eta\_2 \\neq
\\eta](https://latex.codecogs.com/png.latex?%5Ceta_2%20%5Cneq%20%5Ceta
"\\eta_2 \\neq \\eta")
- flattening for long durations: allowing ![\\tau
\\neq 0](https://latex.codecogs.com/png.latex?%5Ctau%20%5Cneq%200
"\\tau \\neq 0").
- simple scaling: using only parameters ![\\tilde{\\mu}, \\sigma\_0, \\xi, \\eta](https://latex.codecogs.com/png.latex?%5Ctilde%7B%5Cmu%7D%2C%20%5Csigma_0%2C%20%5Cxi%2C%20%5Ceta "\tilde{\mu}, \sigma_0, \xi, \eta")
- curvature for small durations: allowing ![\\theta \\neq 0](https://latex.codecogs.com/png.latex?%5Ctheta%20%5Cneq%200 "\theta \neq 0") (default)
- multi-scaling: allowing ![\\eta\_2 \\neq \\eta](https://latex.codecogs.com/png.latex?%5Ceta_2%20%5Cneq%20%5Ceta "\eta_2 \neq \eta")
- flattening for long durations: allowing ![\\tau \\neq 0](https://latex.codecogs.com/png.latex?%5Ctau%20%5Cneq%200 "\tau \neq 0").
Example:
......@@ -189,7 +147,7 @@ set.seed(42)
# durations
ds <- 1/60*2^(seq(0,13,1))
# random data for each duration
xdat <- sapply(ds,rgev.d,n = 20,mut = 2,sigma0 =3,xi = 0.2,theta = 0.1,eta = 0.6,tau = 0.1,eta2 = 0.8)
xdat <- sapply(ds,rgev.d,n = 20,mut = 2,sigma0 =3,xi = 0.2,theta = 0.1,eta = 0.6,tau = 0.1,eta2 = 0.2)
# transform to data.frame
example <- data.frame(xdat=as.numeric(xdat),ds=rep(ds,each=dim(xdat)[1]))
......@@ -223,11 +181,36 @@ for(i.fit in 1:length(all.fits)){
}
for(i.p in 1:length(idf.probs)){
# plotting IDF curves for each model (different colors) and probability (different lty)
print(gev.d.params(all.fits[[i.fit]]))
IDF.plot(1/60*2^(seq(0,13,0.5)),gev.d.params(all.fits[[i.fit]]),probs = idf.probs[i.p]
,add = TRUE,legend = FALSE,lty = i.p,cols = fit.cols[i.fit])
}
mtext(fit.labels[i.fit],3,-1.25)
}
#> mut sigma0 xi theta eta eta2 tau
#> 1 1.754974 2.897875 -0.004559432 0 0.4978471 0 0
#> mut sigma0 xi theta eta eta2 tau
#> 1 1.754974 2.897875 -0.004559432 0 0.4978471 0 0
#> mut sigma0 xi theta eta eta2 tau
#> 1 1.754974 2.897875 -0.004559432 0 0.4978471 0 0
#> mut sigma0 xi theta eta eta2 tau
#> 1 1.782777 3.187756 -0.01039414 0.02949747 0.5376237 0 0
#> mut sigma0 xi theta eta eta2 tau
#> 1 1.782777 3.187756 -0.01039414 0.02949747 0.5376237 0 0
#> mut sigma0 xi theta eta eta2 tau
#> 1 1.782777 3.187756 -0.01039414 0.02949747 0.5376237 0 0
#> mut sigma0 xi theta eta eta2 tau
#> 1 1.896368 2.826447 0.1358611 0.03034914 0.4971133 0.1577612 0
#> mut sigma0 xi theta eta eta2 tau
#> 1 1.896368 2.826447 0.1358611 0.03034914 0.4971133 0.1577612 0
#> mut sigma0 xi theta eta eta2 tau
#> 1 1.896368 2.826447 0.1358611 0.03034914 0.4971133 0.1577612 0
#> mut sigma0 xi theta eta eta2 tau
#> 1 2.038391 2.712067 0.1469354 0.09235748 0.6237206 0.2642782 0.1287691
#> mut sigma0 xi theta eta eta2 tau
#> 1 2.038391 2.712067 0.1469354 0.09235748 0.6237206 0.2642782 0.1287691
#> mut sigma0 xi theta eta eta2 tau
#> 1 2.038391 2.712067 0.1469354 0.09235748 0.6237206 0.2642782 0.1287691
legend('topright',lty=rev(1:3),legend = rev(idf.probs),title = 'p-quantile')
```
......
......@@ -4,7 +4,7 @@
\alias{dgev.d}
\title{d-GEV probability density function}
\usage{
dgev.d(q, mut, sigma0, xi, theta, eta, d, eta2 = NULL, tau = 0, ...)
dgev.d(q, mut, sigma0, xi, theta, eta, d, eta2 = 0, tau = 0, ...)
}
\arguments{
\item{q}{vector of quantiles}
......@@ -18,7 +18,7 @@ shape parameter \eqn{\xi}.}
\item{d}{positive numeric value, giving duration}
\item{eta2}{numeric value, giving a second duration exponent \eqn{\eta_2} (needed for multiscaling). Default: NULL, treated as \eqn{\eta_2=\eta}.}
\item{eta2}{numeric value, giving a second duration exponent \eqn{\eta_{2}} (needed for multiscaling). Default: \eqn{\eta_{2}=0}.}
\item{tau}{numeric value, giving intensity offset \eqn{\tau} (defining flattening of the IDF curve). Default: \eqn{\tau=0}.}
......
man/figures/README-annualmax-1.png

9.68 KB | W: | H:

man/figures/README-annualmax-1.png

9.7 KB | W: | H:

man/figures/README-annualmax-1.png
man/figures/README-annualmax-1.png
man/figures/README-annualmax-1.png
man/figures/README-annualmax-1.png
  • 2-up
  • Swipe
  • Onion skin
man/figures/README-fit-1.png

44.6 KB | W: | H:

man/figures/README-fit-1.png

44.7 KB | W: | H:

man/figures/README-fit-1.png
man/figures/README-fit-1.png
man/figures/README-fit-1.png
man/figures/README-fit-1.png
  • 2-up
  • Swipe
  • Onion skin
man/figures/README-fit-2.png

22.8 KB | W: | H:

man/figures/README-fit-2.png

23.6 KB | W: | H:

man/figures/README-fit-2.png
man/figures/README-fit-2.png
man/figures/README-fit-2.png
man/figures/README-fit-2.png
  • 2-up
  • Swipe
  • Onion skin
man/figures/README-idf-1.png

35.5 KB | W: | H:

man/figures/README-idf-1.png

35.7 KB | W: | H:

man/figures/README-idf-1.png
man/figures/README-idf-1.png
man/figures/README-idf-1.png
man/figures/README-idf-1.png
  • 2-up
  • Swipe
  • Onion skin
......@@ -109,7 +109,7 @@ For details on the d-GEV and the parameter definitions, see \link{IDF-package}.
# xi = 0.5
# theta = 0
# eta = 0.5
# eta2 = 0.5
# eta2 = 0
# tau = 0
data('example',package ='IDF')
......
......@@ -14,7 +14,7 @@ gev.d.lik(
eta,
log = FALSE,
tau = 0,
eta2 = NULL
eta2 = 0
)
}
\arguments{
......
......@@ -4,7 +4,7 @@
\alias{pgev.d}
\title{d-GEV cumulative distribution function}
\usage{
pgev.d(q, mut, sigma0, xi, theta, eta, d, tau = 0, eta2 = NULL, ...)
pgev.d(q, mut, sigma0, xi, theta, eta, d, tau = 0, eta2 = 0, ...)
}
\arguments{
\item{q}{vector of quantiles}
......@@ -19,7 +19,7 @@ pgev.d(q, mut, sigma0, xi, theta, eta, d, tau = 0, eta2 = NULL, ...)
\item{tau}{numeric value, giving intensity offset \eqn{\tau} (defining flattening of the IDF curve). Default: \eqn{\tau=0}.}
\item{eta2}{numeric value, giving a second duration exponent \eqn{\eta_2} (needed for multiscaling). Default: NULL, treated as \eqn{\eta_2=\eta}.}
\item{eta2}{numeric value, giving a second duration exponent \eqn{\eta_{2}} (needed for multiscaling). Default: \eqn{\eta_{2}=0}.}
\item{...}{additional parameters passed to \code{\link[evd]{pgev}}}
}
......
......@@ -4,7 +4,7 @@
\alias{qgev.d}
\title{d-GEV quantile function}
\usage{
qgev.d(p, mut, sigma0, xi, theta, eta, d, tau = 0, eta2 = NULL, ...)
qgev.d(p, mut, sigma0, xi, theta, eta, d, tau = 0, eta2 = 0, ...)
}
\arguments{
\item{p}{vector of probabilities}
......@@ -19,7 +19,7 @@ qgev.d(p, mut, sigma0, xi, theta, eta, d, tau = 0, eta2 = NULL, ...)
\item{tau}{numeric value, giving intensity offset \eqn{\tau} (defining flattening of the IDF curve). Default: \eqn{\tau=0}.}
\item{eta2}{numeric value, giving a second duration exponent \eqn{\eta_2} (needed for multiscaling). Default: NULL, treated as \eqn{\eta_2=\eta}.}
\item{eta2}{numeric value, giving a second duration exponent \eqn{\eta_{2}} (needed for multiscaling). Default: \eqn{\eta_{2}=0}.}
\item{...}{additional parameters passed to \code{\link[evd]{qgev}}}
}
......
......@@ -4,7 +4,7 @@
\alias{rgev.d}
\title{Generation of random variables from d-GEV}
\usage{
rgev.d(n, mut, sigma0, xi, theta, eta, d, tau = 0, eta2 = NULL)
rgev.d(n, mut, sigma0, xi, theta, eta, d, tau = 0, eta2 = 0)
}
\arguments{
\item{n}{number of random variables per duration}
......@@ -19,7 +19,7 @@ rgev.d(n, mut, sigma0, xi, theta, eta, d, tau = 0, eta2 = NULL)
\item{tau}{numeric value, giving intensity offset \eqn{\tau} (defining flattening of the IDF curve). Default: \eqn{\tau=0}.}
\item{eta2}{numeric value, giving a second duration exponent \eqn{\eta_2} (needed for multiscaling). Default: NULL, treated as \eqn{\eta_2=\eta}.}
\item{eta2}{numeric value, giving a second duration exponent \eqn{\eta_{2}} (needed for multiscaling). Default: \eqn{\eta_{2}=0}.}
}
\value{
list containing vectors of random variables.
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment