Commit 132b50e7 authored by Felix Fauer's avatar Felix Fauer
Browse files

adding tau to init values

parent c9422e57
......@@ -203,6 +203,7 @@ NULL
#' @param cols vector of colors for IDF curves. Should have same length as \code{probs}
#' @param add logical indicating if plot should be added to existing plot, default is FALSE
#' @param legend logical indicating if legend should be plotted (TRUE, the default)
#' @param tau_zero Logical value, indicating, whether fitparams contains the parameter tau.
#' @param ... additional parameters passed on to the \code{plot} function
#'
#' @export
......@@ -226,7 +227,7 @@ IDF.plot <- function(durations,fitparams,probs=c(0.5,0.9,0.99),
# if cols is to short, make longer
if(length(cols)!=length(probs))cols <- rep_len(cols,length.out=length(probs))
if(!tau_zero){
print("WARNING in IDF.plot: this might work now, but is not robust any more when multiscaling is added")
#print("WARNING in IDF.plot: this might work now, but is not robust any more when multiscaling is added")
tau=fitparams[6]
}else{
tau=NULL
......
......@@ -105,8 +105,9 @@ 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 tau numeric value, giving intensity offset (defining flattening of the IDF curve for long durations)
#' @param d positive numeric value, giving duration
#' @param ... additional parameters passed to \code{\link[evd]{qgev}}
#'
......@@ -129,11 +130,11 @@ 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, tau=0.1,d=1)
#'
#' # calculate quantiles for sequence of durations
#' ds <- 2^seq(0,4,0.1)
#' qs <- lapply(ds,qgev.d,p=p,mut=4,sigma0=2,xi=0,theta=0.1,eta=0.3)
#' qs <- lapply(ds,qgev.d,p=p,mut=4,sigma0=2,xi=0,theta=0.1,eta=0.3, tau=0.1)
#' qs <- simplify2array(qs)
#'
#' plot(ds,qs[1,],ylim=c(3,20),type='l',log = 'xy',ylab='Intensity',xlab = 'Duration')
......@@ -142,21 +143,7 @@ pgev.d <- function(q,mut,sigma0,xi,theta,eta,d,...) {
#' }
#' legend('topright',title = 'p-quantile',
#' legend = p,lty=1:3,bty = 'n')
qgev.d.old <- function(p,mut,sigma0,xi,theta,eta,d,...) { # cannot deal with tau
if(any(c(length(mut),length(sigma0),length(xi),length(theta),length(eta))>1)){
message('One of the parameters mut, sigma0, xi, theta, eta is a vector. ',
'This is not intended and might cause an error.')}
if (d<=0) {stop('The duration d has to be positive.')}
if(any(d+theta<=0)){
warning('Some shape parameters are negative, resulting from a negativ theta '
,theta, ' this will prododuce NAs.')}
# 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)
,scale=as.numeric(sigma.d),shape=as.numeric(xi),...))}
} # can deal with tau!
qgev.d <- function(p,mut,sigma0,xi,theta,eta,d, tau=NULL, ...) {
qgev.d <- function(p,mut,sigma0,xi,theta,eta,tau,d, ...) {
if(any(c(length(mut),length(sigma0),length(xi),length(theta),length(eta), length(tau))>1)){
message('One of the parameters mut, sigma0, xi, theta, eta is a vector. ',
'This is not intended and might cause an error.')}
......
......@@ -18,18 +18,20 @@
#' @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 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,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
#' 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 Logical value, indicating whether tau should be estimated (TRUE, the default) or
#' should stay zero.
#' @param show Logical; if TRUE (the default), print details of the fit.
#' @param method The optimization method used in \code{\link{optim}}.
......@@ -105,31 +107,41 @@ gev.d.fit<-
if(any(npar>ncol(ydat),npar>0 & is.null(ydat)))stop("Not enough columns in covariates matrix 'ydat'.")
# initial values
init.necessary.length = 4 + as.numeric(!theta_zero) + as.numeric(!tau_zero) # mut, sigma0, xi, theta, eta1, eta2, tau
init.necessary.length = 4 + as.numeric(!theta_zero) + as.numeric(!tau_zero) # mut, sigma0, xi, theta, eta, 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)) {
warning(paste0('Parameter init.vals is not used, because it is no list of length ',init.necessary.length,'.'))
init.vals <- as.list(rep(NA,init.necessary.length))
}
if(!any(is.na(init.vals))){ #all initial values are given
names1=c('mu','sigma','xi','theta','eta')
if (!tau_zero){names1=c(names1,'tau')}
names(init.vals) <- names1 #c('mu','sigma','xi','theta','eta') #old
init.vals <- gev.d.init(xdat,ds,z$link)
}else if(any(!is.na(init.vals))) { #some initial values are given
if (!tau_zero){print("WARNING. The automatic estimation of init.vals is not implentended yet for tau_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[[i]])) {
init.vals[[i]]<-init.vals.user[[i]]
}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 (!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 (tau_zero) init.vals$tau = 0
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 (!tau_zero){print("WARNING. The automatic estimation of init.vals is not implentended yet for tau_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
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
......@@ -172,26 +184,21 @@ gev.d.fit<-
etmat <- cbind(rep(1, length(xdat)), ydat[, etal])
etainit <- c(init.vals$eta, rep(0, length(etal)))[1:npet]
}
if (!tau_zero){
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 (!tau_zero){
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(!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)
}
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 (!tau_zero) {init = c(init,tauinit)} # add tau init (in case)
if (!theta_zero) {init <- c(init,thetainit)} # add theta init (in case)
init <- c(init,etainit) # add eta init (always)
if (!tau_zero) {init <- c(init,tauinit)} # add tau init (in case)
# function to calculate neg log-likelihood:
gev.lik <- function(a) {
......@@ -205,19 +212,23 @@ gev.d.fit<-
if(!tau_zero) {tau <- taulink$linkinv( tamat %*% (a[seq(npmu + npsc + npsh + npth + npet + 1, length = npta)]))}
ifelse(!theta_zero, ds.t <- ds+theta, ds.t <- ds) #Don't use theta if user requested not to have it.
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.
#sigma.d <- sigma/(ds.t^eta)
y <- xdat/sigma.d - mu
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){ #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(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.
sum(log(sigma.d)) + sum(y^(-1/xi)) + sum(log(y) * (1/xi + 1))
}
......@@ -240,13 +251,14 @@ gev.d.fit<-
eta <- etalink$linkinv(etmat %*% (x$par[seq(npmu + npsc + npsh + npth + 1, length = npet)]))
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)]))
} # do nothing when user requests tau to be zero
}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) # old
ifelse(!tau_zero, sc.d <- sc0/((ds+theta)^eta)+tau, sc.d <- sc0/((ds+theta)^eta))
sc.d <- sc0/((ds+theta)^eta)+tau
z$data <- - log(as.vector((1 + xi * (xdat/sc.d-mut))^(-1/xi)))
z$mle <- x$par
test <- try( # catch error
......@@ -257,7 +269,7 @@ 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)
'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
......@@ -269,14 +281,16 @@ gev.d.fit<-
}else{ # if tau is not returned
z$vals <- cbind(mut, sc0, xi, eta)
}
}
}'
z$vals <- cbind(mut, sc0, xi, theta, eta, tau)
z$init.vals <- unlist(init.vals)
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
'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","tau")
z$ds <- ds
z$theta_zero <- theta_zero #Indicates if theta parameter was set to zero by user.
......@@ -286,11 +300,11 @@ gev.d.fit<-
print(z[c(2, 4)]) # print model, link (3) , conv
# print names of link functions:
cat('$link\n')
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))
} else{
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$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
......@@ -348,8 +362,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, tau=tauinit))
}
#### gev.d.lik ####
......
......@@ -11,6 +11,7 @@ IDF.plot(
cols = 4:2,
add = FALSE,
legend = TRUE,
tau_zero = TRUE,
...
)
}
......@@ -30,6 +31,8 @@ as obtained from \code{\link{gev.d.fit}}
\item{legend}{logical indicating if legend should be plotted (TRUE, the default)}
\item{tau_zero}{Logical value, indicating, whether fitparams contains the parameter tau.}
\item{...}{additional parameters passed on to the \code{plot} function}
}
\description{
......
......@@ -13,13 +13,16 @@ gev.d.fit(
xil = NULL,
thetal = NULL,
etal = NULL,
taul = 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)),
taulink = make.link("identity"),
init.vals = NULL,
theta_zero = FALSE,
tau_zero = TRUE,
show = TRUE,
method = "Nelder-Mead",
maxit = 10000,
......@@ -37,12 +40,12 @@ This can be obtained from \code{\link{IDF.agg}}.}
(or NULL (the default) for stationary fitting). The number of rows should be the same as the
length of xdat.}
\item{mutl, sigma0l, xil, thetal, etal}{Numeric vectors of integers, giving the columns of ydat that contain
\item{mutl, sigma0l, xil, thetal, etal, taul}{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.}
\item{mutlink, sigma0link, xilink, thetalink, etalink}{Link functions for generalized linear
\item{mutlink, sigma0link, xilink, thetalink, etalink, taulink}{Link functions for generalized linear
modeling of the parameters, created with \code{\link{make.link}}. The default is \code{make.link("identity")}.}
\item{init.vals}{list of length 5, giving initial values for all or some parameters
......@@ -51,7 +54,10 @@ internally by fitting the GEV separately for each duration and applying a linear
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 whether theta should be estimated (FALSE, the default) or
should stay zero.}
\item{tau_zero}{Logical value, indicating whether tau should be estimated (TRUE, the default) or
should stay zero.}
\item{show}{Logical; if TRUE (the default), print details of the fit.}
......
......@@ -4,17 +4,19 @@
\alias{qgev.d}
\title{d-GEV quantile function}
\usage{
qgev.d(p, mut, sigma0, xi, theta, eta, d, ...)
qgev.d(p, mut, sigma0, xi, theta, eta, tau, d, ...)
}
\arguments{
\item{p}{vector of probabilities}
\item{mut, sigma0, xi}{numeric value, giving modified location, modified scale and shape parameter}
\item{theta}{numeric value, giving duration offset (defining curvature of the IDF curve)}
\item{theta}{numeric value, giving duration offset (defining curvature of the IDF curve for short durations)}
\item{eta}{numeric value, giving duration exponent (defining slope of the IDF curve)}
\item{tau}{numeric value, giving intensity offset (defining flattening of the IDF curve for long durations)}
\item{d}{positive numeric value, giving duration}
\item{...}{additional parameters passed to \code{\link[evd]{qgev}}}
......@@ -38,11 +40,11 @@ For details on the d-GEV and the parameter definitions, see \link{IDF-package}.
\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, tau=0.1,d=1)
# calculate quantiles for sequence of durations
ds <- 2^seq(0,4,0.1)
qs <- lapply(ds,qgev.d,p=p,mut=4,sigma0=2,xi=0,theta=0.1,eta=0.3)
qs <- lapply(ds,qgev.d,p=p,mut=4,sigma0=2,xi=0,theta=0.1,eta=0.3, tau=0.1)
qs <- simplify2array(qs)
plot(ds,qs[1,],ylim=c(3,20),type='l',log = 'xy',ylab='Intensity',xlab = 'Duration')
......
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