Skip to content
GitLab
Commit f7c3a7d0 authored by Jana Ulrich's avatar Jana Ulrich
Browse files

Tidied up. Moved functions to seperate R-files for better clarity. Removed obsolete functions.

parent efbb8850
Expand all Hide whitespace changes
Inline Side-by-side
.Rhistory 0 → 100644
View file @ f7c3a7d0
(requireNamespace('IDF',quietly = TRUE))
?IDF.short
IDF.short(int.vec = test.data, durs = rep(c(1,2,3),each=100), n.y = 100)
install.packages("~/Downloads/IDF-covariates-6c1b09f1f9fc97a0b518b9ae025ba535095f63ea.tar.gz", repos = NULL, type = "source")
library(ismev)
library(IDF)
test.data <- matrix(NA,ncol=3,nrow=100)
for(i in 1:3) test.data[,i] <- rgev.d(n = 100,d=i)
test.data <- as.vector(test.data)
?rgev.d
for(i in 1:3) test.data[,i] <- rgev.d(n = 100,d=i)
?rgev.d
?rgev
?IDF.agg
IDF.agg(test.data,agg.lev=c(1,2,3))
detach("package:IDF", unload=TRUE)
remove.packages("IDF", lib="~/R/x86_64-pc-linux-gnu-library/3.3")
library(IDF)
install.packages("~/Downloads/IDF-covariates-6c1b09f1f9fc97a0b518b9ae025ba535095f63ea.tar.gz", repos = NULL, type = "source")
library(IDF)
IDF.agg(test.data,agg.lev=c(1,2,3))
detach("package:IDF", unload=TRUE)
remove.packages("IDF", lib="~/R/x86_64-pc-linux-gnu-library/3.3")
library(IDF)
install.packages("~/Downloads/IDF-covariates-6c1b09f1f9fc97a0b518b9ae025ba535095f63ea(1).tar.gz", repos = NULL, type = "source")
library(ismev)
library(IDF)
IDF.agg(test.data,agg.lev=c(1,2,3))
?IDF.short
detach("package:IDF", unload=TRUE)
remove.packages("IDF", lib="~/R/x86_64-pc-linux-gnu-library/3.3")
?IDF.short
?install.packages
devtools::install_git('https://gitlab.met.fu-berlin.de/Rpackages/IDF/tree/covariates')
devtools::install_git('https://gitlab.met.fu-berlin.de/Rpackages/IDF')
devtools::install_git('https://gitlab.met.fu-berlin.de/Rpackages/IDF/tree/covariates.git')
DESCRIPTION
View file @ f7c3a7d0
Package: IDF
Type: Package
Title: Estimation and Plotting of IDF Curves
Version: 1.3
Date: 2018-02-06
Author: Christoph Ritschel, Carola Detring, Sarah Joedicke
Maintainer: Christoph Ritschel <christoph.ritschel@met.fu-berlin.de>
Description: Intensity-duration-frequency (IDF) curves are a widely used analysis-tool in hydrology to assess extreme values of precipitation [e.g. Mailhot et al., 2007, <doi:10.1016/j.jhydrol.2007.09.019>]. The package 'IDF' provides a function to read precipitation data from German weather service (DWD) 'webwerdis' <http://www.dwd.de/EN/ourservices/webwerdis/webwerdis.html> files and Berlin station data from 'Stadtmessnetz' <http://www.geo.fu-berlin.de/en/met/service/stadtmessnetz/index.html> files, and additionally IDF parameters can be estimated also from a given data.frame containing a precipitation time series. The data is aggregated to given levels yearly intensity maxima are calculated either for the whole year or given months. From these intensity maxima IDF parameters are estimated on the basis of a duration-dependent generalised extreme value distribution [Koutsoyannis et al., 1998, <doi:10.1016/S0022-1694(98)00097-3>]. IDF curves based on these estimated parameters can be plotted.
Depends:
stats4,
evd,
ismev
Version: 0.0.2
Date: 2019-02-06
Authors@R: c(person("Jana", "Ulrich", email = "jana.ulrich@fu-berlin.de", role = c("aut", "cre")),
person("Christoph", "Ritschel",email= "christoph.ritschel@met.fu-berlin.de", role = "aut"),
person("Carola", "Detring", role = "ctb"),
person("Sarah", "Joedicke", role = "ctb"))
Description: Intensity-duration-frequency (IDF) curves are a widely used analysis-tool
in hydrology to assess extreme values of precipitation
[e.g. Mailhot et al., 2007, <doi:10.1016/j.jhydrol.2007.09.019>].
The package 'IDF' provides functions to estimate IDF parameters for given
precipitation time series on the basis of a duration-dependent
generalised extreme value distribution
[Koutsoyannis et al., 1998, <doi:10.1016/S0022-1694(98)00097-3>].
Imports: stats4,
evd,
ismev,
zoo,
pbapply
License: GPL (>=2)
RoxygenNote: 5.0.1
Encoding: UTF-8
LazyData: true
RoxygenNote: 6.1.1
NAMESPACE
View file @ f7c3a7d0
export(IDF.read,IDF,IDF.short,IDF.plot,TS.acc,gev.d.fit,IDF.agg,rgev.d)
import(ismev,
stats4,
evd)
importFrom("grDevices", "rgb")
importFrom("graphics", "axis", "legend", "points")
importFrom("stats", "filter", "lm", "optim")
importFrom("utils", "read.csv2", "read.table", "str")
# Generated by roxygen2: do not edit by hand
export(IDF.agg)
export(IDF.plot)
export(dgev.d)
export(gev.d.diag)
export(gev.d.fit)
export(gev.d.params)
export(gev.d.rl)
export(pgev.d)
export(qgev.d)
export(rgev.d)
importFrom(evd,dgev)
importFrom(evd,pgev)
importFrom(evd,qgev)
importFrom(evd,rgev)
importFrom(grDevices,rainbow)
importFrom(grDevices,rgb)
importFrom(graphics,abline)
importFrom(graphics,axis)
importFrom(graphics,box)
importFrom(graphics,lines)
importFrom(graphics,par)
importFrom(graphics,plot)
importFrom(graphics,points)
importFrom(graphics,title)
importFrom(ismev,gev.fit)
importFrom(pbapply,pblapply)
importFrom(stats,lm)
importFrom(stats,optim)
importFrom(zoo,rollapply)
R/IDF.R
View file @ f7c3a7d0
This diff is collapsed.
R/d-gev.R 0 → 100644
View file @ f7c3a7d0
# This file contains the functions dgev.d, pgev.d, qgev.d, rgev.d for the duration-dependent-gev.
#### dgev.d() ####
#' Density function of duration dependent GEV distribution
#'
#' @param q vector of quantiles
#' @param mut,sigma0,xi numeric value, giving modified location, modified scale and shape parameter
#' @param theta numeric value, giving duration offset (defining curvature of the IDF curve)
#' @param eta numeric value, giving duration exponent (defining slope of the IDF curve)
#' @param d numeric value, giving duration
#' @param ... additional parameters passed to \code{\link[evd]{dgev}}
#'
#' @details The duration dependent GEV distribution is defined after
#' [Koutsoyannis et al., 1998]:
#' \deqn{ G(x)= exp[-{ 1+\xi(x/\sigma(d)-\mu_t) }^(-1/\xi)] }
#' with the duration dependent scale \eqn{\sigma(d)=\sigma_0/(d+\theta)^\eta} and
#' modified location parameter \eqn{\mu_t=\mu/\sigma(d)}.
#'
#' @return list containing vectors of density values for given quantiles
#' the first element of the list are the dens. values for the first given duration and so on
#'
#' @seealso \code{\link{pgev.d}}, \code{\link{qgev.d}}, \code{\link{rgev.d}}
#' @references Koutsoyannis et al., 1998, doi:10.1016/S0022-1694(98)00097-3
#'
#' @export
#' @importFrom evd dgev
#'
#' @examples
#' x <- seq(4,20,0.1)
#' dens <- dgev.d(q=x,mut=4,sigma0=2,xi=0,theta=0.1,eta=0.1,d=1:4)
#'
#' plot(x,dens[[1]],type='l',ylim = c(0,0.21),ylab = 'Probability Density')
#' for(i in 2:4){
#' lines(x,dens[[i]],lty=i)
#' }
#' legend('topright',title = 'duration',legend = 1:4,lty=1:4)
dgev.d <- function(q,mut,sigma0,xi,theta,eta,d,...) {
if(any(c(length(mut),length(sigma0),length(xi),length(theta),length(eta))>1)){
message('One of the parameters mut, sigma0, xi, theta, eta is a vector. ',
'This is not intended and might cause an error.')}
if(any(d+theta<=0)){
warning('Some shape parameters are negative,resulting from a negativ theta '
,theta, ' this will prododuce NAs.')}
# if denominator is negative NAs will be returned
densfun <- function(d){
if(d+theta<=0){return(rep(NA,length(q)))}else{
sigma.d <-sigma0/(d+theta)^eta
return(dgev(q,loc=mut*sigma.d,scale=sigma.d,shape=xi,...))}
}
# calculate for each duration in d
return(sapply(as.list(d),densfun,simplify = FALSE))
}
#### pgev.d() ####
#' Distribution function of duration dependent GEV distribution
#'
#' @param q vector of quantiles
#' @param mut,sigma0,xi numeric value, giving modified location, modified scale and shape parameter
#' @param theta numeric value, giving duration offset (defining curvature of the IDF curve)
#' @param eta numeric value, giving duration exponent (defining slope of the IDF curve)
#' @param d numeric value, giving duration
#' @param ... additional parameters passed to \code{\link[evd]{pgev}}
#'
#' @details The duration dependent GEV distribution is defined after
#' [Koutsoyannis et al., 1998]:
#' \deqn{ G(x)= exp[-{ 1+\xi(x/\sigma(d)-\mu_t) }^(-1/\xi)] }
#' with the duration dependent scale \eqn{\sigma(d)=\sigma_0/(d+\theta)^\eta} and
#' modified location parameter \eqn{\mu_t=\mu/\sigma(d)}.
#'
#' @return list containing vectors of probability values for given quantiles
#' the first element of the list are the prob. values for the first given duration and so on
#'
#' @seealso \code{\link{dgev.d}}, \code{\link{qgev.d}}, \code{\link{rgev.d}}
#' @references Koutsoyannis et al., 1998, doi:10.1016/S0022-1694(98)00097-3
#'
#' @export
#' @importFrom evd pgev
#'
#' @examples
#' x <- seq(4,20,0.1)
#' prob <- pgev.d(q=x,mu=4,sigma=2,xi=0,theta=0.1,eta=0.1,d=c(1,2))
pgev.d <- function(q,mut,sigma0,xi,theta,eta,d,...) {
if(any(c(length(mut),length(sigma0),length(xi),length(theta),length(eta))>1)){
message('One of the parameters mut, sigma0, xi, theta, eta is a vector. ',
'This is not intended and might cause an error.')}
if(any(d+theta<=0)){
warning('Some shape parameters are negative,resulting from a negativ theta '
,theta, ' this will prododuce NAs.')}
# if denominator is negative NAs will be returned
pfun <- function(d){
if(d+theta<=0){return(rep(NA,length(q)))}else{
sigma.d <-sigma0/(d+theta)^eta
return(pgev(q,loc=mut*sigma.d,scale=sigma.d,shape=xi,...))}
}
# calculate for each duration in d
return(sapply(as.list(d),pfun,simplify = FALSE))
}
#### qgev.d() ####
#' Quantile function of duration dependent GEV distribution
#'
#' @param p vector of probabilities
#' @param mut,sigma0,xi numeric value, giving modified location, modified scale and shape parameter
#' @param theta numeric value, giving duration offset (defining curvature of the IDF curve)
#' @param eta numeric value, giving duration exponent (defining slope of the IDF curve)
#' @param d numeric value, giving duration
#' @param ... additional parameters passed to \code{\link[evd]{qgev}}
#'
#' @details The duration dependent GEV distribution is defined after
#' [Koutsoyannis et al., 1998]:
#' \deqn{ G(x)= exp[-{ 1+\xi(x/\sigma(d)-\mu_t) }^(-1/\xi)] }
#' with the duration dependent scale \eqn{\sigma(d)=\sigma_0/(d+\theta)^\eta} and
#' modified location parameter \eqn{\mu_t=\mu/\sigma(d)}.
#'
#' @return list containing vectors of quantile values for given probabilities
#' the first element of the list are the q. values for the first given duration and so on
#'
#' @seealso \code{\link{pgev.d}}, \code{\link{dgev.d}}, \code{\link{rgev.d}}
#' @references Koutsoyannis et al., 1998, doi:10.1016/S0022-1694(98)00097-3
#'
#' @export
#' @importFrom evd qgev
#'
#' @examples
#' p <- c(0.5,0.9,0.99)
#' d <- 2^seq(0,4,0.1)
#' qs <- qgev.d(p=p,mu=4,sigma=2,xi=0,theta=0.1,eta=0.3,d=d)
#' qs <- simplify2array(qs)
#'
#' plot(d,qs[1,],ylim=c(3,20),type='l',log = 'xy',ylab='Intensity',xlab = 'Duration')
#' for(i in 2:3){
#' lines(d,qs[i,],lty=i)
#' }
#' legend(8,19,title = 'annual frequency \n of exceedance 1/T',
#' legend = 1-p,lty=1:3,bty = 'n')
qgev.d <- function(p,mut,sigma0,xi,theta,eta,d,...) {
if(any(c(length(mut),length(sigma0),length(xi),length(theta),length(eta))>1)){
message('One of the parameters mut, sigma0, xi, theta, eta is a vector. ',
'This is not intended and might cause an error.')}
if(any(d+theta<=0)){
warning('Some shape parameters are negative,resulting from a negativ theta '
,theta, ' this will prododuce NAs.')}
# if denominator is negative NAs will be returned
qfun <- function(d){
if(d+theta<=0){return(rep(NA,length(p)))}else{
sigma.d <-sigma0/(d+theta)^eta
return(qgev(p,loc=mut*sigma.d,scale=sigma.d,shape=xi,...))}
}
# calculate for each duration in d
return(sapply(as.list(d),qfun,simplify = FALSE))
}
#### rgev.d() ####
#' random generation of variables from duration dependent GEV distribution
#'
#' @param n number of observations.
#' @param mut,sigma0,xi numeric value, giving modified location, modified scale and shape parameter
#' @param theta numeric value, giving duration offset (defining curvature of the IDF curve)
#' @param eta numeric value, giving duration exponent (defining slope of the IDF curve)
#' @param d numeric value, giving duration
#'
#' @details The duration dependent GEV distribution is defined after
#' [Koutsoyannis et al., 1998]:
#' \deqn{ G(x)= exp[-{ 1+\xi(x/\sigma(d)-\mu_t) }^(-1/\xi)] }
#' with the duration dependent scale \eqn{\sigma(d)=\sigma_0/(d+\theta)^\eta} and
#' modified location parameter \eqn{\mu_t=\mu/\sigma(d)}.
#'
#' @return list containing vectors of random values
#' the first element of the list are the random values for the first given duration and so on
#'
#' @seealso \code{\link{pgev.d}}, \code{\link{qgev.d}}, \code{\link{dgev.d}}
#' @references Koutsoyannis et al., 1998, doi:10.1016/S0022-1694(98)00097-3
#'
#' @export
#' @importFrom evd rgev
#'
#' @examples
#' samp <- rgev.d(100,4,2,0,0,0.5,c(1,2))
#'
#' hist(samp[[1]],breaks = 10,col=rgb(1,0,0,0.5),freq = FALSE
#' ,ylim=c(0,0.3),xlab='x',main = 'd-GEV samples for two different durations')
#' hist(samp[[2]],breaks = 10,add=TRUE,col=rgb(0,0,1,0.5),freq = FALSE)
rgev.d <- function(n,mut,sigma0,xi,theta,eta,d) {
if(any(c(length(mut),length(sigma0),length(xi),length(theta),length(eta))>1)){
message('One of the parameters mut, sigma0, xi, theta, eta is a vector. ',
'This is not intended and might cause an error.')}
if(any(d+theta<=0)){
warning('Some shape parameters are negative,resulting from a negativ theta '
,theta, ' this will prododuce NAs.')}
# if denominator is negative NAs will be returned
rfun <- function(d){
if(d+theta<=0){return(rep(NA,n))}else{
sigma.d <-sigma0/(d+theta)^eta
return(rgev(n,loc=mut*sigma.d,scale=sigma.d,shape=xi))}
}
# calculate for each duration in d
return(sapply(as.list(d),rfun,simplify = FALSE))
}
R/gevdfit.R 0 → 100644
View file @ f7c3a7d0
# This file contains the functions:
# - gev.d.fit, gev.d.init for fitting
# - gev.d.diag for diagnostic plots
# - gev.d.params, gev.d.rl for calculation of parameters/ return levels
# and the documentation of the example data
#### gev.d.fit ####
#' @title Maximum-likelihood Fitting of the duration dependent GEV Distribution
#' @description Modified \code{\link[ismev]{gev.fit}} function for Maximum-likelihood fitting
#' for the duration dependent generalized extreme
#' value distribution, following Koutsoyiannis et al. (1988), including generalized linear
#' modelling of each parameter.
#' @param xdat A vector containing maxima for different durations.
#' This can be obtained from \code{\link{IDF.agg}}.
#' @param ds A vector of aggregation levels corresponding to the maxima in xdat.
#' @param ydat A matrix of covariates for generalized linear modelling of the parameters
#' (or NULL (the default) for stationary fitting). The number of rows should be the same as the
#' length of xdat.
#' @param mul,sigl,shl,thetal,etal Numeric vectors of integers, giving the columns of ydat that contain
#' covariates for generalized linear modelling of the parameters (or NULL (the default)
#' if the corresponding parameter is stationary).
#' Parameters are: modified location, scale_0, shape, duration offset, duration exponent repectively.
#' @param mulink,siglink,shlink,thetalink,etalink Inverse link functions for generalized linear
#' modelling of the parameters.
#' @param muinit,siginit,shinit,thetainit,etainit initial values as numeric of length
#' equal to total number of parameters
#' used to model the parameters. If NULL (the default) is given, initial parameters are obtained
#' internally by fitting the GEV seperately for each duration and applying a linear model to optain the
#' duration dependency of the location and shape parameter.
#' @param show Logical; if TRUE (the default), print details of the fit.
#' @param method The optimization method used in \code{\link{optim}}.
#' @param maxit The maximum number of iterations.
#' @param ... Other control parameters for the optimization.
#' @return A list containing the following components.
#' A subset of these components are printed after the fit.
#' If show is TRUE, then assuming that successful convergence is indicated,
#' the components nllh, mle and se are always printed.
#' \item{nllh}{single numeric giving the negative log-likelihood value.}
#' \item{mle}{numeric vector giving the MLE's for the modified location, scale_0, shape,
#' duration offset and duration exponent, resp.}
#' \item{se}{numeric vector giving the standard errors for the MLE's (in the same order).}
#' \item{trans}{An logical indicator for a non-stationary fit.}
#' \item{model}{A list with components mul, sigl, shl, thetal and etal.}
#' \item{link}{A character vector giving inverse 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 standart Gumbel.}
#' \item{cov}{The covariance matrix.}
#' @seealso \code{\link{dgev.d}}, \code{\link{IDF.agg}}, \code{\link{gev.fit}}, \code{\link{optim}}
#' @export
#' @importFrom stats optim
#'
#' @examples
#' # sampled random data from d-gev with covariates
#' # GEV parameters:
#' # mu = 4 + 0.2*cov1 +0.5*cov2
#' # sigma = 2+0.5*cov1
#' # xi = 0.5
#' # theta = 0
#' # eta = 0.5
#'
#' data('example',package ='IDF')
#'
#' gev.d.fit(xdat=example$dat,ds = example$d,ydat=as.matrix(example[,c('cov1','cov2')])
#' ,mul=c(1,2),sigl=1)
gev.d.fit<-
function(xdat, ds, ydat = NULL, mul = NULL, sigl = NULL, shl = NULL, thetal = NULL, etal = NULL,
mulink = identity, siglink = identity, shlink = identity, thetalink = identity, etalink = identity,
muinit = NULL, siginit = NULL, shinit = NULL, thetainit = NULL, etainit = NULL,
show = TRUE, method = "Nelder-Mead", maxit = 10000, ...)
{
#
# obtains mles etc for d-gev distn
#
# test for NA values:
if(any(is.na(xdat))) stop('xdat contains NA values. NA values need to be removed first.')
z <- list()
# number of parameters (betas) to estimate for each parameter:
npmu <- length(mul) + 1
npsc <- length(sigl) + 1
npsh <- length(shl) + 1
npth <- length(thetal) + 1
npet <- length(etal) + 1
z$trans <- FALSE # indicates if fit is non-stationary
# calculate initial values for mu.d, sigma_0, xi, eta using IDF.init: (thetainit=0)
init.vals <- gev.d.init(xdat,ds,ifelse(is.null(thetainit),0,thetainit[1]))
# generate covariates matrices:
if (is.null(mul)) {
mumat <- as.matrix(rep(1, length(xdat)))
if (is.null(muinit))
muinit <- init.vals$mu
}else {
z$trans <- TRUE
mumat <- cbind(rep(1, length(xdat)), ydat[, mul])
if (is.null(muinit))
muinit <- c(init.vals$mu, rep(0, length(mul)))
}
if (is.null(sigl)) {
sigmat <- as.matrix(rep(1, length(xdat)))
if (is.null(siginit))
siginit <- init.vals$sigma
}else {
z$trans <- TRUE
sigmat <- cbind(rep(1, length(xdat)), ydat[, sigl])
if (is.null(siginit))
siginit <- c(init.vals$sigma, rep(0, length(sigl)))
}
if (is.null(shl)) {
shmat <- as.matrix(rep(1, length(xdat)))
if (is.null(shinit))
shinit <- init.vals$xi
}else {
z$trans <- TRUE
shmat <- cbind(rep(1, length(xdat)), ydat[, shl])
if (is.null(shinit))
shinit <- c(init.vals$xi, rep(0, length(shl)))
}
if (is.null(thetal)) {
thmat <- as.matrix(rep(1, length(xdat)))
if (is.null(thetainit))
thetainit <- 0
}else {
z$trans <- TRUE
thmat <- cbind(rep(1, length(xdat)), ydat[, thetal])
if (is.null(thetainit))
thetainit <- c(0, rep(0, length(thetal)))
}
if (is.null(etal)) {
etmat <- as.matrix(rep(1, length(xdat)))
if (is.null(etainit))
etainit <- init.vals$eta
}else {
z$trans <- TRUE
etmat <- cbind(rep(1, length(xdat)), ydat[, etal])
if (is.null(etainit))
etainit <- c(init.vals$eta, rep(0, length(etal)))
}
z$model <- list(mul, sigl, shl, thetal, etal)
z$link <- deparse(substitute(c(mulink, siglink, shlink, thetalink, etalink)))
init <- c(muinit, siginit, shinit, thetainit, etainit)
# function to calculate neg log-likelihood:
gev.lik <- function(a) {
# computes neg log lik of d-gev model
mu <- mulink(mumat %*% (a[1:npmu]))
sigma <- siglink(sigmat %*% (a[seq(npmu + 1, length = npsc)]))
xi <- shlink(shmat %*% (a[seq(npmu + npsc + 1, length = npsh)]))
theta <- thetalink(thmat %*% (a[seq(npmu + npsc + npsh + 1, length = npth)]))
eta <- etalink(etmat %*% (a[seq(npmu + npsc + npsh + npth + 1, length = npet)]))
ds.t <- ds+theta
sigma.d <- sigma/(ds.t^eta)
y <- xdat/sigma.d - mu
y <- 1 + xi * y
if(any(eta <= 0) ||any(theta <= -0.01) || any(sigma.d <= 0) || any(y <= 0)) return(10^6)
sum(log(sigma.d)) + sum(y^(-1/xi)) + sum(log(y) * (1/xi + 1))
}
#####################################################################################
# derivations of nll after d-gev-parameters (for boosting):
# get parameters from covariates and a (vector containing predictors)
# mu <- mulink(mumat %*% (a[1:npmu]))
# sigma <- siglink(sigmat %*% (a[seq(npmu + 1, length = npsc)]))
# xi <- shlink(shmat %*% (a[seq(npmu + npsc + 1, length = npsh)]))
# theta <- thetalink(thmat %*% (a[seq(npmu + npsc + npsh + 1, length = npth)]))
# eta <- etalink(etmat %*% (a[seq(npmu + npsc + npsh + npth + 1, length = npet)]))
# xd <- xdat*(ds+theta)^eta
# y <- 1 + xi * (xd/sigma - mu)
#
# nll <- log(sigma/(ds+theta)^eta) + y^(-1/xi) + log(y) * (1/xi + 1)
# dnll.mu <- -xi/y
# dnll.sigma <- 1/(sigma+xi*xd/(1-mu*xi))
# dnll.xi <- 1/(xi+sigma/(xd-mu*sigma))
# dnll.theta <- - eta*sigma*(mu*xi-1)/(ds+theta)/(-xi*xd+mu*xi*sigma-sigma)
# dnll.eta <- -sigma*(mu*xi-1)*log(ds+theta)/(-xi*xd+mu*xi*sigma-sigma)
#####################################################################################
# finding minimum of log-likelihood:
x <- optim(init, gev.lik, hessian = TRUE, method = method,
control = list(maxit = maxit, ...))
# saving output parameters:
z$conv <- x$convergence
mut <- mulink(mumat %*% (x$par[1:npmu]))
sc0 <- siglink(sigmat %*% (x$par[seq(npmu + 1, length = npsc)]))
xi <- shlink(shmat %*% (x$par[seq(npmu + npsc + 1, length = npsh)]))
theta <- thetalink(thmat %*% (x$par[seq(npmu + npsc + npsh + 1, length = npth)]))
eta <- etalink(etmat %*% (x$par[seq(npmu + npsc + npsh + npth + 1, length = npet)]))
z$nllh <- x$value
# normalize data to standart gumbel:
sc.d <- sc0/((ds+theta)^eta)
z$data <- - log(as.vector((1 + xi * (xdat/sc.d-mut))^(-1/xi)))
z$mle <- x$par
z$cov <- solve(x$hessian) # invert hessian to get estimation on var-covar-matrix
z$se <- sqrt(diag(z$cov)) # sqrt(digonal entries) = standart error of mle's
z$vals <- cbind(mut, sc0, xi, theta, eta)
colnames(z$vals) <- c('mut','sigma0','xi','theta','eta')
z$ds <- ds
if(show) {
if(z$trans) # for nonstationary fit
print(z[c(2, 3, 4)]) # print model, link, conv
else print(z[4]) # for stationary fit print only conv
if(!z$conv) # if fit converged
print(z[c(5, 7, 9)]) # print nll, mle, se
}
class( z) <- "gev.d.fit"
invisible(z)
}
#### gev.d.init ####
# function to get initial values for gev.d.fit:
# obtain initial values
# by fitting every duration seperately
# possible ways to improve:
# take given initial values into accout, if there are any
# xi -> mean vs. median ... how do we improve that?
# mu_tilde -> is not very good for small sample sizes yet
# improved inital value for eta, by fitting both mu~d and sigma~d in log-log scale
#' @title get initial values for gev.d.fit
#' @description obtain initial values by fitting every duration seperately
#' @param xdat vector of maxima for differnt durations
#' @param ds vector of durations belonging to maxima in xdat
#' @param thetainit initial parameter for theta
#' @return list of initail values for mu_tilde, sigma_0, xi, eta
#' @importFrom stats lm
#' @importFrom ismev gev.fit
#' @keywords internal
gev.d.init <- function(xdat,ds,thetainit){
durs <- unique(ds)
mles <- matrix(NA, nrow=length(durs), ncol= 3)
for(i in 1:length(durs)){
mles[i,] <- gev.fit(xdat[ds==durs[i]],show = FALSE)$mle
}
# get values for sig0 and eta (also mu_0) from linear model in log-log scale
lmsig <- lm(log(mles[,2])~log(durs+thetainit))
lmmu <- lm(log(mles[,1])~log(durs+thetainit))
# sig0 <- exp Intercept
siginit <- exp(lmsig$coefficients[[1]])
# eta <- mean of negativ slopes
etainit <- mean(c(-lmsig$coefficients[[2]],-lmmu$coefficients[[2]]))
# mean of mu_d/sig_d
# could try:
# mu0/sig0 is also an estimate but needs to be weighted in mean
muinit <- mean(c(mles[,1]/mles[,2])) #exp(lmmu$coefficients[[1]])/exp(lmsig$coefficients[[1]])
# mean of shape parameters
shinit <- mean(mles[,3])
return(list(mu=muinit,sigma=siginit,xi=shinit,eta=etainit))
}
#### gev.d.diag ####
#' Diagnostic Plots for d-gev Models
#'
#' @description Produces diagnostic plots for d-gev models using
#' the output of the function \code{\link{gev.d.fit}}. Values for different durations can be plotted in
#' different colors of with different symbols.
#' @param fit object returned by \code{\link{gev.d.fit}}
#' @param subset an optional vector specifying a subset of observations to be used in the plot
#' @param cols optional either one value or vector of same length as \code{unique(durations)} to
#' specify the colors of plotting points.
#' The default uses the \code{rainbow} function.