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

Merge branch 'covariates' of gitlab.met.fu-berlin.de:Rpackages/IDF into covariates

parents 7b78ee3a e34403ab
...@@ -7,59 +7,59 @@ ...@@ -7,59 +7,59 @@
#### gev.d.fit #### #### gev.d.fit ####
#' @title Maximum-likelihood Fitting of the duration-dependent GEV Distribution #' @title Maximum-likelihood Fitting of the duration-dependent GEV Distribution
#' @description Modified \code{\link[ismev]{gev.fit}} function for Maximum-likelihood fitting #' @description Modified \code{\link[ismev]{gev.fit}} function for Maximum-likelihood fitting
#' for the duration-dependent generalized extreme #' for the duration-dependent generalized extreme
#' value distribution, following Koutsoyiannis et al. (1998), including generalized linear #' value distribution, following Koutsoyiannis et al. (1998), including generalized linear
#' modeling of each parameter. #' modeling of each parameter.
#' @param xdat A vector containing maxima for different durations. #' @param xdat A vector containing maxima for different durations.
#' This can be obtained from \code{\link{IDF.agg}}. #' 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. #' 1/60 corresponds to 1 minute, 1 corresponds to 1 hour.
#' @param ydat A matrix of covariates for generalized linear modeling of the parameters #' @param ydat A matrix of covariates for generalized linear modeling of the parameters
#' (or NULL (the default) for stationary fitting). The number of rows should be the same as the #' (or NULL (the default) for stationary fitting). The number of rows should be the same as the
#' length of xdat. #' length of xdat.
#' @param mutl,sigma0l,xil,thetal,etal Numeric vectors of integers, giving the columns of ydat that contain #' @param mutl,sigma0l,xil,thetal,etal Numeric vectors of integers, giving the columns of ydat that contain
#' covariates for generalized linear modeling of the parameters (or NULL (the default) #' covariates for generalized linear modeling of the parameters (or NULL (the default)
#' if the corresponding parameter is stationary). #' if the corresponding parameter is stationary).
#' Parameters are: modified location, scale offset, shape, duration offset, duration exponent, respectively. #' Parameters are: modified location, scale offset, shape, duration offset, duration exponent, respectively.
#' @param mutlink,sigma0link,xilink,thetalink,etalink Link functions for generalized linear #' @param mutlink,sigma0link,xilink,thetalink,etalink Link functions for generalized linear
#' modeling of the parameters, created with \code{\link{make.link}}. The default is \code{make.link("identity")}. #' modeling of the parameters, created with \code{\link{make.link}}. The default is \code{make.link("identity")}.
#' @param init.vals list of length 5, giving initial values for all or some parameters #' @param init.vals list of length 5, giving initial values for all or some parameters
#' (order: mut, sigma0, xi, theta, eta). If as.list(rep(NA,5)) (the default) is given, initial parameters are obtained #' (order: mut, sigma0, xi, theta, eta). If as.list(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 #' 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. #' duration dependency of the location and shape parameter.
#' Initial values for covariate parameters are assumed as 0 if not given. #' Initial values for covariate parameters are assumed as 0 if not given.
#' @param theta_zero Logical value, indicating if theta should be estimated (FALSE, the default) or #' @param theta_zero Logical value, indicating if theta should be estimated (FALSE, the default) or
#' should stay zero. #' should stay zero.
#' @param show Logical; if TRUE (the default), print details of the fit. #' @param show Logical; if TRUE (the default), print details of the fit.
#' @param method The optimization method used in \code{\link{optim}}. #' @param method The optimization method used in \code{\link{optim}}.
#' @param maxit The maximum number of iterations. #' @param maxit The maximum number of iterations.
#' @param ... Other control parameters for the optimization. #' @param ... Other control parameters for the optimization.
#' @return A list containing the following components. #' @return A list containing the following components.
#' A subset of these components are printed after the fit. #' A subset of these components are printed after the fit.
#' If \code{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. #' the components nllh, mle and se are always printed.
#' \item{nllh}{single numeric giving the negative log-likelihood value} #' \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, #' \item{mle}{numeric vector giving the MLE's for the modified location, scale_0, shape,
#' duration offset and duration exponent, resp.} #' duration offset and duration exponent, resp.}
#' \item{se}{numeric vector giving the standard errors for the MLE's (in the same order)} #' \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{trans}{A logical indicator for a non-stationary fit.}
#' \item{model}{A list with components mutl, sigma0l, xil, thetal and etal.} #' \item{model}{A list with components mutl, sigma0l, xil, thetal and etal.}
#' \item{link}{A character vector giving inverse link functions.} #' \item{link}{A character vector giving inverse link functions.}
#' \item{conv}{The convergence code, taken from the list returned by \code{\link{optim}}. #' \item{conv}{The convergence code, taken from the list returned by \code{\link{optim}}.
#' A zero indicates successful convergence.} #' A zero indicates successful convergence.}
#' \item{data}{data is standardized to standard Gumbel.} #' \item{data}{data is standardized to standard Gumbel.}
#' \item{cov}{The covariance matrix.} #' \item{cov}{The covariance matrix.}
#' \item{vals}{Parameter values for every data point.} #' \item{vals}{Parameter values for every data point.}
#' \item{init.vals}{Initial values that were used.} #' \item{init.vals}{Initial values that were used.}
#' \item{ds}{Durations for every data point.} #' \item{ds}{Durations for every data point.}
#' @details For details on the d-GEV and the parameter definitions, see \link{IDF-package}. #' @details For details on the d-GEV and the parameter definitions, see \link{IDF-package}.
#' @seealso \code{\link{IDF-package}}, \code{\link{IDF.agg}}, \code{\link{gev.fit}}, \code{\link{optim}} #' @seealso \code{\link{IDF-package}}, \code{\link{IDF.agg}}, \code{\link{gev.fit}}, \code{\link{optim}}
#' @export #' @export
#' @importFrom stats optim #' @importFrom stats optim
#' @importFrom stats make.link #' @importFrom stats make.link
#' #'
#' @examples #' @examples
#' # sampled random data from d-gev with covariates #' # sampled random data from d-gev with covariates
#' # GEV parameters: #' # GEV parameters:
#' # mut = 4 + 0.2*cov1 +0.5*cov2 #' # mut = 4 + 0.2*cov1 +0.5*cov2
...@@ -67,24 +67,24 @@ ...@@ -67,24 +67,24 @@
#' # xi = 0.5 #' # xi = 0.5
#' # theta = 0 #' # theta = 0
#' # eta = 0.5 #' # eta = 0.5
#' #'
#' data('example',package ='IDF') #' data('example',package ='IDF')
#' #'
#' gev.d.fit(xdat=example$dat,ds = example$d,ydat=as.matrix(example[,c('cov1','cov2')]) #' gev.d.fit(xdat=example$dat,ds = example$d,ydat=as.matrix(example[,c('cov1','cov2')])
#' ,mutl=c(1,2),sigma0l=1) #' ,mutl=c(1,2),sigma0l=1)
gev.d.fit<- gev.d.fit<-
function(xdat, ds, ydat = NULL, mutl = NULL, sigma0l = NULL, xil = NULL, thetal = NULL, etal = NULL, function(xdat, ds, ydat = NULL, mutl = NULL, sigma0l = NULL, xil = NULL, thetal = NULL, etal = NULL,
mutlink = make.link("identity"), sigma0link = make.link("identity"), xilink = make.link("identity"), mutlink = make.link("identity"), sigma0link = make.link("identity"), xilink = make.link("identity"),
thetalink = make.link("identity"), etalink = make.link("identity"), thetalink = make.link("identity"), etalink = make.link("identity"),
init.vals = as.list(rep(NA,5)), theta_zero = FALSE, init.vals = as.list(rep(NA,5)), theta_zero = FALSE,
show = TRUE, method = "Nelder-Mead", maxit = 10000, ...) show = TRUE, method = "Nelder-Mead", maxit = 10000, ...)
{ {
if (length(xdat) != length(ds)) { if (length(xdat) != length(ds)) {
stop(paste0('The length of xdat is ',length(xdat),', but the length of ds is ',length(ds),'.')) stop(paste0('The length of xdat is ',length(xdat),', but the length of ds is ',length(ds),'.'))
} }
z <- list() z <- list()
# number of parameters (betas) to estimate for each parameter: # number of parameters (betas) to estimate for each parameter:
npmu <- length(mutl) + 1 npmu <- length(mutl) + 1
npsc <- length(sigma0l) + 1 npsc <- length(sigma0l) + 1
npsh <- length(xil) + 1 npsh <- length(xil) + 1
...@@ -93,16 +93,16 @@ gev.d.fit<- ...@@ -93,16 +93,16 @@ gev.d.fit<-
z$trans <- FALSE # indicates if fit is non-stationary z$trans <- FALSE # indicates if fit is non-stationary
z$model <- list(mutl, sigma0l, xil, thetal, etal) z$model <- list(mutl, sigma0l, xil, thetal, etal)
z$link <- list(mutlink=mutlink, sigma0link=sigma0link, xilink=xilink, thetalink=thetalink, etalink=etalink) z$link <- list(mutlink=mutlink, sigma0link=sigma0link, xilink=xilink, thetalink=thetalink, etalink=etalink)
# test for NA values: # test for NA values:
if(any(is.na(xdat))) stop('xdat contains NA values. NA values need to be removed first.') if(any(is.na(xdat))) stop('xdat contains NA values. NA values need to be removed first.')
# test for finite values: # test for finite values:
if(any(is.infinite(xdat))) stop('xdat contains non finite values. Inf and -Inf need to be removed first.') if(any(is.infinite(xdat))) stop('xdat contains non finite values. Inf and -Inf need to be removed first.')
# test if covariates matrix is given correctly # test if covariates matrix is given correctly
npar <- max(sapply(z$model,function(x){return(ifelse(is.null(x),0,max(x)))})) 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'.")
# initial values # initial values
if(length(init.vals)!=5 | !is.list(init.vals)) { if(length(init.vals)!=5 | !is.list(init.vals)) {
warning('Parameter init.vals is not used, because it is no list of length 5.') warning('Parameter init.vals is not used, because it is no list of length 5.')
...@@ -118,11 +118,11 @@ gev.d.fit<- ...@@ -118,11 +118,11 @@ gev.d.fit<-
init.vals[[i]]<-init.vals.user[[i]] init.vals[[i]]<-init.vals.user[[i]]
} }
} }
}else{ #no initial values are given }else{ #no initial values are given
init.vals <- gev.d.init(xdat,ds,z$link) init.vals <- gev.d.init(xdat,ds,z$link)
} }
# generate covariates matrices: # generate covariates matrices:
if (is.null(mutl)) { #stationary if (is.null(mutl)) { #stationary
mumat <- as.matrix(rep(1, length(xdat))) mumat <- as.matrix(rep(1, length(xdat)))
muinit <- init.vals$mu muinit <- init.vals$mu
...@@ -141,7 +141,7 @@ gev.d.fit<- ...@@ -141,7 +141,7 @@ gev.d.fit<-
} }
if (is.null(xil)) { if (is.null(xil)) {
shmat <- as.matrix(rep(1, length(xdat))) shmat <- as.matrix(rep(1, length(xdat)))
shinit <- init.vals$xi shinit <- init.vals$xi
}else { }else {
z$trans <- TRUE z$trans <- TRUE
shmat <- cbind(rep(1, length(xdat)), ydat[, xil]) shmat <- cbind(rep(1, length(xdat)), ydat[, xil])
...@@ -163,7 +163,7 @@ gev.d.fit<- ...@@ -163,7 +163,7 @@ gev.d.fit<-
etmat <- cbind(rep(1, length(xdat)), ydat[, etal]) etmat <- cbind(rep(1, length(xdat)), ydat[, etal])
etainit <- c(init.vals$eta, rep(0, length(etal)))[1:npet] etainit <- c(init.vals$eta, rep(0, length(etal)))[1:npet]
} }
if(!theta_zero){#When theta parameter is not included (default) if(!theta_zero){#When theta parameter is not included (default)
init <- c(muinit, siginit, shinit, thetainit, etainit) init <- c(muinit, siginit, shinit, thetainit, etainit)
...@@ -178,29 +178,29 @@ gev.d.fit<- ...@@ -178,29 +178,29 @@ gev.d.fit<-
mu <- mutlink$linkinv(mumat %*% (a[1:npmu])) mu <- mutlink$linkinv(mumat %*% (a[1:npmu]))
sigma <- sigma0link$linkinv(sigmat %*% (a[seq(npmu + 1, length = npsc)])) sigma <- sigma0link$linkinv(sigmat %*% (a[seq(npmu + 1, length = npsc)]))
xi <- xilink$linkinv(shmat %*% (a[seq(npmu + npsc + 1, length = npsh)])) 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. # 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)]))} 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)])) eta <- etalink$linkinv(etmat %*% (a[seq(npmu + npsc + npsh + npth + 1, length = npet)]))
ifelse(!theta_zero, ds.t <- ds+theta, ds.t <- ds) #Don't use theta if user requested not to have it. ifelse(!theta_zero, ds.t <- ds+theta, ds.t <- ds) #Don't use theta if user requested not to have it.
sigma.d <- sigma/(ds.t^eta) sigma.d <- sigma/(ds.t^eta)
y <- xdat/sigma.d - mu y <- xdat/sigma.d - mu
y <- 1 + xi * y y <- 1 + xi * y
if(!theta_zero){ #When user wants to estimate theta parameter (default) if(!theta_zero){ #When user wants to estimate theta parameter (default)
if(any(eta <= 0) || any(theta < 0) || any(sigma.d <= 0) || any(y <= 0)) return(10^6) if(any(eta <= 0) || any(theta < 0) || any(sigma.d <= 0) || any(y <= 0)) return(10^6)
}else{ }else{
if(any(eta <= 0) || any(sigma.d <= 0) || any(y <= 0)) return(10^6) if(any(eta <= 0) || any(sigma.d <= 0) || any(y <= 0)) return(10^6)
} }
sum(log(sigma.d)) + sum(y^(-1/xi)) + sum(log(y) * (1/xi + 1)) sum(log(sigma.d)) + sum(y^(-1/xi)) + sum(log(y) * (1/xi + 1))
} }
# finding minimum of log-likelihood: # finding minimum of log-likelihood:
x <- optim(init, gev.lik, hessian = TRUE, method = method, x <- optim(init, gev.lik, hessian = TRUE, method = method,
control = list(maxit = maxit, ...)) control = list(maxit = maxit, ...))
# saving output parameters: # saving output parameters:
z$conv <- x$convergence z$conv <- x$convergence
mut <- mutlink$linkinv(mumat %*% (x$par[1:npmu])) mut <- mutlink$linkinv(mumat %*% (x$par[1:npmu]))
...@@ -215,16 +215,16 @@ gev.d.fit<- ...@@ -215,16 +215,16 @@ gev.d.fit<-
z$nllh <- x$value z$nllh <- x$value
# normalize data to standard Gumbel: # normalize data to standard Gumbel:
sc.d <- sc0/((ds+theta)^eta) sc.d <- sc0/((ds+theta)^eta)
z$data <- - log(as.vector((1 + xi * (xdat/sc.d-mut))^(-1/xi))) z$data <- - log(as.vector((1 + xi * (xdat/sc.d-mut))^(-1/xi)))
z$mle <- x$par z$mle <- x$par
test <- try( # catch error test <- try( # catch error
z$cov <- solve(x$hessian) # invert hessian to get estimation on var-covar-matrix z$cov <- solve(x$hessian) # invert hessian to get estimation on var-covar-matrix
,silent = TRUE ) ,silent = TRUE )
if("try-error" %in% class(test)){ if("try-error" %in% class(test)){
warning("Hessian could not be inverted. NAs were produced.") warning("Hessian could not be inverted. NAs were produced.")
z$cov <- matrix(NA,length(z$mle),length(z$mle)) z$cov <- matrix(NA,length(z$mle),length(z$mle))
} }
z$se <- sqrt(diag(z$cov)) # sqrt(digonal entries) = standart error of mle's z$se <- sqrt(diag(z$cov)) # sqrt(digonal entries) = standart error of mle's
if (!theta_zero) {#When theta parameter is returned (default) if (!theta_zero) {#When theta parameter is returned (default)
z$vals <- cbind(mut, sc0, xi, theta, eta) z$vals <- cbind(mut, sc0, xi, theta, eta)
} else {#When theta parameter is not returned, asked by user } else {#When theta parameter is not returned, asked by user
...@@ -237,16 +237,16 @@ gev.d.fit<- ...@@ -237,16 +237,16 @@ gev.d.fit<-
colnames(z$vals) <-c('mut','sigma0','xi','eta') colnames(z$vals) <-c('mut','sigma0','xi','eta')
} }
z$ds <- ds z$ds <- ds
z$theta_zero <- theta_zero #Indicates if theta parameter was set to zero by user. z$theta_zero <- theta_zero #Indicates if theta parameter was set to zero by user.
if(show) { if(show) {
if(z$trans) { # for nonstationary fit if(z$trans) { # for nonstationary fit
print(z[c(2, 4)]) # print model, link (3) , conv print(z[c(2, 4)]) # print model, link (3) , conv
# print names of link functions: # print names of link functions:
cat('$link\n') cat('$link\n')
print(c(z$link$mutlink$name,z$link$sigma0link$name,z$link$xilink$name,z$link$thetalink$name,z$link$etalink$name)) print(c(z$link$mutlink$name,z$link$sigma0link$name,z$link$xilink$name,z$link$thetalink$name,z$link$etalink$name))
cat('\n') cat('\n')
}else{print(z[4])} # for stationary fit print only conv }else{print(z[4])} # for stationary fit print only conv
if(!z$conv){ # if fit converged if(!z$conv){ # if fit converged
print(z[c(5, 7, 9)]) # print nll, mle, se print(z[c(5, 7, 9)]) # print nll, mle, se
} }
} }
...@@ -258,7 +258,7 @@ gev.d.fit<- ...@@ -258,7 +258,7 @@ gev.d.fit<-
#### gev.d.init #### #### gev.d.init ####
# function to get initial values for gev.d.fit: # function to get initial values for gev.d.fit:
# obtain initial values # obtain initial values
# by fitting every duration separately # by fitting every duration separately
# possible ways to improve: # possible ways to improve:
...@@ -273,10 +273,10 @@ gev.d.fit<- ...@@ -273,10 +273,10 @@ gev.d.fit<-
#' @param ds vector of durations belonging to maxima in xdat #' @param ds vector of durations belonging to maxima in xdat
#' @param link list of 5, link functions for parameters, created with \code{\link{make.link}} #' @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 #' @return list of initail values for mu_tilde, sigma_0, xi, eta
#' @importFrom stats lm #' @importFrom stats lm
#' @importFrom stats median #' @importFrom stats median
#' @importFrom ismev gev.fit #' @importFrom ismev gev.fit
#' @keywords internal #' @keywords internal
gev.d.init <- function(xdat,ds,link){ gev.d.init <- function(xdat,ds,link){
durs <- unique(ds) durs <- unique(ds)
...@@ -289,19 +289,19 @@ gev.d.init <- function(xdat,ds,link){ ...@@ -289,19 +289,19 @@ gev.d.init <- function(xdat,ds,link){
# get values for sig0 and eta (also mu_0) from linear model in log-log scale # get values for sig0 and eta (also mu_0) from linear model in log-log scale
lmsig <- lm(log(mles[,2])~log(durs)) lmsig <- lm(log(mles[,2])~log(durs))
lmmu <- lm(log(mles[,1])~log(durs)) lmmu <- lm(log(mles[,1])~log(durs))
# sig0 <- exp Intercept # sig0 <- exp Intercept
siginit <- link$sigma0link$linkfun(exp(lmsig$coefficients[[1]])) siginit <- link$sigma0link$linkfun(exp(lmsig$coefficients[[1]]))
# eta <- mean of negativ slopes # eta <- mean of negativ slopes
etainit <- link$etalink$linkfun(mean(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 # mean of mu_d/sig_d
# could try: # could try:
# mu0/sig0 = exp(lmmu$coefficients[[1]])/exp(lmsig$coefficients[[1]]) # mu0/sig0 = exp(lmmu$coefficients[[1]])/exp(lmsig$coefficients[[1]])
muinit <- link$mutlink$linkfun(median(c(mles[,1]/mles[,2]),na.rm = TRUE)) muinit <- link$mutlink$linkfun(median(c(mles[,1]/mles[,2]),na.rm = TRUE))
# mean of shape parameters # mean of shape parameters
shinit <- link$xilink$linkfun(median(mles[,3],na.rm = TRUE)) shinit <- link$xilink$linkfun(median(mles[,3],na.rm = TRUE))
thetainit <- link$thetalink$linkfun(0) thetainit <- link$thetalink$linkfun(0)
return(list(mu=muinit,sigma=siginit,xi=shinit,theta=thetainit,eta=etainit)) return(list(mu=muinit,sigma=siginit,xi=shinit,theta=thetainit,eta=etainit))
} }
...@@ -315,7 +315,7 @@ gev.d.init <- function(xdat,ds,link){ ...@@ -315,7 +315,7 @@ gev.d.init <- function(xdat,ds,link){
#' @param mut,sigma0,xi,theta,eta numeric vectors containing corresponding estimates for each of the parameters #' @param mut,sigma0,xi,theta,eta numeric vectors containing corresponding estimates for each of the parameters
#' @param log Logical; if TRUE, the log likelihood is returned. #' @param log Logical; if TRUE, the log likelihood is returned.
#' #'
#' @return single value containing (log) likelihood #' @return single value containing (log) likelihood
#' @export #' @export
#' #'
#' @examples #' @examples
...@@ -329,35 +329,35 @@ gev.d.init <- function(xdat,ds,link){ ...@@ -329,35 +329,35 @@ gev.d.init <- function(xdat,ds,link){
#' ,theta = params[,4],eta = params[,5],log=TRUE) #' ,theta = params[,4],eta = params[,5],log=TRUE)
gev.d.lik <- function(xdat,ds,mut,sigma0,xi,theta,eta,log=FALSE) { gev.d.lik <- function(xdat,ds,mut,sigma0,xi,theta,eta,log=FALSE) {
if(any(xi==0)){stop('Function is not defined for shape parameter of zero.')} if(any(xi==0)){stop('Function is not defined for shape parameter of zero.')}
if(any(! c(length(ds),length(mut),length(sigma0),length(xi),length(theta),length(eta)) %in% if(any(! c(length(ds),length(mut),length(sigma0),length(xi),length(theta),length(eta)) %in%
c(1,length(xdat)))){ c(1,length(xdat)))){
stop('Input vectors differ in length, but must have the same length.') stop('Input vectors differ in length, but must have the same length.')
} }
ds.t <- ds+theta ds.t <- ds+theta
sigma.d <- sigma0/(ds.t^eta) sigma.d <- sigma0/(ds.t^eta)
y <- xdat/sigma.d - mut y <- xdat/sigma.d - mut
y <- 1 + xi * y y <- 1 + xi * y
if(log){ if(log){
return(sum(log(sigma.d) + y^(-1/xi) + log(y) * (1/xi + 1))) return(sum(log(sigma.d) + y^(-1/xi) + log(y) * (1/xi + 1)))
}else{ }else{
return(prod(sigma.d * exp(y^(-1/xi)) * y ^ (1/xi + 1))) return(prod(sigma.d * exp(y^(-1/xi)) * y ^ (1/xi + 1)))
} }
} }
#### gev.d.diag #### #### gev.d.diag ####
#' Diagnostic Plots for d-gev Models #' Diagnostic Plots for d-gev Models
#' #'
#' @description Produces diagnostic plots for d-gev models using #' @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 #' the output of the function \code{\link{gev.d.fit}}. Values for different durations can be plotted in
#' different colors of with different symbols. #' different colors of with different symbols.
#' @param fit object returned by \code{\link{gev.d.fit}} #' @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 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(fit$ds)} to #' @param cols optional either one value or vector of same length as \code{unique(fit$ds)} to
#' specify the colors of plotting points. #' specify the colors of plotting points.
#' The default uses the \code{rainbow} function. #' The default uses the \code{rainbow} function.
#' @param pch optional either one value or vector of same length as \code{unique(fit$ds)} containing #' @param pch optional either one value or vector of same length as \code{unique(fit$ds)} containing
#' integers or symbols to specify the plotting points. #' integers or symbols to specify the plotting points.
...@@ -367,7 +367,7 @@ gev.d.lik <- function(xdat,ds,mut,sigma0,xi,theta,eta,log=FALSE) { ...@@ -367,7 +367,7 @@ gev.d.lik <- function(xdat,ds,mut,sigma0,xi,theta,eta,log=FALSE) {
#' @param legend logical indicating if legends should be plotted #' @param legend logical indicating if legends should be plotted
#' @param title character vector of length 2, giving the titles for the pp- and the qq-plot #' @param title character vector of length 2, giving the titles for the pp- and the qq-plot
#' @param emp.lab,mod.lab character string containing names for empirical and model axis #' @param emp.lab,mod.lab character string containing names for empirical and model axis
#' @param ... additional parameters passed on to the plotting function #' @param ... additional parameters passed on to the plotting function
#' #'
#' @export #' @export
#' @importFrom graphics plot abline par title #' @importFrom graphics plot abline par title
...@@ -375,12 +375,12 @@ gev.d.lik <- function(xdat,ds,mut,sigma0,xi,theta,eta,log=FALSE) { ...@@ -375,12 +375,12 @@ gev.d.lik <- function(xdat,ds,mut,sigma0,xi,theta,eta,log=FALSE) {
#' #'
#' @examples #' @examples
#' data('example',package ='IDF') #' data('example',package ='IDF')
#' #'
#' fit <- gev.d.fit(xdat=example$dat,ds = example$d,ydat=as.matrix(example[,c('cov1','cov2')]) #' fit <- gev.d.fit(xdat=example$dat,ds = example$d,ydat=as.matrix(example[,c('cov1','cov2')])
#' ,mutl=c(1,2),sigma0l=1) #' ,mutl=c(1,2),sigma0l=1)
#' # diagnostic plots for complete data #' # diagnostic plots for complete data
#' gev.d.diag(fit,pch=1) #' gev.d.diag(fit,pch=1)
#' # diagnostic plots for subset of data (e.g. one station) #' # diagnostic plots for subset of data (e.g. one station)
#' gev.d.diag(fit,subset = example$cov1==1,pch=1) #' gev.d.diag(fit,subset = example$cov1==1,pch=1)
gev.d.diag <- function(fit,subset=NULL,cols=NULL,pch=NULL,which='both',mfrow=c(1,2),legend=TRUE, gev.d.diag <- function(fit,subset=NULL,cols=NULL,pch=NULL,which='both',mfrow=c(1,2),legend=TRUE,
title=c('Residual Probability Plot','Residual Quantile Plot'), title=c('Residual Probability Plot','Residual Quantile Plot'),
...@@ -391,7 +391,7 @@ gev.d.diag <- function(fit,subset=NULL,cols=NULL,pch=NULL,which='both',mfrow=c(1 ...@@ -391,7 +391,7 @@ gev.d.diag <- function(fit,subset=NULL,cols=NULL,pch=NULL,which='both',mfrow=c(1
# subset data # subset data
df <- data.frame(data=fit$data,ds=fit$ds) df <- data.frame(data=fit$data,ds=fit$ds)
if(!is.null(subset)){ if(!is.null(subset)){
if(dim(df)[1]!=length(subset)){stop("Length of 'subset' does not match length of data if(dim(df)[1]!=length(subset)){stop("Length of 'subset' does not match length of data
'xdat' used for fitting.")} 'xdat' used for fitting.")}
df <- df[subset,] df <- df[subset,]
} }
...@@ -400,9 +400,9 @@ gev.d.diag <- function(fit,subset=NULL,cols=NULL,pch=NULL,which='both',mfrow=c(1 ...@@ -400,9 +400,9 @@ gev.d.diag <- function(fit,subset=NULL,cols=NULL,pch=NULL,which='both',mfrow=c(1
# rescale durations to assign colors # rescale durations to assign colors
df$cval <- sapply(df$ds,function(d){which(durs==d)}) df$cval <- sapply(df$ds,function(d){which(durs==d)})
# sort data # sort data
df <- df[order(df$data),] df <- df[order(df$data),]
# plotting position # plotting position
n <- length(df$data) n <- length(df$data)
px <- (1:n)/(n + 1) px <- (1:n)/(n + 1)
...@@ -412,7 +412,7 @@ gev.d.diag <- function(fit,subset=NULL,cols=NULL,pch=NULL,which='both',mfrow=c(1 ...@@ -412,7 +412,7 @@ gev.d.diag <- function(fit,subset=NULL,cols=NULL,pch=NULL,which='both',mfrow=c(1
# colors and symbols # colors and symbols
if(is.null(cols))cols <- rainbow(length(durs)) if(is.null(cols))cols <- rainbow(length(durs))
if(is.null(pch))pch <- df$cval if(is.null(pch))pch <- df$cval
if(which=='both'|which=='pp'){ if(which=='both'|which=='pp'){
# pp # pp
plot(px, exp( - exp( - df$data)), xlab = plot(px, exp( - exp( - df$data)), xlab =
...@@ -438,15 +438,15 @@ gev.d.diag <- function(fit,subset=NULL,cols=NULL,pch=NULL,which='both',mfrow=c(1 ...@@ -438,15 +438,15 @@ gev.d.diag <- function(fit,subset=NULL,cols=NULL,pch=NULL,which='both',mfrow=c(1
#' Calculate gev(d) parameters from \code{gev.d.fit} output #' Calculate gev(d) parameters from \code{gev.d.fit} output
#' #'
#' @description function to calculate mut, sigma0, xi, theta, eta #' @description function to calculate mut, sigma0, xi, theta, eta
#' (modified location, scale offset, shape, duration offset, duration exponent) #' (modified location, scale offset, shape, duration offset, duration exponent)
#' from results of \code{\link{gev.d.fit}} with covariates or link funktions other than identity. #' from results of \code{\link{gev.d.fit}} with covariates or link funktions other than identity.
#' @param fit fit object returned by \code{\link{gev.d.fit}} or \code{\link{gev.fit}} #' @param fit fit object returned by \code{\link{gev.d.fit}} or \code{\link{gev.fit}}
#' @param ydat A matrix containing the covariates in the same order as used in \code{gev.d.fit}. #' @param ydat A matrix containing the covariates in the same order as used in \code{gev.d.fit}.
#' @seealso \code{\link{IDF-package}} #' @seealso \code{\link{IDF-package}}
#' @return data.frame containing mu_tilde, sigma0, xi, theta, eta (or mu, sigma, xi for gev.fit objects) #' @return data.frame containing mu_tilde, sigma0, xi, theta, eta (or mu, sigma, xi for gev.fit objects)
#' @export #' @export
#' #'
#' @examples #' @examples
#' data('example',package = 'IDF') #' data('example',package = 'IDF')
#' fit <- gev.d.fit(example$dat,example$d,ydat = as.matrix(example[,c("cov1","cov2")]) #' fit <- gev.d.fit(example$dat,example$d,ydat = as.matrix(example[,c("cov1","cov2")])
...@@ -457,7 +457,7 @@ gev.d.diag <- function(fit,subset=NULL,cols=NULL,pch=NULL,which='both',mfrow=c(1 ...@@ -457,7 +457,7 @@ gev.d.diag <- function(fit,subset=NULL,cols=NULL,pch=NULL,which='both',mfrow=c(1
gev.d.params <- function(fit,ydat=NULL){ gev.d.params <- function(fit,ydat=NULL){
if(!class(fit)%in%c("gev.d.fit","gev.fit"))stop("'fit' must be an object returned by 'gev.d.fit' or 'gev.fit'.") if(!class(fit)%in%c("gev.d.fit","gev.fit"))stop("'fit' must be an object returned by 'gev.d.fit' or 'gev.fit'.")