Commit 203c03bf authored by jana ulrich's avatar jana ulrich
Browse files

implemented analytical likelihood gradient -> faster optimization

parent 313e7813
^.*\.Rproj$
^\.Rproj\.user$
^\.Rhistory$
^\IDF.Rproj$
^cran-comments\.md$
^README\.Rmd$
^CRAN-RELEASE$
......@@ -2,7 +2,7 @@ Package: IDF
Type: Package
Title: Estimation and Plotting of IDF Curves
Version: 2.1.0.9000
Date: 2021-05-12
Date: 2021-09-02
Authors@R: c(person("Jana", "Ulrich", email = "jana.ulrich@met.fu-berlin.de", role = c("aut", "cre")),
person("Laura","Mack", email= "laura.mack@fu-berlin.de",role="ctb"),
person("Oscar E.","Jurado", email= "jurado@zedat.fu-berlin.de",role="ctb"),
......
# IDF 2.1.0.9000
- added confidence intervals for diagnostic plots (gev.d.diag)
- changed defintion of eta2
- implemented analytical likelihood gradient for faster oprimization
# IDF 2.1.0
......
......@@ -48,7 +48,7 @@
#' \item{se}{numeric vector giving the standard errors for the MLE's (in the same order)}
#' \item{trans}{A logical indicator for a non-stationary fit.}
#' \item{model}{A list with components mutl, sigma0l, xil, thetal and etal. If requested, contains also eta2l and taul}
#' \item{link}{A character vector giving inverse link functions.}
#' \item{link}{A character vector giving link functions.}
#' \item{conv}{The convergence code, taken from the list returned by \code{\link{optim}}.
#' A zero indicates successful convergence.}
#' \item{data}{data is standardized to standard Gumbel.}
......@@ -78,12 +78,11 @@
#' gev.d.fit(xdat=example$dat,ds = example$d,ydat=as.matrix(example[,c('cov1','cov2')])
#' ,mutl=c(1,2),sigma0l=1)
gev.d.fit<-
function(xdat, ds, ydat = NULL, mutl = NULL, sigma0l = NULL, xil = NULL, thetal = NULL, etal = NULL, taul = NULL, eta2l = NULL,
gev.d.fit <- function(xdat, ds, ydat = NULL, mutl = NULL, sigma0l = NULL, xil = NULL, thetal = NULL, etal = NULL, taul = NULL, eta2l = NULL,
mutlink = make.link("identity"), sigma0link = make.link("identity"), xilink = make.link("identity"),
thetalink = make.link("identity"), etalink = make.link("identity"), taulink = make.link("identity"), eta2link = make.link("identity"),
init.vals = NULL, theta_zero = FALSE, tau_zero=TRUE, eta2_zero=TRUE,
show = TRUE, method = "Nelder-Mead", maxit = 10000, ...)
show = TRUE, method = "Nelder-Mead", maxit = 10000,...)
{
if (length(xdat) != length(ds)) {
stop(paste0('The length of xdat is ',length(xdat),', but the length of ds is ',length(ds),'.'))
......@@ -221,57 +220,99 @@ gev.d.fit<-
if (!eta2_zero) {init <- c(init,eta2init)} # add eta2 init (in case)
if (!tau_zero) {init <- c(init,tauinit)} # add tau init (in case)
# function to calculate neg log-likelihood:
# functions to calculate neg log-likelihood and gradient:
gev.lik <- function(a) {
# computes neg log lik of d-gev model
mu <- mutlink$linkinv(mumat %*% (a[1:npmu]))
sigma <- sigma0link$linkinv(sigmat %*% (a[seq(npmu + 1, length = npsc)]))
xi <- xilink$linkinv(shmat %*% (a[seq(npmu + npsc + 1, length = npsh)]))
# Next line will set the theta likelihood as non-existent in case user requested it. (same for tau)
if(!theta_zero) {theta <- thetalink$linkinv(thmat %*% (a[seq(npmu + npsc + npsh + 1, length = npth)]))}
theta <-if(!theta_zero){thetalink$linkinv(thmat %*% (a[seq(npmu + npsc + npsh + 1, length = npth)]))}else{0}
eta <- etalink$linkinv(etmat %*% (a[seq(npmu + npsc + npsh + npth + 1, length = npet)]))
if(!eta2_zero) {eta2 <- eta2link$linkinv( e2mat %*% (a[seq(npmu + npsc + npsh + npth + npet + 1, length = npe2)]))}
if(!tau_zero) {tau <- taulink$linkinv( tamat %*% (a[seq(npmu + npsc + npsh + npth + npet + npe2 + 1, length = npta)]))}
eta2 <-if(!eta2_zero){eta2link$linkinv( e2mat %*% (a[seq(npmu + npsc + npsh + npth + npet + 1, length = npe2)]))}else{eta}
tau <-if(!tau_zero){taulink$linkinv( tamat %*% (a[seq(npmu + npsc + npsh + npth + npet + npe2 + 1, length = npta)]))}else{0}
ifelse(!theta_zero, ds.t <- ds+theta, ds.t <- ds) #Don't use theta if user requested not to have it.
#ifelse(!tau_zero, sigma.d <- sigma/(ds.t^eta)+tau, sigma.d <- sigma/(ds.t^eta)) # don't use tau if user requested not to have it.
if (tau_zero){ # don't use tau if user requested not to have it.
if (eta2_zero){ # don't use eta2
sigma.d <- sigma/(ds.t^eta)
mu.d <- mu*sigma.d
}else{ # use eta2 (and no tau)
sigma.d <- sigma/(ds.t^eta2)
mu.d <- mu*sigma/(ds.t^eta)
}
}else{ # use tau
if (eta2_zero){ # don't use eta2
sigma.d <- sigma/(ds.t^eta)+tau
mu.d <- mu*sigma.d
}else{ # use eta2 (and tau)
sigma.d <- sigma/(ds.t^eta2)+tau
mu.d <- mu*(sigma/(ds.t^eta)+tau)
}
}
#sigma.d <- sigma/(ds.t^eta) # old
y = (xdat - mu.d) / sigma.d # new
#y = (xdat - mu*sigma.d) / sigma.d # derivation
#y <- xdat/sigma.d - mu # old
ds.t <- ds+theta
sigma.d <- sigma/(ds.t^eta2)+tau
mu.d <- mu*(sigma/(ds.t^eta)+tau)
y = (xdat - mu.d) / sigma.d
y <- 1 + xi * y
if(!theta_zero) {if(any(theta < 0)) {return(10^6)} } # check definition condition for theta
if(any(eta <= 0) || any(sigma.d <= 0) || any(y <= 0)) return(10^6)
if(!tau_zero) {if(any(tau < 0)) {return(10^6)} } # check definition condition for tau.
if(!eta2_zero) {if(any(eta2 < 0)) {return(10^6)} } # check definition condition for eta2.
sum(log(sigma.d)) + sum(y^(-1/xi)) + sum(log(y) * (1/xi + 1))
return(sum(log(sigma.d)) + sum(y^(-1/xi)) + sum(log(y) * (1/xi + 1)))
}
gev.grad <- function(a){
# computes gradient of neg log lik of d-gev model
mu <- mutlink$linkinv(mumat %*% (a[1:npmu]))
sigma <- sigma0link$linkinv(sigmat %*% (a[seq(npmu + 1, length = npsc)]))
xi <- xilink$linkinv(shmat %*% (a[seq(npmu + npsc + 1, length = npsh)]))
theta <-if(!theta_zero){thetalink$linkinv(thmat %*% (a[seq(npmu + npsc + npsh + 1, length = npth)]))}else{0}
eta <- etalink$linkinv(etmat %*% (a[seq(npmu + npsc + npsh + npth + 1, length = npet)]))
eta2 <-if(!eta2_zero){eta2link$linkinv( e2mat %*% (a[seq(npmu + npsc + npsh + npth + npet + 1, length = npe2)]))}else{eta}
tau <-if(!tau_zero){taulink$linkinv( tamat %*% (a[seq(npmu + npsc + npsh + npth + npet + npe2 + 1, length = npta)]))}else{0}
sigma.d <- sigma/((ds + theta)^eta2) + tau
mu.d <- mu * (sigma/((ds + theta)^eta) + tau)
y <- 1 + xi * (xdat - (mu.d))/(sigma.d)
dnll.mu.out <- -(xi * (sigma/((ds + theta)^eta) + tau)/(sigma.d)/(y) * (1/xi + 1) + (y)^((-1/xi) - 1) * ((-1/xi) * (xi * (sigma/((ds + theta)^eta) + tau)/(sigma.d))))
dnll.mu <- apply(mumat,2,function(x)sum(dnll.mu.out*mutlink$mu.eta(mumat %*% (a[1:npmu]))*x))
dnll.sigma.out <- 1/((ds + theta)^eta2)/(sigma.d) - (y)^((-1/xi) - 1) *
((-1/xi) * (xi * (mu * (1/((ds + theta)^eta)))/(sigma.d) + xi * (xdat - (mu.d)) * (1/((ds + theta)^eta2))/(sigma.d)^2)) -
(xi * (mu * (1/((ds + theta)^eta)))/(sigma.d) + xi * (xdat - (mu.d)) * (1/((ds + theta)^eta2))/(sigma.d)^2)/(y) * (1/xi + 1)
dnll.sigma <- apply(sigmat,2,function(x)sum(dnll.sigma.out*sigma0link$mu.eta(sigmat %*% (a[seq(npmu + 1, length = npsc)]))*x))
dnll.xi.out <- (y)^((-1/xi) - 1) * ((-1/xi) * ((xdat - (mu.d))/(sigma.d))) + (y)^(-1/xi) * (log((y)) * (1/xi^2)) +
((xdat - (mu.d))/(sigma.d)/(y) * (1/xi + 1) - log(y) * (1/xi^2))
dnll.xi <- apply(shmat,2,function(x)sum(dnll.xi.out*xilink$mu.eta(shmat %*% (a[seq(npmu + npsc + 1, length = npsh)]))*x))
if(!theta_zero){
dnll.theta.out <- (y)^((-1/xi) - 1) * ((-1/xi) * (xi * (mu * (sigma * ((ds + theta)^(eta - 1) * eta)/((ds + theta)^eta)^2))/
(sigma.d) + xi * (xdat - (mu.d)) * (sigma * ((ds + theta)^(eta2 - 1) * eta2)/((ds + theta)^eta2)^2)/(sigma.d)^2)) -
sigma * ((ds + theta)^(eta2 - 1) * eta2)/((ds + theta)^eta2)^2/(sigma.d) +
(xi * (mu * (sigma * ((ds + theta)^(eta - 1) * eta)/((ds + theta)^eta)^2))/(sigma.d) + xi * (xdat - (mu.d)) *
(sigma * ((ds + theta)^(eta2 - 1) * eta2)/((ds + theta)^eta2)^2)/(sigma.d)^2)/(y) * (1/xi + 1)
dnll.theta <- apply(thmat,2,function(x)sum(dnll.theta.out*thetalink$mu.eta(thmat %*% (a[seq(npmu + npsc + npsh + 1, length = npth)]))*x))
}
if(eta2_zero){
dnll.eta.out <- (y)^((-1/xi) - 1) * ((-1/xi) * (xi * (mu * (sigma * ((ds + theta)^eta * log((ds + theta)))/((ds + theta)^eta)^2))/(sigma.d) + xi * (xdat - (mu.d)) *
(sigma * ((ds + theta)^eta * log((ds + theta)))/((ds + theta)^eta)^2)/(sigma.d)^2)) - sigma *
((ds + theta)^eta * log((ds + theta)))/((ds + theta)^eta)^2/(sigma.d) +
(xi * (mu * (sigma * ((ds + theta)^eta * log((ds + theta)))/((ds + theta)^eta)^2))/(sigma.d) + xi * (xdat - (mu.d)) *
(sigma * ((ds + theta)^eta * log((ds + theta)))/((ds + theta)^eta)^2)/(sigma.d)^2)/(y) * (1/xi + 1)
dnll.eta <- apply(etmat,2,function(x)sum(dnll.eta.out*etalink$mu.eta(etmat %*% (a[seq(npmu + npsc + npsh + npth + 1, length = npet)]))*x))
}else{
dnll.eta.out <- (y)^((-1/xi) - 1) * ((-1/xi) * (xi * (mu * (sigma * ((ds + theta)^eta * log((ds + theta)))/((ds + theta)^eta)^2))/(sigma.d))) + xi *
(mu * (sigma * ((ds + theta)^eta * log((ds + theta)))/((ds + theta)^eta)^2))/(sigma.d)/(y) * (1/xi + 1)
dnll.eta <- apply(etmat,2,function(x)sum(dnll.eta.out*etalink$mu.eta(etmat %*% (a[seq(npmu + npsc + npsh + npth + 1, length = npet)]))*x))
dnll.eta2.out <- (y)^((-1/xi) - 1) * ((-1/xi) * (xi * (xdat - (mu.d)) * (sigma * ((ds + theta)^eta2 * log((ds + theta)))/((ds + theta)^eta2)^2)/(sigma.d)^2)) -
sigma * ((ds + theta)^eta2 * log((ds + theta)))/((ds + theta)^eta2)^2/(sigma.d) + xi * (xdat - (mu.d)) *
(sigma * ((ds + theta)^eta2 * log((ds + theta)))/((ds + theta)^eta2)^2)/(sigma.d)^2/(y) * (1/xi + 1)
dnll.eta2 <- apply(e2mat,2,function(x)sum(dnll.eta2.out*eta2link$mu.eta( e2mat %*% (a[seq(npmu + npsc + npsh + npth + npet + 1, length = npe2)]))*x))
}
if(!tau_zero){
dnll.tau.out <- 1/(sigma.d) - (y)^((-1/xi) - 1) * ((-1/xi) * (xi * mu/(sigma.d) + xi * (xdat - (mu.d))/(sigma.d)^2)) -
(xi * mu/(sigma.d) + xi * (xdat - (mu.d))/(sigma.d)^2)/(y) * (1/xi + 1)
dnll.tau <- apply(tamat,2,function(x)sum(dnll.tau.out*taulink$mu.eta(tamat %*% (a[seq(npmu + npsc + npsh + npth + npet + npe2 + 1, length = npta)]))*x))
}
grad.nll <- c(dnll.mu,dnll.sigma,dnll.xi,if(!theta_zero){dnll.theta},dnll.eta,if(!eta2_zero){dnll.eta2},if(!tau_zero){dnll.tau})
if(any(theta < 0)) {return(rep(0,length(grad.nll)))} # check definition condition for theta
if(any(eta <= 0) || any(sigma.d <= 0) || any(y <= 0)) return(rep(0,length(grad.nll)))
if(any(tau < 0)) {return(rep(0,length(grad.nll)))} # check definition condition for tau.
if(any(eta2 <= 0)) {return(rep(0,length(grad.nll)))} # check definition condition for eta2.
return(grad.nll)
}
# finding minimum of log-likelihood:
x <- optim(init, gev.lik, hessian = TRUE, method = method,
x <- optim(init, gev.lik, hessian = TRUE, method = method, gr = gev.grad,
control = list(maxit = maxit, ...))
# saving output parameters:
......@@ -321,32 +362,11 @@ gev.d.fit<-
z$cov <- matrix(NA,length(z$mle),length(z$mle))
}
z$se <- sqrt(diag(z$cov)) # sqrt(digonal entries) = standart error of mle's
'if (!theta_zero) {#When theta parameter is returned (default)
if (!tau_zero){ # when tau is returned
z$vals <- cbind(mut, sc0, xi, theta, eta, tau)
}else{ # when tau is not returned
z$vals <- cbind(mut, sc0, xi, theta, eta)
}
} else {#When theta parameter is not returned, asked by user
if (!tau_zero){ # if tau is returned
z$vals <- cbind(mut, sc0, xi, eta, tau)
}else{ # if tau is not returned
z$vals <- cbind(mut, sc0, xi, eta)
}
}'
z$vals <- cbind(mut, sc0, xi, theta, eta, eta2, tau)
colnames(z$vals) <- c("mut","sigma0","xi","theta","eta","eta2","tau")
z$init.vals <- unlist(init.vals)
# transform eta2 to zero-oriented
z$init.vals[6] = z$init.vals[6] + z$init.vals[4]
'names2 = c("mut","sigma0","xi") # fixed part for set of names
if(!theta_zero){names2=c(names2,"theta")} # add theta (in case)
names2 = c(names2, "eta") # add eta (always)
if(!tau_zero){names2=c(names2, "tau")} # add tau (in case)
colnames(z$vals) <- names2'
colnames(z$vals) <- c("mut","sigma0","xi","theta","eta","eta2","tau")
z$ds <- ds
z$theta_zero <- theta_zero #Indicates if theta parameter was set to zero by user.
z$tau_zero <- tau_zero #Indicates if tau parameter was set to zero by user. (default)
......
......@@ -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,cl=TRUE)
gev.d.diag(fit,pch=1,ci=TRUE)
# parameter estimates
params <- gev.d.params(fit)
print(params)
......@@ -102,8 +102,7 @@ IDF.plot(durations,params,add=TRUE)
## 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).
This Example depicts the different features that can be used to model the IDF curves, see Fauer et. al (2021, https://hess.copernicus.org/preprints/hess-2021-334/).
Here we assume, that the block maxima of each duration can be modeled with the GEV distribution ($\xi\neq0$):
$$G(z;\mu,\sigma,\xi)=\exp \left\lbrace -\left[
1+\xi \left( \frac{z-\mu}{\sigma} \right)
......@@ -146,10 +145,12 @@ fit.simple <- gev.d.fit(xdat=example$xdat,ds = example$ds,theta_zero = TRUE,show
fit.theta <- gev.d.fit(xdat=example$xdat,ds = example$ds,show=FALSE)
fit.eta2 <- gev.d.fit(xdat=example$xdat,ds = example$ds,eta2_zero = FALSE,show=FALSE)
fit.tau <- gev.d.fit(xdat=example$xdat,ds = example$ds,eta2_zero = FALSE,tau_zero = FALSE,show=FALSE)
# group fits
all.fits <- list(simple=fit.simple,curvature=fit.theta,multiscaling=fit.eta2,flattening=fit.tau)
# compare parameter estimates:
print(t(sapply(all.fits,gev.d.params)))
### compare resulting idf-curves
# group fits
all.fits <- list(fit.simple,fit.theta,fit.eta2,fit.tau)
fit.cols <- c('red','purple','blue','darkgreen')
fit.labels <- c('simple-scaling','theta!=0','eta2!=eta','tau!=0')
# plotting probabilities
......@@ -171,7 +172,6 @@ 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>).
Installation
------------
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
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")
......@@ -22,12 +33,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
- Step 0: sample 20 years of example hourly precipitation data
``` r
set.seed(999)
......@@ -65,7 +76,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,cl=TRUE)
gev.d.diag(fit,pch=1,ci=TRUE)
```
<img src="man/figures/README-fit-1.png" width="100%" />
......@@ -103,10 +114,14 @@ IDF.plot(durations,params,add=TRUE)
<img src="man/figures/README-idf-1.png" width="100%" />
IDF Features
------------
## 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")):
This Example depicts the different features that can be used to model
the IDF curves, see Fauer et. al (2021,
<https://hess.copernicus.org/preprints/hess-2021-334/>). 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)
......@@ -128,16 +143,24 @@ where the GEV parameters depend on duration according to:
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:
......@@ -156,10 +179,22 @@ fit.simple <- gev.d.fit(xdat=example$xdat,ds = example$ds,theta_zero = TRUE,show
fit.theta <- gev.d.fit(xdat=example$xdat,ds = example$ds,show=FALSE)
fit.eta2 <- gev.d.fit(xdat=example$xdat,ds = example$ds,eta2_zero = FALSE,show=FALSE)
fit.tau <- gev.d.fit(xdat=example$xdat,ds = example$ds,eta2_zero = FALSE,tau_zero = FALSE,show=FALSE)
# group fits
all.fits <- list(simple=fit.simple,curvature=fit.theta,multiscaling=fit.eta2,flattening=fit.tau)
# compare parameter estimates:
print(t(sapply(all.fits,gev.d.params)))
#> mut sigma0 xi theta eta eta2
#> simple 1.754974 2.897875 -0.004559432 0 0.4978471 0
#> curvature 1.782777 3.187756 -0.01039414 0.02949747 0.5376237 0
#> multiscaling 1.896368 2.826447 0.1358611 0.03034914 0.4971133 0.1577612
#> flattening 2.038391 2.712067 0.1469354 0.09235748 0.6237206 0.2642782
#> tau
#> simple 0
#> curvature 0
#> multiscaling 0
#> flattening 0.1287691
### compare resulting idf-curves
# group fits
all.fits <- list(fit.simple,fit.theta,fit.eta2,fit.tau)
fit.cols <- c('red','purple','blue','darkgreen')
fit.labels <- c('simple-scaling','theta!=0','eta2!=eta','tau!=0')
# plotting probabilities
......@@ -181,36 +216,11 @@ 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')
```
......
......@@ -28,7 +28,7 @@ and describe the slope and curvature in the resulting IDF curves, respectively.
\item 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)}
\item A useful introduction to \strong{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
......
......@@ -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: \eqn{\eta_{2}=0}.}
\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.7 KB | W: | H:

man/figures/README-annualmax-1.png

9.68 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.7 KB | W: | H:

man/figures/README-fit-1.png

51.5 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

23.6 KB | W: | H:

man/figures/README-fit-2.png

23.2 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.7 KB | W: | H:

man/figures/README-idf-1.png

38.3 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
......@@ -83,7 +83,7 @@ duration offset and duration exponent, resp. If requested, contains also second
\item{se}{numeric vector giving the standard errors for the MLE's (in the same order)}
\item{trans}{A logical indicator for a non-stationary fit.}
\item{model}{A list with components mutl, sigma0l, xil, thetal and etal. If requested, contains also eta2l and taul}
\item{link}{A character vector giving inverse link functions.}
\item{link}{A character vector giving link functions.}
\item{conv}{The convergence code, taken from the list returned by \code{\link{optim}}.
A zero indicates successful convergence.}
\item{data}{data is standardized to standard Gumbel.}
......
......@@ -19,7 +19,7 @@ pgev.d(q, mut, sigma0, xi, theta, eta, d, tau = 0, eta2 = 0, ...)
\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: \eqn{\eta_{2}=0}.}
\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}}}
}
......
......@@ -19,7 +19,7 @@ qgev.d(p, mut, sigma0, xi, theta, eta, d, tau = 0, eta2 = 0, ...)
\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: \eqn{\eta_{2}=0}.}
\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}}}
}
......
......@@ -19,7 +19,7 @@ rgev.d(n, mut, sigma0, xi, theta, eta, d, tau = 0, eta2 = 0)
\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: \eqn{\eta_{2}=0}.}
\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