Commit c9422e57 authored by Felix Fauer's avatar Felix Fauer
Browse files

added intensity offset to gev.d.fit and IDF.plot

parent 53c55a35
...@@ -221,14 +221,19 @@ NULL ...@@ -221,14 +221,19 @@ NULL
#' points(example[example$cov1==1,]$d,example[example$cov1==1,]$dat) #' points(example[example$cov1==1,]$d,example[example$cov1==1,]$dat)
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,tau_zero=TRUE,...){
# if cols is to 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))
if(!tau_zero){
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
}
## 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]) theta=fitparams[4],eta=fitparams[5], tau=tau)
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
...@@ -246,7 +251,6 @@ IDF.plot <- function(durations,fitparams,probs=c(0.5,0.9,0.99), ...@@ -246,7 +251,6 @@ IDF.plot <- function(durations,fitparams,probs=c(0.5,0.9,0.99),
# empty plot # empty plot
plot(NA,xlim=xlim,ylim=ylim,xlab="Duration [h]",ylab="Intensity [mm/h]",log="xy",main=main) plot(NA,xlim=xlim,ylim=ylim,xlab="Duration [h]",ylab="Intensity [mm/h]",log="xy",main=main)
} }
## plot IDF curves ## plot IDF curves
for(i in 1:length(probs)){ for(i in 1:length(probs)){
lines(durations,idf.array[i,],col=cols[i],...) lines(durations,idf.array[i,],col=cols[i],...)
......
...@@ -142,7 +142,7 @@ pgev.d <- function(q,mut,sigma0,xi,theta,eta,d,...) { ...@@ -142,7 +142,7 @@ pgev.d <- function(q,mut,sigma0,xi,theta,eta,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,...) { 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)){ 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. ', message('One of the parameters mut, sigma0, xi, theta, eta is a vector. ',
'This is not intended and might cause an error.')} 'This is not intended and might cause an error.')}
...@@ -155,6 +155,21 @@ qgev.d <- function(p,mut,sigma0,xi,theta,eta,d,...) { ...@@ -155,6 +155,21 @@ qgev.d <- function(p,mut,sigma0,xi,theta,eta,d,...) {
sigma.d <-sigma0/(d+theta)^eta sigma.d <-sigma0/(d+theta)^eta
return(qgev(p,loc=as.numeric(mut*sigma.d) return(qgev(p,loc=as.numeric(mut*sigma.d)
,scale=as.numeric(sigma.d),shape=as.numeric(xi),...))} ,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, ...) {
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.')}
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
ifelse(!is.null(tau), sigma.d <-sigma0/(d+theta)^eta+tau, sigma.d <-sigma0/(d+theta)^eta)
return(qgev(p,loc=as.numeric(mut*sigma.d)
,scale=as.numeric(sigma.d),shape=as.numeric(xi),...))}
} }
......
...@@ -74,10 +74,10 @@ ...@@ -74,10 +74,10 @@
#' ,mutl=c(1,2),sigma0l=1) #' ,mutl=c(1,2),sigma0l=1)
gev.d.fit<- gev.d.fit<-
function(xdat, ds, ydat = NULL, mutl = NULL, sigma0l = NULL, xil = NULL, thetal = NULL, etal = NULL, function(xdat, ds, ydat = NULL, mutl = NULL, sigma0l = NULL, xil = NULL, thetal = NULL, etal = NULL, taul = NULL,
mutlink = make.link("identity"), sigma0link = make.link("identity"), xilink = make.link("identity"), mutlink = make.link("identity"), sigma0link = make.link("identity"), xilink = make.link("identity"),
thetalink = make.link("identity"), etalink = make.link("identity"), thetalink = make.link("identity"), etalink = make.link("identity"), taulink = make.link("identity"),
init.vals = as.list(rep(NA,5)), theta_zero = FALSE, init.vals = NULL, theta_zero = FALSE, tau_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)) {
...@@ -90,9 +90,10 @@ gev.d.fit<- ...@@ -90,9 +90,10 @@ gev.d.fit<-
npsh <- length(xil) + 1 npsh <- length(xil) + 1
npth <- ifelse(!theta_zero,length(thetal) + 1,0) npth <- ifelse(!theta_zero,length(thetal) + 1,0)
npet <- length(etal) + 1 npet <- length(etal) + 1
npta <- ifelse(!tau_zero, length(taul) + 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) z$model <- list(mutl, sigma0l, xil, thetal, etal, taul)
z$link <- list(mutlink=mutlink, sigma0link=sigma0link, xilink=xilink, thetalink=thetalink, etalink=etalink) z$link <- list(mutlink=mutlink, sigma0link=sigma0link, xilink=xilink, thetalink=thetalink, etalink=etalink, 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.')
...@@ -104,13 +105,21 @@ gev.d.fit<- ...@@ -104,13 +105,21 @@ 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 | !is.list(init.vals)) { init.necessary.length = 4 + as.numeric(!theta_zero) + as.numeric(!tau_zero) # mut, sigma0, xi, theta, eta1, eta2, tau
warning('Parameter init.vals is not used, because it is no list of length 5.') if (is.null(init.vals)) {init.vals = as.list(rep(NA,init.necessary.length))
init.vals <- as.list(rep(NA,5)) }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 if(!any(is.na(init.vals))){ #all initial values are given
names(init.vals) <- c('mu','sigma','xi','theta','eta') 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
}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 (!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.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
...@@ -163,27 +172,43 @@ gev.d.fit<- ...@@ -163,27 +172,43 @@ gev.d.fit<-
etmat <- cbind(rep(1, length(xdat)), ydat[, etal]) etmat <- cbind(rep(1, length(xdat)), ydat[, etal])
etainit <- c(init.vals$eta, rep(0, length(etal)))[1:npet] etainit <- c(init.vals$eta, rep(0, length(etal)))[1:npet]
} }
if (!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) 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)
} }
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)
# 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
mu <- mutlink$linkinv(mumat %*% (a[1:npmu])) mu <- mutlink$linkinv(mumat %*% (a[1:npmu]))
sigma <- sigma0link$linkinv(sigmat %*% (a[seq(npmu + 1, length = npsc)])) sigma <- sigma0link$linkinv(sigmat %*% (a[seq(npmu + 1, length = npsc)]))
xi <- xilink$linkinv(shmat %*% (a[seq(npmu + npsc + 1, length = npsh)])) xi <- xilink$linkinv(shmat %*% (a[seq(npmu + npsc + 1, length = npsh)]))
# Next line will set the theta likelihood as non-existent in case user requested it. # Next line will set the theta likelihood as non-existent in case user requested it. (same for tau)
if(!theta_zero) {theta <- thetalink$linkinv(thmat %*% (a[seq(npmu + npsc + npsh + 1, length = npth)]))} if(!theta_zero) {theta <- thetalink$linkinv(thmat %*% (a[seq(npmu + npsc + npsh + 1, length = npth)]))}
eta <- etalink$linkinv(etmat %*% (a[seq(npmu + npsc + npsh + npth + 1, length = npet)])) eta <- etalink$linkinv(etmat %*% (a[seq(npmu + npsc + npsh + npth + 1, length = npet)]))
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(!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)
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 <- xdat/sigma.d - mu
y <- 1 + xi * y y <- 1 + xi * y
...@@ -192,7 +217,8 @@ gev.d.fit<- ...@@ -192,7 +217,8 @@ gev.d.fit<-
}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)
} }
if(!tau_zero) {if(any(tau<=0)) return(10^6)} # check condition for tau as well.
sum(log(sigma.d)) + sum(y^(-1/xi)) + sum(log(y) * (1/xi + 1)) sum(log(sigma.d)) + sum(y^(-1/xi)) + sum(log(y) * (1/xi + 1))
} }
...@@ -212,9 +238,15 @@ gev.d.fit<- ...@@ -212,9 +238,15 @@ 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(!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
z$nllh <- x$value z$nllh <- x$value
# normalize data to standard Gumbel: # normalize data to standard Gumbel:
sc.d <- sc0/((ds+theta)^eta) #sc.d <- sc0/((ds+theta)^eta) # old
ifelse(!tau_zero, sc.d <- sc0/((ds+theta)^eta)+tau, sc.d <- sc0/((ds+theta)^eta))
z$data <- - log(as.vector((1 + xi * (xdat/sc.d-mut))^(-1/xi))) z$data <- - log(as.vector((1 + xi * (xdat/sc.d-mut))^(-1/xi)))
z$mle <- x$par z$mle <- x$par
test <- try( # catch error test <- try( # catch error
...@@ -226,24 +258,39 @@ gev.d.fit<- ...@@ -226,24 +258,39 @@ gev.d.fit<-
} }
z$se <- sqrt(diag(z$cov)) # sqrt(digonal entries) = standart error of mle's z$se <- sqrt(diag(z$cov)) # sqrt(digonal entries) = standart error of mle's
if (!theta_zero) {#When theta parameter is returned (default) if (!theta_zero) {#When theta parameter is returned (default)
z$vals <- cbind(mut, sc0, xi, theta, eta) 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 } 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$init.vals <- unlist(init.vals) z$init.vals <- unlist(init.vals)
if(!theta_zero){ #When theta parameter is returned (default)
colnames(z$vals) <-c('mut','sigma0','xi','theta','eta') names2 = c('mut','sigma0','xi') # fixed part for set of names
} else { #When theta parameter is not returned, asked by user if(!theta_zero){names2=c(names2,'theta')} # add theta (in case)
colnames(z$vals) <-c('mut','sigma0','xi','eta') names2 = c(names2, 'eta') # add eta (always)
} if(!tau_zero){names2=c(names2, 'tau')} # add tau (in case)
colnames(z$vals) <- names2
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.
if(show) { if(show) {
if(z$trans) { # for nonstationary fit if(z$trans) { # for nonstationary fit
print(z[c(2, 4)]) # print model, link (3) , conv print(z[c(2, 4)]) # print model, link (3) , conv
# print names of link functions: # print names of link functions:
cat('$link\n') cat('$link\n')
print(c(z$link$mutlink$name,z$link$sigma0link$name,z$link$xilink$name,z$link$thetalink$name,z$link$etalink$name)) 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') cat('\n')
}else{print(z[4])} # for stationary fit print only conv }else{print(z[4])} # for stationary fit print only conv
if(!z$conv){ # if fit converged if(!z$conv){ # if fit converged
......
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