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

finished adding tau parameter to all functions. all parameters are allways...

finished adding tau parameter to all functions. all parameters are allways stored in fit (except in fit) .also created IDF.plot.fit() (but it is very slow)
parent 132b50e7
...@@ -2,6 +2,7 @@ ...@@ -2,6 +2,7 @@
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)
......
...@@ -203,7 +203,6 @@ NULL ...@@ -203,7 +203,6 @@ NULL
#' @param cols vector of colors for IDF curves. Should have same length as \code{probs} #' @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 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 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 #' @param ... additional parameters passed on to the \code{plot} function
#' #'
#' @export #' @export
...@@ -222,19 +221,13 @@ NULL ...@@ -222,19 +221,13 @@ 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,tau_zero=TRUE,...){ legend=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], tau=tau) theta=fitparams[4],eta=fitparams[5], tau=fitparams[6])
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
...@@ -270,3 +263,26 @@ IDF.plot <- function(durations,fitparams,probs=c(0.5,0.9,0.99), ...@@ -270,3 +263,26 @@ IDF.plot <- function(durations,fitparams,probs=c(0.5,0.9,0.99),
col=cols,lty=lty,lwd=lwd) col=cols,lty=lty,lwd=lwd)
} }
} }
#### IDF.plot.fit ####
#' This is a faster usable version of \code{\link{IDF.plot}}. The only argument needed is a fit object returned by \code{\link{gev.d.fit}}
#'
#' @param fit A fit object returned by \code{\link{gev.d.fit}}
#' @param ... Options to be passed to \code{\link{IDF.plot}}
#' @export
#' @importFrom
#' @example
#' #' data('example',package = 'IDF')
#' # fit d-gev
#' fit <- gev.d.fit(example$dat,example$d,ydat = as.matrix(example[,c("cov1","cov2")])
#' ,mutl = c(1,2),sigma0l = 1)
#' # plot quantiles
#' IDF.plot.fit(fit)
#' # add data points
#' points(example[example$cov1==1,]$d,example[example$cov1==1,]$dat)
IDF.plot.fit <- function(fit){
fitted_params=gev.d.params(fit,...)
ds = fit$ds
IDF.plot(ds, fitted_params,...)
}
\ No newline at end of file
...@@ -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 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}}
#' #'
...@@ -26,20 +27,20 @@ ...@@ -26,20 +27,20 @@
#' @examples #' @examples
#' x <- seq(4,20,0.1) #' x <- seq(4,20,0.1)
#' # calculate probability density for one duration #' # calculate probability density for one duration
#' dgev.d(q=x,mut=4,sigma0=2,xi=0,theta=0.1,eta=0.1,d=1) #' dgev.d(q=x,mut=4,sigma0=2,xi=0,theta=0.1,eta=0.1,tau=0.1,d=1)
#' #'
#' # calculate probability density for different durations #' # calculate probability density for different durations
#' ds <- 1:4 #' ds <- 1:4
#' dens <- lapply(ds,dgev.d,q=x,mut=4,sigma0=2,xi=0,theta=0.1,eta=0.1) #' dens <- lapply(ds,dgev.d,q=x,mut=4,sigma0=2,xi=0,theta=0.1,eta=0.1, tau=0.1)
#' #'
#' plot(x,dens[[1]],type='l',ylim = c(0,0.21),ylab = 'Probability Density') #' plot(x,dens[[1]],type='l',ylim = c(0,0.21),ylab = 'Probability Density')
#' for(i in 2:4){ #' for(i in 2:4){
#' 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,...) { dgev.d <- function(q,mut,sigma0,xi,theta,eta,tau,d,...) {
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), length(tau))>1)){
message('One of the parameters mut, sigma0, xi, theta, eta 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.')}
if(any(d+theta<=0)){ if(any(d+theta<=0)){
...@@ -47,7 +48,7 @@ dgev.d <- function(q,mut,sigma0,xi,theta,eta,d,...) { ...@@ -47,7 +48,7 @@ dgev.d <- function(q,mut,sigma0,xi,theta,eta,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 sigma.d <-sigma0/(d+theta)^eta+ tau
return(dgev(q,loc=mut*sigma.d,scale=sigma.d,shape=xi,...))} return(dgev(q,loc=mut*sigma.d,scale=sigma.d,shape=xi,...))}
} }
...@@ -61,6 +62,7 @@ dgev.d <- function(q,mut,sigma0,xi,theta,eta,d,...) { ...@@ -61,6 +62,7 @@ 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 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 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}}
#' #'
...@@ -82,10 +84,10 @@ dgev.d <- function(q,mut,sigma0,xi,theta,eta,d,...) { ...@@ -82,10 +84,10 @@ dgev.d <- function(q,mut,sigma0,xi,theta,eta,d,...) {
#' #'
#' @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,tau=0.1,d=1)
pgev.d <- function(q,mut,sigma0,xi,theta,eta,d,...) { pgev.d <- function(q,mut,sigma0,xi,theta,eta,tau,d,...) {
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),length(tau))>1)){
message('One of the parameters mut, sigma0, xi, theta, eta 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.')}
if(any(d+theta<=0)){ if(any(d+theta<=0)){
...@@ -93,7 +95,7 @@ pgev.d <- function(q,mut,sigma0,xi,theta,eta,d,...) { ...@@ -93,7 +95,7 @@ pgev.d <- function(q,mut,sigma0,xi,theta,eta,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 sigma.d <-sigma0/(d+theta)^eta+tau
return(pgev(q,loc=mut*sigma.d,scale=sigma.d,shape=xi,...))} return(pgev(q,loc=mut*sigma.d,scale=sigma.d,shape=xi,...))}
} }
...@@ -154,7 +156,7 @@ qgev.d <- function(p,mut,sigma0,xi,theta,eta,tau,d, ...) { ...@@ -154,7 +156,7 @@ qgev.d <- function(p,mut,sigma0,xi,theta,eta,tau,d, ...) {
# 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
ifelse(!is.null(tau), sigma.d <-sigma0/(d+theta)^eta+tau, sigma.d <-sigma0/(d+theta)^eta) sigma.d <-sigma0/(d+theta)^eta+tau
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),...))}
} }
...@@ -169,6 +171,7 @@ qgev.d <- function(p,mut,sigma0,xi,theta,eta,tau,d, ...) { ...@@ -169,6 +171,7 @@ qgev.d <- function(p,mut,sigma0,xi,theta,eta,tau,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) #' @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 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
#' #'
#' @details For details on the d-GEV and the parameter definitions, see \link{IDF-package} #' @details For details on the d-GEV and the parameter definitions, see \link{IDF-package}
...@@ -184,20 +187,20 @@ qgev.d <- function(p,mut,sigma0,xi,theta,eta,tau,d, ...) { ...@@ -184,20 +187,20 @@ qgev.d <- function(p,mut,sigma0,xi,theta,eta,tau,d, ...) {
#' #'
#' @examples #' @examples
#' # random sample for one duration #' # random sample for one duration
#' rgev.d(n=100,mut=4,sigma0=2,xi=0,theta=0.1,eta=0.3,d=1) #' rgev.d(n=100,mut=4,sigma0=2,xi=0,theta=0.1,eta=0.3,tau=0.1,d=1)
#' #'
#' # compare randomn samples for different durations #' # compare randomn samples for different durations
#' ds <- c(1,4) #' ds <- c(1,4)
#' samp <- lapply(ds,rgev.d,n=100,mut=4,sigma0=2,xi=0,theta=0.1,eta=0.3) #' samp <- lapply(ds,rgev.d,n=100,mut=4,sigma0=2,xi=0,theta=0.1,eta=0.3,tau=0.1)
#' #'
#' hist(samp[[1]],breaks = 10,col=rgb(1,0,0,0.5),freq = FALSE #' hist(samp[[1]],breaks = 10,col=rgb(1,0,0,0.5),freq = FALSE
#' ,ylim=c(0,0.3),xlim=c(3,20),xlab='x',main = 'Random d-GEV samples') #' ,ylim=c(0,0.3),xlim=c(3,20),xlab='x',main = 'Random d-GEV samples')
#' 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) { rgev.d <- function(n,mut,sigma0,xi,theta,eta,tau,d) {
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),length(tau))>1)){
message('One of the parameters mut, sigma0, xi, theta, eta 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.')}
if(any(d+theta<=0)){ if(any(d+theta<=0)){
...@@ -205,7 +208,7 @@ rgev.d <- function(n,mut,sigma0,xi,theta,eta,d) { ...@@ -205,7 +208,7 @@ rgev.d <- function(n,mut,sigma0,xi,theta,eta,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,n))}else{ if(d+theta<=0){return(rep(NA,n))}else{
sigma.d <-sigma0/(d+theta)^eta sigma.d <-sigma0/(d+theta)^eta+tau
return(rgev(n,loc=mut*sigma.d,scale=sigma.d,shape=xi))} return(rgev(n,loc=mut*sigma.d,scale=sigma.d,shape=xi))}
} }
......
...@@ -130,7 +130,6 @@ gev.d.fit<- ...@@ -130,7 +130,6 @@ 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 (!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
...@@ -374,7 +373,7 @@ gev.d.init <- function(xdat,ds,link){ ...@@ -374,7 +373,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 numeric vectors containing corresponding estimates for each of the parameters #' @param mut,sigma0,xi,theta,eta,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
...@@ -388,16 +387,16 @@ gev.d.init <- function(xdat,ds,link){ ...@@ -388,16 +387,16 @@ gev.d.init <- function(xdat,ds,link){
#' ,ydat = as.matrix(train.set[c('cov1','cov2')])) #' ,ydat = as.matrix(train.set[c('cov1','cov2')]))
#' 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],tau=params[,6],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,tau,log=FALSE) {
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)) %in% if(any(! c(length(ds),length(mut),length(sigma0),length(xi),length(theta),length(eta),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) sigma.d <- sigma0/(ds.t^eta) + tau
y <- xdat/sigma.d - mut y <- xdat/sigma.d - mut
y <- 1 + xi * y y <- 1 + xi * y
...@@ -533,9 +532,17 @@ gev.d.params <- function(fit,ydat=NULL){ ...@@ -533,9 +532,17 @@ gev.d.params <- function(fit,ydat=NULL){
npsc <- length(fit$model[[2]]) + 1 npsc <- length(fit$model[[2]]) + 1
npsh <- length(fit$model[[3]]) + 1 npsh <- length(fit$model[[3]]) + 1
if(class(fit)=="gev.d.fit"){ if(class(fit)=="gev.d.fit"){
if(!fit$theta_zero){npth <- length(fit$model[[4]]) + 1 #Including theta parameter (default)] if(!fit$theta_zero){
}else{npth <- 0}#With no theta parameter, asked by user npth <- length(fit$model[[4]]) + 1 #Including theta parameter (default)]
}else{
npth <- 0
}#With no theta parameter, asked by user
npet <- length(fit$model[[5]]) + 1 npet <- length(fit$model[[5]]) + 1
if(!fit$tau_zero){
npta <- length(fit$model[[6]]) + 1 #Including tau parameter (not default)]
}else{
npta <- 0
}#With no tau parameter, asked by user
} }
# inverse link functions # inverse link functions
...@@ -545,6 +552,7 @@ gev.d.params <- function(fit,ydat=NULL){ ...@@ -545,6 +552,7 @@ gev.d.params <- function(fit,ydat=NULL){
shlink <- fit$link$xilink$linkinv shlink <- fit$link$xilink$linkinv
if(!fit$theta_zero) thetalink <- fit$link$thetalink$linkinv if(!fit$theta_zero) thetalink <- fit$link$thetalink$linkinv
etalink <- fit$link$etalink$linkinv etalink <- fit$link$etalink$linkinv
taulink <- fit$link$taulink$linkinv
}else{ }else{
mulink <- eval(parse(text=fit$link))[[1]] mulink <- eval(parse(text=fit$link))[[1]]
siglink <- eval(parse(text=fit$link))[[2]] siglink <- eval(parse(text=fit$link))[[2]]
...@@ -558,20 +566,29 @@ gev.d.params <- function(fit,ydat=NULL){ ...@@ -558,20 +566,29 @@ gev.d.params <- function(fit,ydat=NULL){
if(class(fit)=="gev.d.fit"){ if(class(fit)=="gev.d.fit"){
if(!fit$theta_zero){thmat <- cbind(rep(1, dim(ydat)[1]), matrix(ydat[, fit$model[[4]]],dim(ydat)[1],npth-1))} if(!fit$theta_zero){thmat <- cbind(rep(1, dim(ydat)[1]), matrix(ydat[, fit$model[[4]]],dim(ydat)[1],npth-1))}
etmat <- cbind(rep(1, dim(ydat)[1]), matrix(ydat[, fit$model[[5]]],dim(ydat)[1],npet-1)) etmat <- cbind(rep(1, dim(ydat)[1]), matrix(ydat[, fit$model[[5]]],dim(ydat)[1],npet-1))
if(!fit$tau_zero) {tamat <- cbind(rep(1, dim(ydat)[1]), matrix(ydat[, fit$model[[6]]],dim(ydat)[1],npta-1))}
} }
# calculate parameters # calculate parameters
mut <- mulink(mumat %*% (fit$mle[1:npmu])) mut <- mulink(mumat %*% (fit$mle[1:npmu]))
sc0 <- siglink(sigmat %*% (fit$mle[seq(npmu + 1, length = npsc)])) sc0 <- siglink(sigmat %*% (fit$mle[seq(npmu + 1, length = npsc)]))
xi <- shlink(shmat %*% (fit$mle[seq(npmu + npsc + 1, length = npsh)])) xi <- shlink(shmat %*% (fit$mle[seq(npmu + npsc + 1, length = npsh)]))
if(class(fit)=="gev.d.fit" ){
if(!fit$theta_zero){theta <- thetalink(thmat %*% (fit$mle[seq(npmu + npsc + npsh + 1, length = npth)]))
}else{theta <- rep(0,dim(ydat)[1])}}
if(class(fit)=="gev.d.fit"){eta <- etalink(etmat %*% (fit$mle[seq(npmu + npsc + npsh + npth + 1, length = npet)]))}
if(class(fit)=="gev.d.fit"){ if(class(fit)=="gev.d.fit"){
return(data.frame(mut=mut,sigma0=sc0,xi=xi,theta=theta,eta=eta)) if(!fit$theta_zero){
}else{return(data.frame(mu=mut,sig=sc0,xi=xi))} theta <- thetalink(thmat %*% (fit$mle[seq(npmu + npsc + npsh + 1, length = npth)]))
}else{
theta <- rep(0,dim(ydat)[1])
}
eta <- etalink(etmat %*% (fit$mle[seq(npmu + npsc + npsh + npth + 1, length = npet)]))
if(!fit$tau_zero){
tau <- taulink(tamat %*% (fit$mle[seq(npmu + npsc + npsh + npth + npet + 1, length = npta)]))
}else{
tau <- rep(0,dim(ydat)[1])
}
return(data.frame(mut=mut,sigma0=sc0,xi=xi,theta=theta,eta=eta, tau=tau))
}else{
return(data.frame(mu=mut,sig=sc0,xi=xi))
}
} }
...@@ -591,7 +608,8 @@ gev.d.params <- function(fit,ydat=NULL){ ...@@ -591,7 +608,8 @@ gev.d.params <- function(fit,ydat=NULL){
#' \item \eqn{\sigma_0 = 2+0.5 cov_1} #' \item \eqn{\sigma_0 = 2+0.5 cov_1}
#' \item \eqn{\xi = 0.5} #' \item \eqn{\xi = 0.5}
#' \item \eqn{\theta = 0} #' \item \eqn{\theta = 0}
#' \item \eqn{\eta = 0.5}} #' \item \eqn{\eta = 0.5}
#' \item \eqn{\tau = 0}}
#' #'
#' #'
#' @docType data #' @docType data
......
...@@ -66,7 +66,7 @@ plot(ann.max$ds,ann.max$xdat,log='xy',xlab = 'Duration [h]',ylab='Intensity [mm/ ...@@ -66,7 +66,7 @@ plot(ann.max$ds,ann.max$xdat,log='xy',xlab = 'Duration [h]',ylab='Intensity [mm/
* Step 2: fit d-GEV to annual maxima * Step 2: fit d-GEV to annual maxima
```{r fit} ```{r fit}
fit <- gev.d.fit(xdat = ann.max$xdat,ds = ann.max$ds,sigma0link = make.link('log')) fit <- gev.d.fit(xdat = ann.max$xdat,ds = ann.max$ds,sigma0link = make.link('log'), tau_zero = FALSE)
# checking the fit # checking the fit
gev.d.diag(fit,pch=1,) gev.d.diag(fit,pch=1,)
# parameter estimates # parameter estimates
...@@ -83,7 +83,7 @@ for(d in 1:2){ # d=1h and d=2h ...@@ -83,7 +83,7 @@ for(d in 1:2){ # d=1h and d=2h
hist(ann.max$xdat[ann.max$ds==d],main = paste('d=',d),q.min:q.max hist(ann.max$xdat[ann.max$ds==d],main = paste('d=',d),q.min:q.max
,freq = FALSE,add=TRUE,border = d) ,freq = FALSE,add=TRUE,border = d)
# etimated prob. density: # etimated prob. density:
lines(q,dgev.d(q,params$mut,params$sigma0,params$xi,params$theta,params$eta,d = d),col=d) lines(q,dgev.d(q,params$mut,params$sigma0,params$xi,params$theta,params$eta,params$tau,d = d),col=d)
} }
legend('topright',col=1:2,lwd=1,legend = paste('d=',1:2,'h'),title = 'Duration') legend('topright',col=1:2,lwd=1,legend = paste('d=',1:2,'h'),title = 'Duration')
``` ```
......
...@@ -11,7 +11,6 @@ IDF.plot( ...@@ -11,7 +11,6 @@ IDF.plot(
cols = 4:2, cols = 4:2,
add = FALSE, add = FALSE,
legend = TRUE, legend = TRUE,
tau_zero = TRUE,
... ...
) )
} }
...@@ -31,8 +30,6 @@ as obtained from \code{\link{gev.d.fit}} ...@@ -31,8 +30,6 @@ as obtained from \code{\link{gev.d.fit}}
\item{legend}{logical indicating if legend should be plotted (TRUE, the default)} \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} \item{...}{additional parameters passed on to the \code{plot} function}
} }
\description{ \description{
......
% Generated by roxygen2: do not edit by hand
% Please edit documentation in R/IDF.R
\name{IDF.plot.fit}
\alias{IDF.plot.fit}
\title{This is a faster usable version of \code{\link{IDF.plot}}. The only argument needed is a fit object returned by \code{\link{gev.d.fit}}}
\usage{
IDF.plot.fit(fit, ...)
}
\arguments{
\item{fit}{A fit object returned by \code{\link{gev.d.fit}}}
\item{...}{Options to be passed to \code{\link{IDF.plot}}}
}
\description{
This is a faster usable version of \code{\link{IDF.plot}}. The only argument needed is a fit object returned by \code{\link{gev.d.fit}}
}
...@@ -4,7 +4,7 @@ ...@@ -4,7 +4,7 @@
\alias{dgev.d} \alias{dgev.d}
\title{d-GEV probability density function} \title{d-GEV probability density function}
\usage{ \usage{
dgev.d(q, mut, sigma0, xi, theta, eta, d, ...) dgev.d(q, mut, sigma0, xi, theta, eta, tau, d, ...)
} }
\arguments{ \arguments{
\item{q}{vector of quantiles} \item{q}{vector of quantiles}
...@@ -16,6 +16,8 @@ shape parameter \eqn{\xi}.} ...@@ -16,6 +16,8 @@ shape parameter \eqn{\xi}.}
\item{eta}{numeric value, giving duration exponent \eqn{\eta} (defining slope of the IDF curve)} \item{eta}{numeric value, giving duration exponent \eqn{\eta} (defining slope of the IDF curve)}
\item{tau}{numeric value, giving intensity offset \eqn{\tau} (defining flattening of the IDF curve)}
\item{d}{positive numeric value, giving duration} \item{d}{positive numeric value, giving duration}
\item{...}{additional parameters passed to \code{\link[evd]{dgev}}} \item{...}{additional parameters passed to \code{\link[evd]{dgev}}}
...@@ -33,11 +35,11 @@ For details on the d-GEV and the parameter definitions, see \link{IDF-package}. ...@@ -33,11 +35,11 @@ For details on the d-GEV and the parameter definitions, see \link{IDF-package}.
\examples{ \examples{
x <- seq(4,20,0.1) x <- seq(4,20,0.1)
# calculate probability density for one duration # calculate probability density for one duration
dgev.d(q=x,mut=4,sigma0=2,xi=0,theta=0.1,eta=0.1,d=1) dgev.d(q=x,mut=4,sigma0=2,xi=0,theta=0.1,eta=0.1,tau=0.1,d=1)
# calculate probability density for different durations # calculate probability density for different durations
ds <- 1:4 ds <- 1:4
dens <- lapply(ds,dgev.d,q=x,mut=4,sigma0=2,xi=0,theta=0.1,eta=0.1) dens <- lapply(ds,dgev.d,q=x,mut=4,sigma0=2,xi=0,theta=0.1,eta=0.1, tau=0.1)
plot(x,dens[[1]],type='l',ylim = c(0,0.21),ylab = 'Probability Density') plot(x,dens[[1]],type='l',ylim = c(0,0.21),ylab = 'Probability Density')
for(i in 2:4){ for(i in 2:4){
......
...@@ -22,6 +22,7 @@ d-GEV parameters used for sampling: ...@@ -22,6 +22,7 @@ d-GEV parameters used for sampling:
\item \eqn{\sigma_0 = 2+0.5 cov_1} \item \eqn{\sigma_0 = 2+0.5 cov_1}
\item \eqn{\xi = 0.5} \item \eqn{\xi = 0.5}
\item \eqn{\theta = 0} \item \eqn{\theta = 0}
\item \eqn{\eta = 0.5}} \item \eqn{\eta = 0.5}
\item \eqn{\tau = 0}}
} }
\keyword{datasets} \keyword{datasets}
...@@ -4,14 +4,14 @@ ...@@ -4,14 +4,14 @@
\alias{gev.d.lik} \alias{gev.d.lik}
\title{d-GEV Likelihood} \title{d-GEV Likelihood}
\usage{ \usage{
gev.d.lik(xdat, ds, mut, sigma0, xi, theta, eta, log = FALSE) gev.d.lik(xdat, ds, mut, sigma0, xi, theta, eta, tau, log = FALSE)
} }
\arguments{ \arguments{
\item{xdat}{numeric vector containing observations} \item{xdat}{numeric vector containing observations}
\item{ds}{numeric vector containing corresponding durations (1/60 corresponds to 1 minute, 1 corresponds to 1 hour)} \item{ds}{numeric vector containing corresponding durations (1/60 corresponds to 1 minute, 1 corresponds to 1 hour)}
\item{mut, sigma0, xi, theta, eta}{numeric vectors containing corresponding estimates for each of the parameters} \item{mut, sigma0, xi, theta, eta, tau}{numeric vectors containing corresponding estimates for each of the parameters}
\item{log}{Logical; if TRUE, the log likelihood is returned.} \item{log}{Logical; if TRUE, the log likelihood is returned.}
} }
...@@ -29,5 +29,5 @@ fit <- gev.d.fit(train.set$dat,train.set$d,mutl = c(1,2),sigma0l = 1 ...@@ -29,5 +29,5 @@ fit <- gev.d.fit(train.set$dat,train.set$d,mutl = c(1,2),sigma0l = 1
,ydat = as.matrix(train.set[c('cov1','cov2')])) ,ydat = as.matrix(train.set[c('cov1','cov2')]))
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],tau=params[,6],log=TRUE)
} }
...@@ -4,7 +4,7 @@ ...@@ -4,7 +4,7 @@
\alias{pgev.d} \alias{pgev.d}
\title{d-GEV cumulative distribution function} \title{d-GEV cumulative distribution function}
\usage{ \usage{
pgev.d(q, mut, sigma0, xi, theta, eta, d, ...) pgev.d(q, mut, sigma0, xi, theta, eta, tau, d, ...)
} }
\arguments{ \arguments{
\item{q}{vector of quantiles} \item{q}{vector of quantiles}
...@@ -15,6 +15,8 @@ pgev.d(q, mut, sigma0, xi, theta, eta, d, ...) ...@@ -15,6 +15,8 @@ pgev.d(q, mut, sigma0, xi, theta, eta, d, ...)
\item{eta}{numeric value, giving duration exponent (defining slope of the IDF curve)} \item{eta}{numeric value, giving duration exponent (defining slope of the IDF curve)}