Commit e99b28d8 authored by Jana Ulrich's avatar Jana Ulrich
Browse files

Merge branch 'branch-001' into 'covariates'

Branch 001

See merge request Rpackages/IDF!3
parents d44032f2 0aff0141
......@@ -17,9 +17,9 @@ Description: Intensity-duration-frequency (IDF) curves are a widely used analysi
Imports: stats4,
evd,
ismev,
zoo,
RcppRoll,
pbapply
License: GPL (>=2)
Encoding: UTF-8
LazyData: true
RoxygenNote: 6.1.1
RoxygenNote: 7.1.0
......@@ -10,6 +10,7 @@ export(gev.d.params)
export(pgev.d)
export(qgev.d)
export(rgev.d)
importFrom(RcppRoll,roll_sum)
importFrom(evd,dgev)
importFrom(evd,pgev)
importFrom(evd,qgev)
......@@ -30,5 +31,3 @@ importFrom(stats,lm)
importFrom(stats,make.link)
importFrom(stats,median)
importFrom(stats,optim)
importFrom(zoo,as.zoo)
importFrom(zoo,rollapplyr)
......@@ -11,7 +11,7 @@
#' the function \code{\link{gev.d.fit}}.
#'
#' @param data list of data.frames containing time series for every station.
#' The data.frame must have the columns 'date' and 'RR' unless onther names
#' The data.frame must have the columns 'date' and 'RR' unless other names
#' are specified in the parameter `names`. The column 'date' must contain strings with
#' standard date format.
#' @param ds numeric vector of aggregation durations.
......@@ -24,8 +24,8 @@
#' @param cl optional, number of cores to be used from \code{\link[pbapply]{pblapply}} for parallelization.
#'
#' @details If data contains stations with different time resolutions that need to be aggregated at
#' different durations, IDF.agg needs to be run seperately for the different groups of stations.
#' Afterwards he results can be joint together using `rbind`.
#' different durations, IDF.agg needs to be run separately for the different groups of stations.
#' Afterwards the results can be joint together using `rbind`.
#'
#' @return data.frame containing the annual intensity maxima [mm/h] in `$xdat`, the corresponding duration in `$ds`
#' and the station id or name in `$station`.
......@@ -34,8 +34,7 @@
#'
#' @export
#' @importFrom pbapply pbsapply
#' @importFrom zoo rollapplyr
#' @importFrom zoo as.zoo
#' @importFrom RcppRoll roll_sum
#'
#' @examples
#' dates <- as.Date("2019-01-01")+0:729
......@@ -48,57 +47,55 @@
#'## 2 0.4112304 24 1
#'## 3 0.1650978 48 1
#'## 4 0.2356849 48 1
IDF.agg <- function(data,ds,na.accept = 0,
which.stations = NULL,which.mon = 0:11,names = c('date','RR'),cl = NULL){
if(!inherits(data, "list"))stop("Argument 'data' must be a list, instead it is a: ", class(data))
# function 2: aggregate station data over durations and find annual maxima:
agg.station <- function(station){
data.s <- data[[station]]
if(!is.data.frame(data.s)){stop("Elements of 'data' must be data.frames. But element "
,station," contains: ", class(data.s))}
if(sum(is.element(names[1:2],names(data.s)))!=2){stop('Dataframe of station ', station
,' does not contain $', names[1]
,' or $', names[2], '.')}
dtime<-as.numeric((data.s[,names[1]][2]-data.s[,names[1]][1]),units="hours")
if(any(ds %% dtime != 0)){
stop('Given aggragation durations are not multiples of the time resolution = '
,dtime,' of station ',station,'.')}
# function 1: aggregate over single durations and find annual maxima:
agg.ts <- function(ds){
runmean <- rollapplyr(as.zoo(data.s[,names[2]]),ds/dtime,FUN=sum,fill =NA,align='right')
runmean <- runmean/ds #intensity per hour
subset <- is.element(as.POSIXlt(data.s[,names[1]])$mon,which.mon)
max <- tapply(runmean[subset],(as.POSIXlt(data.s[,names[1]])$year+1900)[subset],
function(vec){
n.na <- sum(is.na(vec))
max <- ifelse(n.na <= na.accept,max(vec,na.rm = TRUE),NA)
return(max)})
return(max) # maxima for single durations
IDF.agg <- function(data,ds,na.accept = 0,
which.stations = NULL,which.mon = 0:11,names = c('date','RR'),cl = NULL){
if(!inherits(data, "list"))stop("Argument 'data' must be a list, instead it is a: ", class(data))
# function 2: aggregate station data over durations and find annual maxima:
agg.station <- function(station){
data.s <- data[[station]]
if(!is.data.frame(data.s)){stop("Elements of 'data' must be data.frames. But element "
,station," contains: ", class(data.s))}
if(sum(is.element(names[1:2],names(data.s)))!=2){stop('Dataframe of station ', station
,' does not contain $', names[1]
,' or $', names[2], '.')}
dtime<-as.numeric((data.s[,names[1]][2]-data.s[,names[1]][1]),units="hours")
if(any(ds %% dtime > 10e-16)){
stop('At least one of the given aggregation durations is not multiple of the time resolution = '
,dtime,' of station ',station,'.')}
# function 1: aggregate over single durations and find annual maxima:
agg.ts <- function(ds){
runsum = RcppRoll::roll_sum(data.s[,names[2]],ds/dtime,fill=NA)
#runmean <- rollapplyr(as.zoo(data.s[,names[2]]),ds/dtime,FUN=sum,fill =NA,align='right')
runsum <- runsum/ds #intensity per hour
subset <- is.element(as.POSIXlt(data.s[,names[1]])$mon,which.mon)
max <- tapply(runsum[subset],(as.POSIXlt(data.s[,names[1]])$year+1900)[subset],
function(vec){
n.na <- sum(is.na(vec))
max <- ifelse(n.na <= na.accept,max(vec,na.rm = TRUE),NA)
return(max)})
df <- data.frame(xdat=max,ds=ds,year=as.numeric(names(max)))
return(df) # maxima for single durations
}
# call function 1 in lapply to aggregate over all durations at single station
data.agg <- pbapply::pblapply(ds,agg.ts,cl=cl)
#browser()
df <- do.call(rbind,data.agg)
return(df) # maxima for all durations at one station
}
# which stations should be used?
if(is.null(which.stations))which.stations <- 1:length(data)
# call function 2 in lapply to aggregate over all durations at all stations
station.list <- lapply(which.stations,agg.station)
return(do.call('rbind',station.list))
}
# call function 1 in lapply to aggregate over all durations at single station
data.agg <- pbsapply(ds,agg.ts,simplify=TRUE,cl=cl)
df <- data.frame(xdat = as.vector(data.agg), ds = rep(ds,each=length(data.agg[,1])))
df$station <- station
df$year <- rep(unique(as.POSIXlt(data.s[,names[1]])$year+1900),length(ds))
return(df) # maxima for all durations at one station
}
# which stations should be used?
if(is.null(which.stations))which.stations <- 1:length(data)
# call function 2 in lapply to aggregate over all durations at all stations
station.list <- lapply(which.stations,agg.station)
return(do.call('rbind',station.list))
}
#### IDF.plot ####
#' Plotting of IDF curves at a chosen station
#'
#' @param durations vector of durations for which to calculate the quantiles.
......@@ -108,8 +105,8 @@ IDF.agg <- function(data,ds,na.accept = 0,
#' (or \code{\link{gev.d.params}} for model with covariates).
#' @param probs vector of exeedance probabilities for which to plot IDF curves (p = 1-1/ReturnPeriod)
#' @param cols vector of colors for IDF curves. Should have same length as \code{probs}
#' @param add logical indicating if plot should be added to existing plot
#' @param legend logical indicating if legend should be plotted
#' @param add logical indicating if plot should be added to existing plot, default is FALSE
#' @param legend logical indicating if legend should be plotted (TRUE, the default)
#' @param ... additional parameters passed on to the \code{plot} function
#'
#' @export
......@@ -130,7 +127,7 @@ IDF.plot <- function(durations,fitparams,probs=c(0.5,0.9,0.99),
cols=4:2,add=FALSE,
legend=TRUE,...){
# if cols is to short, make longer
# if cols is to short, make longer
if(length(cols)!=length(probs))cols <- rep_len(cols,length.out=length(probs))
## calculate IDF values for given probability and durations
......
......@@ -8,7 +8,7 @@
#' @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 d positive numeric value, giving duration
#' @param ... additional parameters passed to \code{\link[evd]{dgev}}
#'
#' @details The duration dependent GEV distribution is defined after
......@@ -44,6 +44,7 @@ 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 (d<=0) {stop('The duration d has to be positive.')}
if(any(d+theta<=0)){
warning('Some shape parameters are negative,resulting from a negativ theta '
,theta, ' this will prododuce NAs.')}
......@@ -62,7 +63,7 @@ dgev.d <- function(q,mut,sigma0,xi,theta,eta,d,...) {
#' @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 d positive numeric value, giving duration
#' @param ... additional parameters passed to \code{\link[evd]{pgev}}
#'
#' @details The duration dependent GEV distribution is defined after
......@@ -87,6 +88,7 @@ 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 (d<=0) {stop('The duration d has to be positive.')}
if(any(d+theta<=0)){
warning('Some shape parameters are negative,resulting from a negativ theta '
,theta, ' this will prododuce NAs.')}
......@@ -105,7 +107,7 @@ pgev.d <- function(q,mut,sigma0,xi,theta,eta,d,...) {
#' @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 d positive numeric value, giving duration
#' @param ... additional parameters passed to \code{\link[evd]{qgev}}
#'
#' @details The duration dependent GEV distribution is defined after
......@@ -143,6 +145,7 @@ 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 (d<=0) {stop('The duration d has to be positive.')}
if(any(d+theta<=0)){
warning('Some shape parameters are negative,resulting from a negativ theta '
,theta, ' this will prododuce NAs.')}
......@@ -162,7 +165,7 @@ qgev.d <- function(p,mut,sigma0,xi,theta,eta,d,...) {
#' @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 d positive numeric value, giving duration
#'
#' @details The duration dependent GEV distribution is defined after
#' [Koutsoyannis et al., 1998]:
......@@ -194,8 +197,9 @@ 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 (d<=0) {stop('The duration d has to be positive.')}
if(any(d+theta<=0)){
warning('Some shape parameters are negative,resulting from a negativ theta '
warning('Some shape parameters are negative, resulting from a negative theta '
,theta, ' this will prododuce NAs.')}
# if denominator is negative NAs will be returned
if(d+theta<=0){return(rep(NA,n))}else{
......
......@@ -9,11 +9,12 @@
#' @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
#' value distribution, following Koutsoyiannis et al. (1998), 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 ds A vector of aggregation levels corresponding to the maxima in xdat.
#' 1/60 corresponds to 1 minute, 1 corresponds to 1 hour.
#' @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.
......@@ -23,33 +24,30 @@
#' Parameters are: modified location, scale_0, shape, duration offset, duration exponent repectively.
#' @param mulink,siglink,shlink,thetalink,etalink Link functions for generalized linear
#' modelling of the parameters, created with \code{\link{make.link}}.
#' @param muinit,siginit,shinit,thetainit,etainit Initial values as numeric of length
#' equal to total number of parameters. Alternatively initial values for only the parameter intercepts
#' can be passed to \code{init.vals}.
#' @param init.vals vector of length 5, giving initial values for parameter intercepts
#' 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
#' used to model the parameters (order: mu, sigma, xi, theta, eta). If rep(NA,5) (the default) is given, initial parameters are obtained
#' internally by fitting the GEV separately for each duration and applying a linear model to obtain the
#' duration dependency of the location and shape parameter.
#' @param theta_zero Logical value, indicating if theta should be estimated (FALSE, the default) or
#' should stay zero.
#' @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 theta_zero Logical value, indicating if theta parameter should be estimated (TRUE, the default) or
#' remain zero.
#' @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,
#' If \code{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{trans}{A 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{data}{data is standardized to standard Gumbel.}
#' \item{cov}{The covariance matrix.}
#' \item{vals}{Parameter values for every data point.}
#' \item{init.vals}{Initial values that where used.}
......@@ -77,10 +75,12 @@ gev.d.fit<-
function(xdat, ds, ydat = NULL, mul = NULL, sigl = NULL, shl = NULL, thetal = NULL, etal = NULL,
mulink = make.link("identity"), siglink = make.link("identity"), shlink = make.link("identity"),
thetalink = make.link("identity"), etalink = make.link("identity"),
muinit = NULL, siginit = NULL, shinit = NULL, thetainit = NULL, etainit = NULL,
show = TRUE, method = "Nelder-Mead", maxit = 10000, init.vals = NULL, theta_zero = FALSE, ...)
init.vals = rep(NA,5), theta_zero = FALSE,
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),'.'))
}
z <- list()
# number of parameters (betas) to estimate for each parameter:
npmu <- length(mul) + 1
......@@ -97,79 +97,76 @@ gev.d.fit<-
if(any(is.na(xdat))) stop('xdat contains NA values. NA values need to be removed first.')
# test if covariates matrix is given correctly
npar <- max(sapply(z$model,function(x){return(ifelse(is.null(x),0,max(x)))}))
if(any(npar>ncol(ydat),npar>0 & is.null(ydat)))stop("Not enough columns in covariates-Matrix 'ydat'.")
if(any(npar>ncol(ydat),npar>0 & is.null(ydat)))stop("Not enough columns in covariates matrix 'ydat'.")
# if no initial values where passed, calculate initial values for mu.d, sigma_0, xi, eta using IDF.init
if(!is.null(init.vals)){
if(length(init.vals)!=5){
warning('Parameter init.vals is not used, because it is not of length 5.')
init.vals <- NULL
}else{
init.vals <- data.frame(mu = init.vals[1], sigma = init.vals[2], xi = init.vals[3]
# initial values
if(length(init.vals)!=5) {
warning('Parameter init.vals is not used, because it is not of length 5.')
init.vals <- rep(NA,5)
}
if(!any(is.na(init.vals))){ #all initial values are given
init.vals <- data.frame(mu = init.vals[1], sigma = init.vals[2], xi = init.vals[3]
,theta = init.vals[4], eta = init.vals[5])
}else if(any(!is.na(init.vals))) { #some initial values are given
init.vals.user <- init.vals
init.vals <- gev.d.init(xdat,ds,z$link) #calculate init.vals using gev.d.init
for (i in 1:length(init.vals)){ #overwrite the calculated initial values with the ones given by the user
if(!is.na(init.vals.user[i])) {
init.vals[i]<-init.vals.user[i]
}
}
if(any(is.null(c(muinit,siginit,shinit,etainit)))& is.null(init.vals)){
# message('Initial values are calculated.')
}else{ #no initial values are given
init.vals <- gev.d.init(xdat,ds,z$link)
}
if(theta_zero==TRUE) { #if theta should stay zero
init.vals[4] = 0
}
# generate covariates matrices:
if (is.null(mul)) {
if (is.null(mul)) { #stationary
mumat <- as.matrix(rep(1, length(xdat)))
if (is.null(muinit))
muinit <- init.vals$mu
}else {
muinit <- init.vals$mu
}else { #non stationary
z$trans <- TRUE
mumat <- cbind(rep(1, length(xdat)), ydat[, mul])
if (is.null(muinit))
muinit <- c(init.vals$mu, rep(0, length(mul)))
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
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)))
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
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)))
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 <- init.vals$theta
thetainit <- init.vals$theta
}else {
z$trans <- TRUE
thmat <- cbind(rep(1, length(xdat)), ydat[, thetal])
if (is.null(thetainit))
thetainit <- c(init.vals$theta, rep(0, length(thetal)))
thetainit <- c(init.vals$theta, rep(0, length(thetal)))
}
if (is.null(etal)) {
etmat <- as.matrix(rep(1, length(xdat)))
if (is.null(etainit))
etainit <- init.vals$eta
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)))
etainit <- c(init.vals$eta, rep(0, length(etal)))
}
if(!theta_zero){#When theta parameter is included (default)
if(!theta_zero){#When theta parameter is not included (default)
init <- c(muinit, siginit, shinit, thetainit, etainit)
}else{ #Do not return initial value for theta if user does not want theta, as Hessian will fail.
}else{ #Do not return initial value for theta, if user does not want theta, as Hessian will fail.
init <- c(muinit, siginit, shinit, etainit)
}
......@@ -179,7 +176,7 @@ gev.d.fit<-
mu <- mulink$linkinv(mumat %*% (a[1:npmu]))
sigma <- siglink$linkinv(sigmat %*% (a[seq(npmu + 1, length = npsc)]))
xi <- shlink$linkinv(shmat %*% (a[seq(npmu + npsc + 1, length = npsh)]))
#Next line will set the theta likelihood as non-existent in case user requested it.
# Next line will set the theta likelihood as non-existent in case user requested it.
if(!theta_zero) {theta <- thetalink$linkinv(thmat %*% (a[seq(npmu + npsc + npsh + 1, length = npth)]))}
eta <- etalink$linkinv(etmat %*% (a[seq(npmu + npsc + npsh + npth + 1, length = npet)]))
......@@ -188,9 +185,9 @@ gev.d.fit<-
y <- xdat/sigma.d - mu
y <- 1 + xi * y
if(!theta_zero){ #When user wants theta parameter (default)
if(!theta_zero){ #When user does not want to estimate theta parameter (default)
if(any(eta <= 0) || any(theta < 0) || any(sigma.d <= 0) || any(y <= 0)) return(10^6)
}else{ #When user did not ask for theta parameter
}else{
if(any(eta <= 0) || any(sigma.d <= 0) || any(y <= 0)) return(10^6)
}
......@@ -240,14 +237,18 @@ gev.d.fit<-
z$ds <- ds
z$theta_zero <- theta_zero #Indicates if theta parameter was set to zero by user.
if(show) {
if(z$trans) # for nonstationary fit
if(z$trans) { # for nonstationary fit
print(z[c(2, 4)]) # print model, link (3) , conv
#TODO: print link function names
else print(z[4]) # for stationary fit print only conv
if(!z$conv) # if fit converged
# print names of link functions:
cat('$link\n')
print(c(z$link$mulink$name,z$link$siglink$name,z$link$shlink$name,z$link$thetalink$name,z$link$etalink$name))
cat('\n')
}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"
class(z) <- "gev.d.fit"
invisible(z)
}
......@@ -256,10 +257,10 @@ gev.d.fit<-
# function to get initial values for gev.d.fit:
# obtain initial values
# by fitting every duration seperately
# by fitting every duration separately
# possible ways to improve:
# take given initial values into accout, if there are any
# take given initial values into account, 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
......@@ -268,7 +269,6 @@ 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
#' @param link list of 5, link functions for parameters, created with \code{\link{make.link}}
#' @return list of initail values for mu_tilde, sigma_0, xi, eta
#' @importFrom stats lm
......@@ -291,7 +291,7 @@ gev.d.init <- function(xdat,ds,link){
# sig0 <- exp Intercept
siginit <- link$siglink$linkfun(exp(lmsig$coefficients[[1]]))
# eta <- mean of negativ slopes
etainit <- link$etalink$linkfun(median(c(-lmsig$coefficients[[2]],-lmmu$coefficients[[2]])))
etainit <- link$etalink$linkfun(mean(c(-lmsig$coefficients[[2]],-lmmu$coefficients[[2]])))
# mean of mu_d/sig_d
# could try:
# mu0/sig0 = exp(lmmu$coefficients[[1]])/exp(lmsig$coefficients[[1]])
......@@ -307,7 +307,7 @@ gev.d.init <- function(xdat,ds,link){
#' computes negative log-likelihood of d-gev model
#'
#' @param xdat numeric vector containing observations
#' @param ds numeric vector containing coresponding durations
#' @param ds numeric vector containing coresponding durations (1/60 corresponds to 1 minute, 1 corresponds to 1 hour)
#' @param mut,sig0,xi,theta,eta numeric vectors containing corresponding mles for each of the parameters
#'
#' @return single value containing negative log likelihood
......
......@@ -4,12 +4,19 @@
\alias{IDF.agg}
\title{Aggregation and annual maxima for choosen durations}
\usage{
IDF.agg(data, ds, na.accept = 0, which.stations = NULL,
which.mon = 0:11, names = c("date", "RR"), cl = NULL)
IDF.agg(
data,
ds,
na.accept = 0,
which.stations = NULL,
which.mon = 0:11,
names = c("date", "RR"),
cl = NULL
)
}
\arguments{
\item{data}{list of data.frames containing time series for every station.
The data.frame must have the columns 'date' and 'RR' unless onther names
The data.frame must have the columns 'date' and 'RR' unless other names
are specified in the parameter `names`. The column 'date' must contain strings with
standard date format.}
......@@ -38,8 +45,8 @@ the function \code{\link{gev.d.fit}}.
}
\details{
If data contains stations with different time resolutions that need to be aggregated at
different durations, IDF.agg needs to be run seperately for the different groups of stations.
Afterwards he results can be joint together using `rbind`.
different durations, IDF.agg needs to be run separately for the different groups of stations.
Afterwards the results can be joint together using `rbind`.
}
\examples{
dates <- as.Date("2019-01-01")+0:729
......
......@@ -4,8 +4,15 @@
\alias{IDF.plot}
\title{Plotting of IDF curves at a chosen station}
\usage{
IDF.plot(durations, fitparams, probs = c(0.5, 0.9, 0.99), cols = 4:2,
add = FALSE, legend = TRUE, ...)
IDF.plot(
durations,
fitparams,
probs = c(0.5, 0.9, 0.99),
cols = 4:2,
add = FALSE,
legend = TRUE,
...
)
}
\arguments{
\item{durations}{vector of durations for which to calculate the quantiles.}
......@@ -19,9 +26,9 @@ as obtained from \code{\link{gev.d.fit}}
\item{cols}{vector of colors for IDF curves. Should have same length as \code{probs}}
\item{add}{logical indicating if plot should be added to existing plot}
\item{add}{logical indicating if plot should be added to existing plot, default is FALSE}
\item{legend}{logical indicating if legend should be plotted}
\item{legend}{logical indicating if legend should be plotted (TRUE, the default)}
\item{...}{additional parameters passed on to the \code{plot} function}
}
......
......@@ -15,7 +15,7 @@ dgev.d(q, mut, sigma0, xi, theta, eta, d, ...)
\item{eta}{numeric value, giving duration exponent (defining slope of the IDF curve)}
\item{d}{numeric value, giving duration}
\item{d}{positive numeric value, giving duration}
\item{...}{additional parameters passed to \code{\link[evd]{dgev}}}
}
......
......@@ -4,10 +4,19 @@
\alias{gev.d.diag}
\title{Diagnostic Plots for d-gev Models}
\usage{
gev.d.diag(fit, subset = NULL, cols = NULL, pch = NULL,
which = "both", mfrow = c(1, 2), legend = TRUE,
gev.d.diag(
fit,
subset = NULL,
cols = NULL,
pch = NULL,
which = "both",
mfrow = c(1, 2),
legend = TRUE,
title = c("Residual Probability Plot", "Residual Quantile Plot"),
emp.lab = "Empirical", mod.lab = "Model", ...)
emp.lab = "Empirical",
mod.lab = "Model",
...
)