Commit 109388cb authored by Felix Fauer's avatar Felix Fauer
Browse files

merge parametrization (with multiscaling and flattening) into master, which...

merge parametrization (with multiscaling and flattening) into master, which had minor bugfixes in the meantime (aggregation function and likelihood sign)

Merge branch 'parametrization' into development
parents 1871bd41 43f461b4
......@@ -27,6 +27,11 @@
#' and shape parameter \eqn{\xi\in R}, \eqn{\xi\neq 0}.
#' The parameters \eqn{\theta \leq 0} and \eqn{0<\eta<1} are duration offset and duration exponent
#' and describe the slope and curvature in the resulting IDF curves, respectively.
#' * 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{\mu(x)=\tilde{\mu}(\sigma_0(d+\theta)^{-\eta_1}+\tau)}
#' * A useful introduction to __Maximum Likelihood Estimation__ for fitting for the
#' generalized extreme value distribution (GEV) is provided by Coles (2001). It should be noted, however, that this method uses
#' the assumption that block maxima (of different durations or stations) are independent of each other.
......@@ -225,13 +230,11 @@ NULL
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 too short, make longer
if(length(cols)!=length(probs))cols <- rep_len(cols,length.out=length(probs))
## calculate IDF values for given probability and durations
qs <- lapply(durations,qgev.d,p=probs,mut=fitparams[1],sigma0=fitparams[2],xi=fitparams[3],
theta=fitparams[4],eta=fitparams[5])
theta=fitparams[4],eta=fitparams[5], eta2=fitparams[6], tau=fitparams[7])
idf.array <- matrix(unlist(qs),length(probs),length(durations)) # array[probs,durs]
if(!add){ #new plot
## initialize plot window
......@@ -249,7 +252,6 @@ IDF.plot <- function(durations,fitparams,probs=c(0.5,0.9,0.99),
# empty plot
plot(NA,xlim=xlim,ylim=ylim,xlab="Duration [h]",ylab="Intensity [mm/h]",log="xy",main=main)
}
## plot IDF curves
for(i in 1:length(probs)){
lines(durations,idf.array[i,],col=cols[i],...)
......
......@@ -10,6 +10,8 @@
#' shape parameter \eqn{\xi}.
#' @param theta numeric value, giving duration offset \eqn{\theta} (defining curvature of the IDF curve)
#' @param eta numeric value, giving duration exponent \eqn{\eta} (defining slope of the IDF curve)
#' @param eta2 numeric value, giving a second duration exponent \eqn{\eta_2} (needed for multiscaling). Default: NULL, treated as \eqn{\eta_2=\eta}.
#' @param tau numeric value, giving intensity offset \eqn{\tau} (defining flattening of the IDF curve). Default: \eqn{\tau=0}.
#' @param d positive numeric value, giving duration
#' @param ... additional parameters passed to \code{\link[evd]{dgev}}
#'
......@@ -37,9 +39,10 @@
#' 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. ',
dgev.d <- function(q,mut,sigma0,xi,theta,eta,d,eta2=NULL,tau=0,...) {
if(is.null(eta2)){eta2=eta}
if(any(c(length(mut),length(sigma0),length(xi),length(theta),length(eta),length(eta2),length(tau))>1)){
message('One of the parameters mut, sigma0, xi, theta, eta, tau 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)){
......@@ -47,8 +50,12 @@ dgev.d <- function(q,mut,sigma0,xi,theta,eta,d,...) {
,theta, ' this will prododuce NAs.')}
# if denominator is negative NAs will be returned
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,...))}
#sigma.d <-sigma0/(d+theta)^eta+ tau # old
sigma.d <- sigma0/(d+theta)^eta2 +tau
mu.d <- mut*(sigma0/(d+theta)^eta +tau)
return(dgev(q,loc=mu.d,scale=sigma.d,shape=xi,...))}
}
......@@ -61,6 +68,8 @@ 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 eta2 numeric value, giving a second duration exponent \eqn{\eta_2} (needed for multiscaling). Default: NULL, treated as \eqn{\eta_2=\eta}.
#' @param tau numeric value, giving intensity offset \eqn{\tau} (defining flattening of the IDF curve). Default: \eqn{\tau=0}.
#' @param d positive numeric value, giving duration
#' @param ... additional parameters passed to \code{\link[evd]{pgev}}
#'
......@@ -83,9 +92,10 @@ dgev.d <- function(q,mut,sigma0,xi,theta,eta,d,...) {
#' @examples
#' x <- seq(4,20,0.1)
#' prob <- pgev.d(q=x,mut=4,sigma0=2,xi=0,theta=0.1,eta=0.1,d=1)
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. ',
pgev.d <- function(q,mut,sigma0,xi,theta,eta,d,tau=0,eta2=NULL, ...) {
if(is.null(eta2)){eta2=eta}
if(any(c(length(mut),length(sigma0),length(xi),length(theta),length(eta),length(eta2),length(tau))>1)){
message('One of the parameters mut, sigma0, xi, theta, eta, tau 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)){
......@@ -93,8 +103,12 @@ pgev.d <- function(q,mut,sigma0,xi,theta,eta,d,...) {
,theta, ' this will prododuce NAs.')}
# if denominator is negative NAs will be returned
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,...))}
#sigma.d <-sigma0/(d+theta)^eta+tau # old
sigma.d <- sigma0/(d+theta)^eta2 +tau
mu.d <- mut*(sigma0/(d+theta)^eta +tau)
return(pgev(q,loc=mu.d,scale=sigma.d,shape=xi,...))}
}
......@@ -105,8 +119,10 @@ pgev.d <- function(q,mut,sigma0,xi,theta,eta,d,...) {
#' @description Quantile function of duration-dependent GEV distribution (inverse of the cumulative probability distribution function)
#' @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 theta numeric value, giving duration offset (defining curvature of the IDF curve for short durations)
#' @param eta numeric value, giving duration exponent (defining slope of the IDF curve)
#' @param eta2 numeric value, giving a second duration exponent \eqn{\eta_2} (needed for multiscaling). Default: NULL, treated as \eqn{\eta_2=\eta}.
#' @param tau numeric value, giving intensity offset \eqn{\tau} (defining flattening of the IDF curve). Default: \eqn{\tau=0}.
#' @param d positive numeric value, giving duration
#' @param ... additional parameters passed to \code{\link[evd]{qgev}}
#'
......@@ -129,7 +145,7 @@ pgev.d <- function(q,mut,sigma0,xi,theta,eta,d,...) {
#' @examples
#' p <- c(0.5,0.9,0.99)
#' # calulate quantiles for one duration
#' qgev.d(p=p,mut=4,sigma0=2,xi=0,theta=0.1,eta=0.3,d=1)
#' qgev.d(p=p,mut=4,sigma0=2,xi=0,theta=0.1,eta=0.3, d=1)
#'
#' # calculate quantiles for sequence of durations
#' ds <- 2^seq(0,4,0.1)
......@@ -142,9 +158,10 @@ pgev.d <- function(q,mut,sigma0,xi,theta,eta,d,...) {
#' }
#' legend('topright',title = 'p-quantile',
#' legend = 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. ',
qgev.d <- function(p,mut,sigma0,xi,theta,eta,d,tau=0,eta2=NULL, ...) {
if (is.null(eta2)){eta2=eta}
if(any(c(length(mut),length(sigma0),length(xi),length(theta),length(eta), length(eta2), length(tau))>1)){
message('One of the parameters mut, sigma0, xi, theta, eta, eta2, tau 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)){
......@@ -152,9 +169,17 @@ qgev.d <- function(p,mut,sigma0,xi,theta,eta,d,...) {
,theta, ' this will prododuce NAs.')}
# if denominator is negative NAs will be returned
if(d+theta<=0){return(rep(NA,length(p)))}else{
sigma.d <-sigma0/(d+theta)^eta
return(qgev(p,loc=as.numeric(mut*sigma.d)
#sigma.d <-sigma0/(d+theta)^eta
#sigma.d <-sigma0/(d+theta)^eta+tau
sigma.d <- sigma0/(d+theta)^eta2 +tau
mu.d <- mut*(sigma0/(d+theta)^eta +tau)
return(qgev(p,loc=as.numeric(mu.d)
,scale=as.numeric(sigma.d),shape=as.numeric(xi),...))}
#return(qgev(p,loc=as.numeric(mut*sigma.d) # old
# ,scale=as.numeric(sigma.d),shape=as.numeric(xi),...))} # old
}
......@@ -167,6 +192,8 @@ 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 eta2 numeric value, giving a second duration exponent \eqn{\eta_2} (needed for multiscaling). Default: NULL, treated as \eqn{\eta_2=\eta}.
#' @param tau numeric value, giving intensity offset \eqn{\tau} (defining flattening of the IDF curve). Default: \eqn{\tau=0}.
#' @param d positive numeric value, giving duration
#'
#' @details For details on the d-GEV and the parameter definitions, see \link{IDF-package}
......@@ -193,9 +220,10 @@ qgev.d <- function(p,mut,sigma0,xi,theta,eta,d,...) {
#' hist(samp[[2]],breaks = 10,add=TRUE,col=rgb(0,0,1,0.5),freq = FALSE)
#' legend('topright',fill = c(rgb(1,0,0,0.5),rgb(0,0,1,0.5)),
#' legend = paste('d=',1:2,'h'),title = 'Duration')
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. ',
rgev.d <- function(n,mut,sigma0,xi,theta,eta,d,tau=0,eta2=NULL) {
if (is.null(eta2)){eta2=eta}
if(any(c(length(mut),length(sigma0),length(xi),length(theta),length(eta),length(eta2),length(tau))>1)){
message('One of the parameters mut, sigma0, xi, theta, eta, tau 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)){
......@@ -203,8 +231,12 @@ rgev.d <- function(n,mut,sigma0,xi,theta,eta,d) {
,theta, ' this will prododuce NAs.')}
# if denominator is negative NAs will be returned
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))}
#sigma.d <-sigma0/(d+theta)^eta+tau # old
sigma.d <- sigma0/(d+theta)^eta2 +tau
mu.d <- mut*(sigma0/(d+theta)^eta +tau)
return(rgev(n,loc=mu.d,scale=sigma.d,shape=xi))}
}
......@@ -18,19 +18,22 @@
#' @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
#' 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,taul,eta2l Numeric vectors of integers, giving the columns of ydat that contain
#' covariates for generalized linear modeling of the parameters (or NULL (the default)
#' if the corresponding parameter is stationary).
#' 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,eta2link,taulink Link functions for generalized linear
#' 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
#' (order: mut, sigma0, xi, theta, eta). If as.list(rep(NA,5)) (the default) is given, initial parameters are obtained
#' @param init.vals list, giving initial values for all or some parameters
#' (order: mut, sigma0, xi, theta, eta, eta2, tau). If one of those parameters shall not be used (see theta_zero, eta2_zero, tau_zero),
#' the number of parameters has to be reduced accordingly. If some or all given values in init.vals are NA or
#' no init.vals at all is declared (the default), 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.
#' 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 whether theta should be estimated (FALSE, the default) or
#' should stay zero.
#' @param tau_zero,eta2_zero Logical values, indicating whether tau,eta2 should be estimated (TRUE, the default).
#' @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.
......@@ -41,10 +44,10 @@
#' 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.}
#' duration offset and duration exponent, resp. If requested, contains also second duration exponent and intensity-offset}
#' \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.}
#' \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{conv}{The convergence code, taken from the list returned by \code{\link{optim}}.
#' A zero indicates successful convergence.}
......@@ -67,6 +70,8 @@
#' # xi = 0.5
#' # theta = 0
#' # eta = 0.5
#' # eta2 = 0.5
#' # tau = 0
#'
#' data('example',package ='IDF')
#'
......@@ -74,10 +79,10 @@
#' ,mutl=c(1,2),sigma0l=1)
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, taul = NULL, eta2l = NULL,
mutlink = make.link("identity"), sigma0link = make.link("identity"), xilink = make.link("identity"),
thetalink = make.link("identity"), etalink = make.link("identity"),
init.vals = as.list(rep(NA,5)), theta_zero = FALSE,
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, ...)
{
if (length(xdat) != length(ds)) {
......@@ -90,9 +95,11 @@ gev.d.fit<-
npsh <- length(xil) + 1
npth <- ifelse(!theta_zero,length(thetal) + 1,0)
npet <- length(etal) + 1
npta <- ifelse(!tau_zero, length(taul) + 1,0)
npe2 <- ifelse(!eta2_zero, length(eta2l) + 1,0)
z$trans <- FALSE # indicates if fit is non-stationary
z$model <- list(mutl, sigma0l, xil, thetal, etal)
z$link <- list(mutlink=mutlink, sigma0link=sigma0link, xilink=xilink, thetalink=thetalink, etalink=etalink)
z$model <- list(mutl, sigma0l, xil, thetal, etal, eta2l, taul)
z$link <- list(mutlink=mutlink, sigma0link=sigma0link, xilink=xilink, thetalink=thetalink, etalink=etalink, eta2link=eta2link, taulink=taulink)
# test for NA values:
if(any(is.na(xdat))) stop('xdat contains NA values. NA values need to be removed first.')
......@@ -104,23 +111,47 @@ gev.d.fit<-
if(any(npar>ncol(ydat),npar>0 & is.null(ydat)))stop("Not enough columns in covariates matrix 'ydat'.")
# initial values
if(length(init.vals)!=5 | !is.list(init.vals)) {
warning('Parameter init.vals is not used, because it is no list of length 5.')
init.vals <- as.list(rep(NA,5))
}
if(!any(is.na(init.vals))){ #all initial values are given
names(init.vals) <- c('mu','sigma','xi','theta','eta')
}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]]
init.necessary.length = 4 + as.numeric(!theta_zero) + as.numeric(!eta2_zero) + as.numeric(!tau_zero) # mut, sigma0, xi, theta, eta, eta2, tau
if (is.null(init.vals)) {init.vals = as.list(rep(NA,init.necessary.length))
}else{init.vals = as.list(init.vals)}
if(length(init.vals)!=init.necessary.length | !is.list(init.vals)) {
print(paste0('Parameter init.vals is not used, because it is no list of length ',init.necessary.length,'.'))
init.vals <- gev.d.init(xdat,ds,z$link)
}else{ # length of given values is correct
# name given initial values
names1=c('mu','sigma','xi') # standard set of parameters
if (!theta_zero){names1=c(names1,'theta')} # add theta (in case)
names1=c(names1,'eta') # add eta (always)
if (!eta2_zero){names1=c(names1,'eta2')} # add eta2 (in case)
if (!tau_zero){names1=c(names1,'tau')} # add tau (in case)
names(init.vals) <- names1
# add missing initial value (fixed internal number of parameters: 7)
if (theta_zero) init.vals$theta = 0
if (eta2_zero) init.vals$eta2 = init.vals$eta
if (tau_zero) init.vals$tau = 0
init.vals=init.vals[c("mu","sigma","xi","theta","eta","eta2","tau")]
iv=init.vals
init.vals = list(mu=iv$mu,sigma=iv$sigma,xi=iv$xi,theta=iv$theta,eta=iv$eta,eta2=iv$eta2,tau=iv$tau)
if(!any(is.na(init.vals))){ #all initial values are given
# do nothing
}else if(any(!is.na(init.vals))) { #some initial values are given
#if (!eta2_zero) print("autmoatic inital value setting not implemented yet for multiscaling (eta2_zero=FALSE)")
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[[names(init.vals.user)[i]]])) {
init.vals[[names(init.vals.user)[i]]]<-init.vals.user[[names(init.vals.user)[i]]]
}
}
}else{ #no initial values are given
#if (!eta2_zero) print("autmoatic inital value setting not implemented yet for multiscaling (eta2_zero=FALSE)")
init.vals <- gev.d.init(xdat,ds,z$link)
}
}else{ #no initial values are given
init.vals <- gev.d.init(xdat,ds,z$link)
}
}
# generate covariates matrices:
if (is.null(mutl)) { #stationary
......@@ -163,37 +194,75 @@ gev.d.fit<-
etmat <- cbind(rep(1, length(xdat)), ydat[, etal])
etainit <- c(init.vals$eta, rep(0, length(etal)))[1:npet]
}
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.
init <- c(muinit, siginit, shinit, etainit)
if (is.null(taul)) {
tamat <- as.matrix(rep(1, length(xdat)))
tauinit <- init.vals$tau
}else {
z$trans <- TRUE
tamat <- cbind(rep(1, length(xdat)), ydat[, taul])
tauinit <- c(init.vals$tau, rep(0, length(taul)))[1:npta]
}
if (is.null(eta2l)) {
e2mat <- as.matrix(rep(1, length(xdat)))
eta2init <- init.vals$eta2
}else {
z$trans <- TRUE
e2mat <- cbind(rep(1, length(xdat)), ydat[, eta2l])
eta2init <- c(init.vals$eta2, rep(0, length(eta2l)))[1:npe2]
}
init <- c(muinit,siginit,shinit)
if (!theta_zero) {init <- c(init,thetainit)} # add theta init (in case)
init <- c(init,etainit) # add eta init (always)
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:
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.
# 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)]))}
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)]))}
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)
y <- xdat/sigma.d - mu
#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
y <- 1 + xi * y
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)
}else{
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))
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))
}
......@@ -212,10 +281,27 @@ gev.d.fit<-
theta <- thetalink$linkinv(thmat %*% (0))
}
eta <- etalink$linkinv(etmat %*% (x$par[seq(npmu + npsc + npsh + npth + 1, length = npet)]))
if(!eta2_zero){ #When user wants to use eta2 parameter
eta2 <- eta2link$linkinv(e2mat %*% (x$par[seq(npmu + npsc + npsh + npth + npet + 1,length = npe2)]))
}else{ #When user requests not to have eta2 parameter (default)
eta2 <- eta
}
if(!tau_zero){ #When user does NOT set tau parameter to zero (not default)
tau <- taulink$linkinv(tamat %*% (x$par[seq(npmu + npsc + npsh + npth + npet + npe2 + 1,length = npta)]))
}else{ #When user requests tau parameter to be zero
tau <- taulink$linkinv(tamat %*% (0))
}
z$nllh <- x$value
# normalize data to standard Gumbel:
sc.d <- sc0/((ds+theta)^eta)
z$data <- - log(as.vector((1 + xi * (xdat/sc.d-mut))^(-1/xi)))
#sc.d <- sc0/((ds+theta)^eta)+tau # old
sc.d <- sc0/((ds+theta)^eta2)+tau # new
mut.d <- mut*(sc0/((ds+theta)^eta )+tau) # new
#z$data <- - log(as.vector((1 + xi * (xdat/sc.d-mut))^(-1/xi))) # old
z$data <- - log(as.vector((1 + xi * ((xdat-mut.d)/sc.d))^(-1/xi))) # new
z$mle <- x$par
test <- try( # catch error
z$cov <- solve(x$hessian) # invert hessian to get estimation on var-covar-matrix
......@@ -225,25 +311,43 @@ 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)
z$vals <- cbind(mut, sc0, xi, theta, eta)
'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
z$vals <- cbind(mut, sc0, xi, eta)
}
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)
z$init.vals <- unlist(init.vals)
if(!theta_zero){ #When theta parameter is returned (default)
colnames(z$vals) <-c('mut','sigma0','xi','theta','eta')
} else { #When theta parameter is not returned, asked by user
colnames(z$vals) <-c('mut','sigma0','xi','eta')
}
'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)
z$eta2_zero <- eta2_zero #Indicates if eta2 parameter was set to zero by user. (default)
if(show) {
if(z$trans) { # for nonstationary fit
print(z[c(2, 4)]) # print model, link (3) , conv
# print names of link functions:
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))
#if(!tau_zero){
print(c(z$link$mutlink$name,z$link$sigma0link$name,z$link$xilink$name,z$link$thetalink$name,z$link$etalink$name,z$link$eta2link$name,z$link$taulink$name))
#} else{
# print(c(z$link$mutlink$name,z$link$sigma0link$name,z$link$xilink$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
......@@ -272,7 +376,7 @@ gev.d.fit<-
#' @param xdat vector of maxima for different durations
#' @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}}
#' @return list of initial values for mu_tilde, sigma_0, xi, eta
#' @return list of initial values for mu_tilde, sigma_0, xi, theta, eta, eta2, tau
#' @importFrom stats lm
#' @importFrom stats median
#' @importFrom ismev gev.fit
......@@ -294,6 +398,7 @@ gev.d.init <- function(xdat,ds,link){
siginit <- link$sigma0link$linkfun(exp(lmsig$coefficients[[1]]))
# eta <- mean of negativ slopes
etainit <- link$etalink$linkfun(mean(c(-lmsig$coefficients[[2]],-lmmu$coefficients[[2]])))
eta2init <- etainit + 0.0
# mean of mu_d/sig_d
# could try:
# mu0/sig0 = exp(lmmu$coefficients[[1]])/exp(lmsig$coefficients[[1]])
......@@ -301,8 +406,9 @@ gev.d.init <- function(xdat,ds,link){
# mean of shape parameters
shinit <- link$xilink$linkfun(median(mles[,3],na.rm = TRUE))
thetainit <- link$thetalink$linkfun(0)
tauinit <- link$taulink$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, eta2=eta2init, tau=tauinit))
}
#### gev.d.lik ####
......@@ -312,7 +418,7 @@ gev.d.init <- function(xdat,ds,link){
#' Computes (log-) likelihood of d-GEV model
#' @param xdat numeric vector containing observations
#' @param ds numeric vector containing corresponding durations (1/60 corresponds to 1 minute, 1 corresponds to 1 hour)
#' @param mut,sigma0,xi,theta,eta numeric vectors containing corresponding estimates for each of the parameters
#' @param mut,sigma0,xi,theta,eta,eta2,tau numeric vectors containing corresponding estimates for each of the parameters
#' @param log Logical; if TRUE, the log likelihood is returned.
#'
#' @return single value containing (log) likelihood
......@@ -327,18 +433,24 @@ gev.d.init <- function(xdat,ds,link){
#' params <- gev.d.params(fit,ydat = as.matrix(test.set[c('cov1','cov2')]))
#' gev.d.lik(xdat = test.set$dat,ds = test.set$d,mut = params[,1],sigma0 = params[,2],xi = params[,3]
#' ,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,tau=0,eta2=NULL) {
if (is.null(eta2)){eta2=eta}
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),length(eta2),length(tau)) %in%
c(1,length(xdat)))){
stop('Input vectors differ in length, but must have the same length.')
}
ds.t <- ds+theta
sigma.d <- sigma0/(ds.t^eta)
y <- xdat/sigma.d - mut
sigma.d <- sigma0/(ds.t^eta2)+tau
mu.d <- mut*(sigma0/(ds.t^eta)+tau)
y = (xdat - mu.d) / sigma.d # new
y <- 1 + xi * y
#sigma.d <- sigma0/(ds.t^eta) + tau # old
#y <- xdat/sigma.d - mut # old
#y <- 1 + xi * y # old
if(log){