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

changed initial values input in gev.d.fit function 'init.val' into a lsit of 5...

changed initial values input in gev.d.fit function 'init.val' into a lsit of 5 so that it can now contain multiple values for each parameter (because i noticed that might be necessary after all)
parent a1d270b8
...@@ -24,10 +24,11 @@ ...@@ -24,10 +24,11 @@
#' Parameters are: modified location, scale_0, shape, duration offset, duration exponent repectively. #' Parameters are: modified location, scale_0, shape, duration offset, duration exponent repectively.
#' @param mulink,siglink,shlink,thetalink,etalink Link functions for generalized linear #' @param mulink,siglink,shlink,thetalink,etalink Link functions for generalized linear
#' modelling of the parameters, created with \code{\link{make.link}}. #' modelling of the parameters, created with \code{\link{make.link}}.
#' @param init.vals vector of length 5, giving initial values for parameter intercepts #' @param init.vals list of length 5, giving initial values for all or some parameters
#' used to model the parameters (order: mu, sigma, xi, theta, eta). If rep(NA,5) (the default) is given, initial parameters are obtained #' (order: mu, sigma, 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.
#' @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.
...@@ -75,7 +76,7 @@ gev.d.fit<- ...@@ -75,7 +76,7 @@ gev.d.fit<-
function(xdat, ds, ydat = NULL, mul = NULL, sigl = NULL, shl = NULL, thetal = NULL, etal = NULL, function(xdat, ds, ydat = NULL, mul = NULL, sigl = NULL, shl = NULL, thetal = NULL, etal = NULL,
mulink = make.link("identity"), siglink = make.link("identity"), shlink = make.link("identity"), mulink = make.link("identity"), siglink = make.link("identity"), shlink = make.link("identity"),
thetalink = make.link("identity"), etalink = make.link("identity"), thetalink = make.link("identity"), etalink = make.link("identity"),
init.vals = 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)) {
...@@ -92,7 +93,6 @@ gev.d.fit<- ...@@ -92,7 +93,6 @@ gev.d.fit<-
z$model <- list(mul, sigl, shl, thetal, etal) z$model <- list(mul, sigl, shl, thetal, etal)
z$link <- list(mulink=mulink, siglink=siglink, shlink=shlink, thetalink=thetalink, etalink=etalink) z$link <- list(mulink=mulink, siglink=siglink, shlink=shlink, 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 if covariates matrix is given correctly # test if covariates matrix is given correctly
...@@ -100,27 +100,23 @@ gev.d.fit<- ...@@ -100,27 +100,23 @@ gev.d.fit<-
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) { if(length(init.vals)!=5 | !is.list(init.vals)) {
warning('Parameter init.vals is not used, because it is not of length 5.') warning('Parameter init.vals is not used, because it is no list of length 5.')
init.vals <- rep(NA,5) init.vals <- as.list(rep(NA,5))
} }
if(!any(is.na(init.vals))){ #all initial values are given if(!any(is.na(init.vals))){ #all initial values are given
init.vals <- data.frame(mu = init.vals[1], sigma = init.vals[2], xi = init.vals[3] names(init.vals) <- c('mu','sigma','xi','theta','eta')
,theta = init.vals[4], eta = init.vals[5])
}else if(any(!is.na(init.vals))) { #some initial values are given }else if(any(!is.na(init.vals))) { #some initial values are given
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
if(!is.na(init.vals.user[i])) { if(!is.na(init.vals.user[[i]])) {
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)
} }
if(theta_zero==TRUE) { #if theta should stay zero
init.vals[4] = 0
}
# generate covariates matrices: # generate covariates matrices:
if (is.null(mul)) { #stationary if (is.null(mul)) { #stationary
...@@ -129,7 +125,7 @@ gev.d.fit<- ...@@ -129,7 +125,7 @@ gev.d.fit<-
}else { #non stationary }else { #non stationary
z$trans <- TRUE z$trans <- TRUE
mumat <- cbind(rep(1, length(xdat)), ydat[, mul]) mumat <- cbind(rep(1, length(xdat)), ydat[, mul])
muinit <- c(init.vals$mu, rep(0, length(mul))) muinit <- c(init.vals$mu, rep(0, length(mul)))[1:npmu] #fill with 0s to length npmu
} }
if (is.null(sigl)) { if (is.null(sigl)) {
sigmat <- as.matrix(rep(1, length(xdat))) sigmat <- as.matrix(rep(1, length(xdat)))
...@@ -137,7 +133,7 @@ gev.d.fit<- ...@@ -137,7 +133,7 @@ gev.d.fit<-
}else { }else {
z$trans <- TRUE z$trans <- TRUE
sigmat <- cbind(rep(1, length(xdat)), ydat[, sigl]) sigmat <- cbind(rep(1, length(xdat)), ydat[, sigl])
siginit <- c(init.vals$sigma, rep(0, length(sigl))) siginit <- c(init.vals$sigma, rep(0, length(sigl)))[1:npsc]
} }
if (is.null(shl)) { if (is.null(shl)) {
shmat <- as.matrix(rep(1, length(xdat))) shmat <- as.matrix(rep(1, length(xdat)))
...@@ -145,7 +141,7 @@ gev.d.fit<- ...@@ -145,7 +141,7 @@ gev.d.fit<-
}else { }else {
z$trans <- TRUE z$trans <- TRUE
shmat <- cbind(rep(1, length(xdat)), ydat[, shl]) shmat <- cbind(rep(1, length(xdat)), ydat[, shl])
shinit <- c(init.vals$xi, rep(0, length(shl))) shinit <- c(init.vals$xi, rep(0, length(shl)))[1:npsh]
} }
if (is.null(thetal)) { if (is.null(thetal)) {
thmat <- as.matrix(rep(1, length(xdat))) thmat <- as.matrix(rep(1, length(xdat)))
...@@ -153,7 +149,7 @@ gev.d.fit<- ...@@ -153,7 +149,7 @@ gev.d.fit<-
}else { }else {
z$trans <- TRUE z$trans <- TRUE
thmat <- cbind(rep(1, length(xdat)), ydat[, thetal]) thmat <- cbind(rep(1, length(xdat)), ydat[, thetal])
thetainit <- c(init.vals$theta, rep(0, length(thetal))) thetainit <- c(init.vals$theta, rep(0, length(thetal)))[1:npth]
} }
if (is.null(etal)) { if (is.null(etal)) {
etmat <- as.matrix(rep(1, length(xdat))) etmat <- as.matrix(rep(1, length(xdat)))
...@@ -161,15 +157,17 @@ gev.d.fit<- ...@@ -161,15 +157,17 @@ gev.d.fit<-
}else { }else {
z$trans <- TRUE z$trans <- TRUE
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))) 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)
}else{ #Do not return initial value for theta, if user does not want theta, as Hessian will fail. }else{ #Do not return initial value for theta, if user does not want theta, as Hessian will fail.
init <- c(muinit, siginit, shinit, etainit) init <- c(muinit, siginit, shinit, etainit)
} }
# function to calculate neg log-likelihood: # function to calculate neg log-likelihood:
gev.lik <- function(a) { gev.lik <- function(a) {
# computes neg log lik of d-gev model # computes neg log lik of d-gev model
...@@ -185,7 +183,7 @@ gev.d.fit<- ...@@ -185,7 +183,7 @@ gev.d.fit<-
y <- xdat/sigma.d - mu y <- xdat/sigma.d - mu
y <- 1 + xi * y y <- 1 + xi * y
if(!theta_zero){ #When user does not want 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)
...@@ -228,7 +226,7 @@ gev.d.fit<- ...@@ -228,7 +226,7 @@ gev.d.fit<-
} else {#When theta parameter is not returned, asked by user } else {#When theta parameter is not returned, asked by user
z$vals <- cbind(mut, sc0, xi, eta) z$vals <- cbind(mut, sc0, xi, eta)
} }
z$init.vals <- as.numeric(init.vals) z$init.vals <- unlist(init.vals)
if(!theta_zero){ #When theta parameter is returned (default) if(!theta_zero){ #When theta parameter is returned (default)
colnames(z$vals) <-c('mut','sigma0','xi','theta','eta') colnames(z$vals) <-c('mut','sigma0','xi','theta','eta')
} else { #When theta parameter is not returned, asked by user } else { #When theta parameter is not returned, asked by user
......
...@@ -9,7 +9,7 @@ IDF.agg( ...@@ -9,7 +9,7 @@ IDF.agg(
ds, ds,
na.accept = 0, na.accept = 0,
which.stations = NULL, which.stations = NULL,
which.mon = 0:11, which.mon = list(0:11),
names = c("date", "RR"), names = c("date", "RR"),
cl = NULL cl = NULL
) )
...@@ -23,7 +23,8 @@ standard date format.} ...@@ -23,7 +23,8 @@ standard date format.}
\item{ds}{numeric vector of aggregation durations. \item{ds}{numeric vector of aggregation durations.
(Must be multiples of time resolution at all stations.)} (Must be multiples of time resolution at all stations.)}
\item{na.accept}{numeric giving maximum number of missing values for which annual max. should still be calculated} \item{na.accept}{numeric in [0,1] giving maximum percentage of missing values
for which block max. should still be calculated}
\item{which.stations}{optional, subset of stations. Either numeric vector or character vector \item{which.stations}{optional, subset of stations. Either numeric vector or character vector
containing names of elements in data. If not given, all elements in `data` will be used.} containing names of elements in data. If not given, all elements in `data` will be used.}
......
...@@ -18,7 +18,7 @@ gev.d.fit( ...@@ -18,7 +18,7 @@ gev.d.fit(
shlink = make.link("identity"), shlink = make.link("identity"),
thetalink = make.link("identity"), thetalink = make.link("identity"),
etalink = make.link("identity"), etalink = make.link("identity"),
init.vals = rep(NA, 5), init.vals = as.list(rep(NA, 5)),
theta_zero = FALSE, theta_zero = FALSE,
show = TRUE, show = TRUE,
method = "Nelder-Mead", method = "Nelder-Mead",
...@@ -45,10 +45,11 @@ Parameters are: modified location, scale_0, shape, duration offset, duration exp ...@@ -45,10 +45,11 @@ Parameters are: modified location, scale_0, shape, duration offset, duration exp
\item{mulink, siglink, shlink, thetalink, etalink}{Link functions for generalized linear \item{mulink, siglink, shlink, thetalink, etalink}{Link functions for generalized linear
modelling of the parameters, created with \code{\link{make.link}}.} modelling of the parameters, created with \code{\link{make.link}}.}
\item{init.vals}{vector of length 5, giving initial values for parameter intercepts \item{init.vals}{list of length 5, giving initial values for all or some parameters
used to model the parameters (order: mu, sigma, xi, theta, eta). If rep(NA,5) (the default) is given, initial parameters are obtained (order: mu, sigma, 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.}
\item{theta_zero}{Logical value, indicating if theta should be estimated (FALSE, the default) or \item{theta_zero}{Logical value, indicating if theta should be estimated (FALSE, the default) or
should stay zero.} should stay zero.}
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment