Commit 9d84b444 authored by Felix Fauer's avatar Felix Fauer
Browse files

finished adding eta2 to all functions. some minor checks are still missing....

finished adding eta2 to all functions. some minor checks are still missing. also readme shows bad diagnosis plots (worse than with only tau)
parent 4c34d9cc
...@@ -2,7 +2,6 @@ ...@@ -2,7 +2,6 @@
export(IDF.agg) export(IDF.agg)
export(IDF.plot) export(IDF.plot)
export(IDF.plot.fit)
export(dgev.d) export(dgev.d)
export(gev.d.diag) export(gev.d.diag)
export(gev.d.fit) export(gev.d.fit)
......
...@@ -222,12 +222,11 @@ NULL ...@@ -222,12 +222,11 @@ NULL
IDF.plot <- function(durations,fitparams,probs=c(0.5,0.9,0.99), IDF.plot <- function(durations,fitparams,probs=c(0.5,0.9,0.99),
cols=4:2,add=FALSE, cols=4:2,add=FALSE,
legend=TRUE,...){ legend=TRUE,...){
# if cols is too short, make longer
# if cols is to short, make longer
if(length(cols)!=length(probs))cols <- rep_len(cols,length.out=length(probs)) if(length(cols)!=length(probs))cols <- rep_len(cols,length.out=length(probs))
## calculate IDF values for given probability and durations ## calculate IDF values for given probability and durations
qs <- lapply(durations,qgev.d,p=probs,mut=fitparams[1],sigma0=fitparams[2],xi=fitparams[3], qs <- lapply(durations,qgev.d,p=probs,mut=fitparams[1],sigma0=fitparams[2],xi=fitparams[3],
theta=fitparams[4],eta=fitparams[5], tau=fitparams[6]) theta=fitparams[4],eta=fitparams[5], eta2=fitparams[6], tau=fitparams[7])
idf.array <- matrix(unlist(qs),length(probs),length(durations)) # array[probs,durs] idf.array <- matrix(unlist(qs),length(probs),length(durations)) # array[probs,durs]
if(!add){ #new plot if(!add){ #new plot
## initialize plot window ## initialize plot window
......
...@@ -10,6 +10,7 @@ ...@@ -10,6 +10,7 @@
#' shape parameter \eqn{\xi}. #' shape parameter \eqn{\xi}.
#' @param theta numeric value, giving duration offset \eqn{\theta} (defining curvature of the IDF curve) #' @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 eta numeric value, giving duration exponent \eqn{\eta} (defining slope of the IDF curve)
#' @param eta2 numeric value, giving a second duration exponent (needed for multiscaling). If multiscaling is not requested, \eqn{eta2=eta} should be used.
#' @param tau numeric value, giving intensity offset \eqn{\tau} (defining flattening of the IDF curve) #' @param tau numeric value, giving intensity offset \eqn{\tau} (defining flattening of the IDF curve)
#' @param d positive numeric value, giving duration #' @param d positive numeric value, giving duration
#' @param ... additional parameters passed to \code{\link[evd]{dgev}} #' @param ... additional parameters passed to \code{\link[evd]{dgev}}
...@@ -38,8 +39,9 @@ ...@@ -38,8 +39,9 @@
#' lines(x,dens[[i]],lty=i) #' lines(x,dens[[i]],lty=i)
#' } #' }
#' legend('topright',title = 'Duration',legend = 1:4,lty=1:4) #' legend('topright',title = 'Duration',legend = 1:4,lty=1:4)
dgev.d <- function(q,mut,sigma0,xi,theta,eta,d,tau=0,...) { dgev.d <- function(q,mut,sigma0,xi,theta,eta,d,eta2=NULL,tau=0,...) {
if(any(c(length(mut),length(sigma0),length(xi),length(theta),length(eta), length(tau))>1)){ 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. ', message('One of the parameters mut, sigma0, xi, theta, eta, tau is a vector. ',
'This is not intended and might cause an error.')} 'This is not intended and might cause an error.')}
if (d<=0) {stop('The duration d has to be positive.')} if (d<=0) {stop('The duration d has to be positive.')}
...@@ -48,8 +50,12 @@ dgev.d <- function(q,mut,sigma0,xi,theta,eta,d,tau=0,...) { ...@@ -48,8 +50,12 @@ dgev.d <- function(q,mut,sigma0,xi,theta,eta,d,tau=0,...) {
,theta, ' this will prododuce NAs.')} ,theta, ' this will prododuce NAs.')}
# if denominator is negative NAs will be returned # if denominator is negative NAs will be returned
if(d+theta<=0){return(rep(NA,length(q)))}else{ if(d+theta<=0){return(rep(NA,length(q)))}else{
sigma.d <-sigma0/(d+theta)^eta+ tau
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,...))}
} }
...@@ -62,6 +68,7 @@ dgev.d <- function(q,mut,sigma0,xi,theta,eta,d,tau=0,...) { ...@@ -62,6 +68,7 @@ dgev.d <- function(q,mut,sigma0,xi,theta,eta,d,tau=0,...) {
#' @param mut,sigma0,xi numeric value, giving modified location, modified scale and shape parameter #' @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)
#' @param eta numeric value, giving duration exponent (defining slope 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 (needed for multiscaling). If multiscaling is not requested, \eqn{eta2=eta} should be used.
#' @param tau numeric value, giving intensity offset (defining flattening of the IDF curve) #' @param tau numeric value, giving intensity offset (defining flattening of the IDF curve)
#' @param d positive numeric value, giving duration #' @param d positive numeric value, giving duration
#' @param ... additional parameters passed to \code{\link[evd]{pgev}} #' @param ... additional parameters passed to \code{\link[evd]{pgev}}
...@@ -85,8 +92,9 @@ dgev.d <- function(q,mut,sigma0,xi,theta,eta,d,tau=0,...) { ...@@ -85,8 +92,9 @@ dgev.d <- function(q,mut,sigma0,xi,theta,eta,d,tau=0,...) {
#' @examples #' @examples
#' x <- seq(4,20,0.1) #' 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) #' 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,tau=0,d, ...) { pgev.d <- function(q,mut,sigma0,xi,theta,eta,d,tau=0,eta2=NULL, ...) {
if(any(c(length(mut),length(sigma0),length(xi),length(theta),length(eta),length(tau))>1)){ 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. ', message('One of the parameters mut, sigma0, xi, theta, eta, tau is a vector. ',
'This is not intended and might cause an error.')} 'This is not intended and might cause an error.')}
if (d<=0) {stop('The duration d has to be positive.')} if (d<=0) {stop('The duration d has to be positive.')}
...@@ -95,8 +103,12 @@ pgev.d <- function(q,mut,sigma0,xi,theta,eta,tau=0,d, ...) { ...@@ -95,8 +103,12 @@ pgev.d <- function(q,mut,sigma0,xi,theta,eta,tau=0,d, ...) {
,theta, ' this will prododuce NAs.')} ,theta, ' this will prododuce NAs.')}
# if denominator is negative NAs will be returned # if denominator is negative NAs will be returned
if(d+theta<=0){return(rep(NA,length(q)))}else{ if(d+theta<=0){return(rep(NA,length(q)))}else{
sigma.d <-sigma0/(d+theta)^eta+tau
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,...))}
} }
...@@ -109,6 +121,7 @@ pgev.d <- function(q,mut,sigma0,xi,theta,eta,tau=0,d, ...) { ...@@ -109,6 +121,7 @@ pgev.d <- function(q,mut,sigma0,xi,theta,eta,tau=0,d, ...) {
#' @param mut,sigma0,xi numeric value, giving modified location, modified scale and shape parameter #' @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 for short durations) #' @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 eta numeric value, giving duration exponent (defining slope of the IDF curve)
#' @param eta2 numeric value, giving a second duration exponent (needed for multiscaling). If multiscaling is not requested, \eqn{eta2=eta} should be used.
#' @param tau numeric value, giving intensity offset (defining flattening of the IDF curve for long durations) #' @param tau numeric value, giving intensity offset (defining flattening of the IDF curve for long durations)
#' @param d positive numeric value, giving duration #' @param d positive numeric value, giving duration
#' @param ... additional parameters passed to \code{\link[evd]{qgev}} #' @param ... additional parameters passed to \code{\link[evd]{qgev}}
...@@ -145,9 +158,10 @@ pgev.d <- function(q,mut,sigma0,xi,theta,eta,tau=0,d, ...) { ...@@ -145,9 +158,10 @@ pgev.d <- function(q,mut,sigma0,xi,theta,eta,tau=0,d, ...) {
#' } #' }
#' legend('topright',title = 'p-quantile', #' legend('topright',title = 'p-quantile',
#' legend = p,lty=1:3,bty = 'n') #' legend = p,lty=1:3,bty = 'n')
qgev.d <- function(p,mut,sigma0,xi,theta,eta,d,tau=0, ...) { qgev.d <- function(p,mut,sigma0,xi,theta,eta,d,tau=0,eta2=NULL, ...) {
if(any(c(length(mut),length(sigma0),length(xi),length(theta),length(eta), length(tau))>1)){ if (is.null(eta2)){eta2=eta}
message('One of the parameters mut, sigma0, xi, theta, eta is a vector. ', 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.')} 'This is not intended and might cause an error.')}
if (d<=0) {stop('The duration d has to be positive.')} if (d<=0) {stop('The duration d has to be positive.')}
if(any(d+theta<=0)){ if(any(d+theta<=0)){
...@@ -156,9 +170,16 @@ qgev.d <- function(p,mut,sigma0,xi,theta,eta,d,tau=0, ...) { ...@@ -156,9 +170,16 @@ qgev.d <- function(p,mut,sigma0,xi,theta,eta,d,tau=0, ...) {
# if denominator is negative NAs will be returned # if denominator is negative NAs will be returned
if(d+theta<=0){return(rep(NA,length(p)))}else{ if(d+theta<=0){return(rep(NA,length(p)))}else{
#sigma.d <-sigma0/(d+theta)^eta #sigma.d <-sigma0/(d+theta)^eta
sigma.d <-sigma0/(d+theta)^eta+tau #sigma.d <-sigma0/(d+theta)^eta+tau
return(qgev(p,loc=as.numeric(mut*sigma.d)
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),...))} ,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
} }
...@@ -171,6 +192,7 @@ qgev.d <- function(p,mut,sigma0,xi,theta,eta,d,tau=0, ...) { ...@@ -171,6 +192,7 @@ qgev.d <- function(p,mut,sigma0,xi,theta,eta,d,tau=0, ...) {
#' @param mut,sigma0,xi numeric value, giving modified location, modified scale and shape parameter #' @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)
#' @param eta numeric value, giving duration exponent (defining slope 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 (needed for multiscaling). If multiscaling is not requested, \eqn{eta2=eta} should be used.
#' @param tau numeric value, giving intensity offset (defining flattening of the IDF curve) #' @param tau numeric value, giving intensity offset (defining flattening of the IDF curve)
#' @param d positive numeric value, giving duration #' @param d positive numeric value, giving duration
#' #'
...@@ -198,8 +220,9 @@ qgev.d <- function(p,mut,sigma0,xi,theta,eta,d,tau=0, ...) { ...@@ -198,8 +220,9 @@ qgev.d <- function(p,mut,sigma0,xi,theta,eta,d,tau=0, ...) {
#' hist(samp[[2]],breaks = 10,add=TRUE,col=rgb(0,0,1,0.5),freq = FALSE) #' 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('topright',fill = c(rgb(1,0,0,0.5),rgb(0,0,1,0.5)),
#' legend = paste('d=',1:2,'h'),title = 'Duration') #' legend = paste('d=',1:2,'h'),title = 'Duration')
rgev.d <- function(n,mut,sigma0,xi,theta,eta,d,tau=0) { rgev.d <- function(n,mut,sigma0,xi,theta,eta,d,tau=0,eta2=NULL) {
if(any(c(length(mut),length(sigma0),length(xi),length(theta),length(eta),length(tau))>1)){ 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. ', message('One of the parameters mut, sigma0, xi, theta, eta, tau is a vector. ',
'This is not intended and might cause an error.')} 'This is not intended and might cause an error.')}
if (d<=0) {stop('The duration d has to be positive.')} if (d<=0) {stop('The duration d has to be positive.')}
...@@ -208,8 +231,12 @@ rgev.d <- function(n,mut,sigma0,xi,theta,eta,d,tau=0) { ...@@ -208,8 +231,12 @@ rgev.d <- function(n,mut,sigma0,xi,theta,eta,d,tau=0) {
,theta, ' this will prododuce NAs.')} ,theta, ' this will prododuce NAs.')}
# if denominator is negative NAs will be returned # if denominator is negative NAs will be returned
if(d+theta<=0){return(rep(NA,n))}else{ if(d+theta<=0){return(rep(NA,n))}else{
sigma.d <-sigma0/(d+theta)^eta+tau
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))}
} }
...@@ -22,17 +22,18 @@ ...@@ -22,17 +22,18 @@
#' 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,taulink,eta2link 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")}. #' 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, 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, 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 #' 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 whether 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. #' should stay zero.
#' @param tau_zero,eta2_zero Logical values, indicating whether tau,eta2 should be estimated (TRUE, the default) or #' @param tau_zero,eta2_zero Logical values, indicating whether tau,eta2 should be estimated (TRUE, the default).
#' 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.
...@@ -69,6 +70,8 @@ ...@@ -69,6 +70,8 @@
#' # xi = 0.5 #' # xi = 0.5
#' # theta = 0 #' # theta = 0
#' # eta = 0.5 #' # eta = 0.5
#' # eta2 = 0.5
#' # tau = 0
#' #'
#' data('example',package ='IDF') #' data('example',package ='IDF')
#' #'
...@@ -79,7 +82,7 @@ gev.d.fit<- ...@@ -79,7 +82,7 @@ gev.d.fit<-
function(xdat, ds, ydat = NULL, mutl = NULL, sigma0l = NULL, xil = NULL, thetal = NULL, etal = NULL, taul = NULL, eta2l = 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"), mutlink = make.link("identity"), sigma0link = make.link("identity"), xilink = make.link("identity"),
thetalink = make.link("identity"), etalink = make.link("identity"), taulink = make.link("identity"), eta2link = make.link("identity"), 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, eta_zero=TRUE init.vals = NULL, theta_zero = FALSE, tau_zero=TRUE, eta2_zero=TRUE,
show = TRUE, method = "Nelder-Mead", maxit = 10000, ...) show = TRUE, method = "Nelder-Mead", maxit = 10000, ...)
{ {
if (length(xdat) != length(ds)) { if (length(xdat) != length(ds)) {
...@@ -95,8 +98,8 @@ gev.d.fit<- ...@@ -95,8 +98,8 @@ gev.d.fit<-
npta <- ifelse(!tau_zero, length(taul) + 1,0) npta <- ifelse(!tau_zero, length(taul) + 1,0)
npe2 <- ifelse(!eta2_zero, length(eta2l) + 1,0) npe2 <- ifelse(!eta2_zero, length(eta2l) + 1,0)
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, taul) z$model <- list(mutl, sigma0l, xil, thetal, etal, eta2l, taul)
z$link <- list(mutlink=mutlink, sigma0link=sigma0link, xilink=xilink, thetalink=thetalink, etalink=etalink, taulink=taulink, eta2link=eta2link) z$link <- list(mutlink=mutlink, sigma0link=sigma0link, xilink=xilink, thetalink=thetalink, etalink=etalink, eta2link=eta2link, taulink=taulink)
# 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.')
...@@ -133,7 +136,7 @@ gev.d.fit<- ...@@ -133,7 +136,7 @@ gev.d.fit<-
if(!any(is.na(init.vals))){ #all initial values are given if(!any(is.na(init.vals))){ #all initial values are given
# do nothing # do nothing
}else if(any(!is.na(init.vals))) { #some initial values are given }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)") #if (!eta2_zero) print("autmoatic inital value setting not implemented yet for multiscaling (eta2_zero=FALSE)")
init.vals.user <- init.vals init.vals.user <- init.vals
init.vals <- gev.d.init(xdat,ds,z$link) #calculate init.vals using gev.d.init 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 for (i in 1:length(init.vals)){ #overwrite the calculated initial values with the ones given by the user
...@@ -142,7 +145,7 @@ gev.d.fit<- ...@@ -142,7 +145,7 @@ gev.d.fit<-
} }
} }
}else{ #no initial values are given }else{ #no initial values are given
if (!eta2_zero) print("autmoatic inital value setting not implemented yet for multiscaling (eta2_zero=FALSE)") #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) init.vals <- gev.d.init(xdat,ds,z$link)
} }
} }
...@@ -242,27 +245,21 @@ gev.d.fit<- ...@@ -242,27 +245,21 @@ gev.d.fit<-
mu.d <- mu*(sigma/(ds.t^eta)+tau) mu.d <- mu*(sigma/(ds.t^eta)+tau)
} }
} }
#sigma.d <- sigma/(ds.t^eta) #sigma.d <- sigma/(ds.t^eta) # old
y = (xdat - mu.d) / sigma.d # new y = (xdat - mu.d) / sigma.d # new
#y = (xdat - mu*sigma.d) / sigma.d # derivation #y = (xdat - mu*sigma.d) / sigma.d # derivation
#y <- xdat/sigma.d - mu # old #y <- xdat/sigma.d - mu # old
y <- 1 + xi * y 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)
#}
#if(!tau_zero) {if(any(tau < 0)) return(10^6)} # check condition for tau as well.
if(!theta_zero) {if(any(theta < 0)) {return(10^6)} } # check definition condition for theta 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(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(!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. 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)) # xxx continue here sum(log(sigma.d)) + sum(y^(-1/xi)) + sum(log(y) * (1/xi + 1))
} }
...@@ -281,17 +278,27 @@ gev.d.fit<- ...@@ -281,17 +278,27 @@ gev.d.fit<-
theta <- thetalink$linkinv(thmat %*% (0)) theta <- thetalink$linkinv(thmat %*% (0))
} }
eta <- etalink$linkinv(etmat %*% (x$par[seq(npmu + npsc + npsh + npth + 1, length = npet)])) 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) 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 + 1,length = npta)])) 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 }else{ #When user requests tau parameter to be zero
tau <- taulink$linkinv(tamat %*% (0)) tau <- taulink$linkinv(tamat %*% (0))
} }
z$nllh <- x$value z$nllh <- x$value
# normalize data to standard Gumbel: # normalize data to standard Gumbel:
sc.d <- sc0/((ds+theta)^eta)+tau #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))) #z$data <- - log(as.vector((1 + xi * (xdat/sc.d-mut))^(-1/xi))) # old
z$data <- - log(as.vector((1 + xi * ((xdat-mut)/sc.d))^(-1/xi))) # new
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
...@@ -314,7 +321,7 @@ gev.d.fit<- ...@@ -314,7 +321,7 @@ gev.d.fit<-
z$vals <- cbind(mut, sc0, xi, eta) z$vals <- cbind(mut, sc0, xi, eta)
} }
}' }'
z$vals <- cbind(mut, sc0, xi, theta, eta, tau) z$vals <- cbind(mut, sc0, xi, theta, eta, eta2, tau)
z$init.vals <- unlist(init.vals) z$init.vals <- unlist(init.vals)
'names2 = c("mut","sigma0","xi") # fixed part for set of names 'names2 = c("mut","sigma0","xi") # fixed part for set of names
...@@ -322,18 +329,19 @@ gev.d.fit<- ...@@ -322,18 +329,19 @@ gev.d.fit<-
names2 = c(names2, "eta") # add eta (always) names2 = c(names2, "eta") # add eta (always)
if(!tau_zero){names2=c(names2, "tau")} # add tau (in case) if(!tau_zero){names2=c(names2, "tau")} # add tau (in case)
colnames(z$vals) <- names2' colnames(z$vals) <- names2'
colnames(z$vals) <- c("mut","sigma0","xi","theta","eta","tau") colnames(z$vals) <- c("mut","sigma0","xi","theta","eta","eta2","tau")
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.
z$tau_zero <- tau_zero #Indicates if tau 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(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')
#if(!tau_zero){ #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$taulink$name)) 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{ #} else{
# 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))
#} #}
...@@ -365,7 +373,7 @@ gev.d.fit<- ...@@ -365,7 +373,7 @@ gev.d.fit<-
#' @param xdat vector of maxima for different durations #' @param xdat vector of maxima for different durations
#' @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 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 lm
#' @importFrom stats median #' @importFrom stats median
#' @importFrom ismev gev.fit #' @importFrom ismev gev.fit
...@@ -387,6 +395,7 @@ gev.d.init <- function(xdat,ds,link){ ...@@ -387,6 +395,7 @@ gev.d.init <- function(xdat,ds,link){
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]])))
eta2init <- etainit + 0.0
# 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]])
...@@ -396,7 +405,7 @@ gev.d.init <- function(xdat,ds,link){ ...@@ -396,7 +405,7 @@ gev.d.init <- function(xdat,ds,link){
thetainit <- link$thetalink$linkfun(0) thetainit <- link$thetalink$linkfun(0)
tauinit <- link$taulink$linkfun(0) tauinit <- link$taulink$linkfun(0)
return(list(mu=muinit,sigma=siginit,xi=shinit,theta=thetainit,eta=etainit, tau=tauinit)) return(list(mu=muinit,sigma=siginit,xi=shinit,theta=thetainit,eta=etainit, eta2=eta2init, tau=tauinit))
} }
#### gev.d.lik #### #### gev.d.lik ####
...@@ -406,7 +415,7 @@ gev.d.init <- function(xdat,ds,link){ ...@@ -406,7 +415,7 @@ gev.d.init <- function(xdat,ds,link){
#' Computes (log-) likelihood of d-GEV model #' Computes (log-) likelihood of d-GEV model
#' @param xdat numeric vector containing observations #' @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 ds numeric vector containing corresponding durations (1/60 corresponds to 1 minute, 1 corresponds to 1 hour)
#' @param mut,sigma0,xi,theta,eta,tau 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. #' @param log Logical; if TRUE, the log likelihood is returned.
#' #'
#' @return single value containing (log) likelihood #' @return single value containing (log) likelihood
...@@ -421,18 +430,24 @@ gev.d.init <- function(xdat,ds,link){ ...@@ -421,18 +430,24 @@ gev.d.init <- function(xdat,ds,link){
#' params <- gev.d.params(fit,ydat = as.matrix(test.set[c('cov1','cov2')])) #' 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] #' 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) #' ,theta = params[,4],eta = params[,5],log=TRUE)
gev.d.lik <- function(xdat,ds,mut,sigma0,xi,theta,eta,log=FALSE,tau=0) { 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(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),length(tau)) %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)))){ 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) + tau sigma.d <- sigma0/(ds.t^eta2)+tau
y <- xdat/sigma.d - mut mu.d <- mut*(sigma0/(ds.t^eta)+tau)
y = (xdat - mu.d) / sigma.d # new
y <- 1 + xi * y 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){ 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{
...@@ -538,7 +553,7 @@ gev.d.diag <- function(fit,subset=NULL,cols=NULL,pch=NULL,which='both',mfrow=c(1 ...@@ -538,7 +553,7 @@ gev.d.diag <- function(fit,subset=NULL,cols=NULL,pch=NULL,which='both',mfrow=c(1
#' @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, eta2, tau (or mu, sigma, xi for gev.fit objects)
#' @export #' @export
#' #'
#' @examples #' @examples
...@@ -571,8 +586,13 @@ gev.d.params <- function(fit,ydat=NULL){ ...@@ -571,8 +586,13 @@ gev.d.params <- function(fit,ydat=NULL){
npth <- 0 npth <- 0
}#With no theta parameter, asked by user }#With no theta parameter, asked by user
npet <- length(fit$model[[5]]) + 1 npet <- length(fit$model[[5]]) + 1
if(!fit$eta2_zero){
npe2 <- length(fit$model[[6]]) + 1 #Including eta2 parameter (not default)]
}else{
npe2 <- 0