Commit d4ef04ee authored by Christoph Ritschel's avatar Christoph Ritschel
Browse files

*Paketversion1.1*

parent d98f1874
Package: IDF
Type: Package
Title: Estimating intensity-duration-frequency parameters
Version: 1.01
Date: 2016-01-28
Author: Christoph Ritschel, Sarah Joedicke
Title: Estimation and Plotting of IDF Curves
Version: 1.1
Date: 2017-09-29
Author: Christoph Ritschel, Carola Detring, Sarah Joedicke
Maintainer: Christoph Ritschel <christoph.ritschel@met.fu-berlin.de>
Description: The package 'IDF' reads precipitation data from webwerdis and Berlin Stadtmessnetz.
It aggregates these data and calculates yearly intensity maxima. From those maxima intensity-duration-frequency
parameters are estimated. IDF curves based on these estimated parameters can be plotted.
Depends: stats4, evd
Description: Intensity-duration-frequency (IDF) curves are a widely used analysis-tool in hydrology to assess extreme values of precipitation [e.g. Mailhot et al., 2007, <doi:10.1016/j.jhydrol.2007.09.019>]. The package 'IDF' provides a function to read precipitation data from German weather service (DWD) 'webwerdis' <http://www.dwd.de/EN/ourservices/webwerdis/webwerdis.html> files and Berlin station data from 'Stadtmessnetz' <http://www.geo.fu-berlin.de/en/met/service/stadtmessnetz/index.html> files, and additionally IDF parameters can be estimated also from a given data.frame containing a precipitation time series. The data is aggregated to given levels yearly intensity maxima are calculated either for the whole year or given months. From these intensity maxima IDF parameters are estimated on the basis of a duration-dependent generalised extreme value distribution [Koutsoyannis et al., 1998, <doi:10.1016/S0022-1694(98)00097-3>]. IDF curves based on these estimated parameters can be plotted.
Depends:
stats4,
evd,
ismev
License: GPL (>=2)
RoxygenNote: 5.0.1
exportPattern(IDF.read,IDF.fit,IDF.plot)
export(IDF.read,IDF.fit,IDF.plot,TS.acc)
import(ismev,
stats4,
evd)
importFrom("grDevices", "rgb")
importFrom("graphics", "axis", "legend", "points")
importFrom("stats", "filter", "lm")
importFrom("utils", "read.csv2", "read.table", "str")
# Generated by roxygen2: do not edit by hand
# Import all packages listed as Imports or Depends
import(
stats4,
evd
)
<
##################################################
## IDF package
## Authors: Sarah Joedicke, Christoph Ritschel
## Update: 29.01.2016
## Authors: Sarah Joedicke, Carola Detring, Christoph Ritschel
## Update: 15.09.2017
###################################################
###############################################
############# Read Data function ##############
###############################################
##' @title Reading precipitation data
##' @description The function \code{IDF.read} reads a file in table format and creates a \code{list} from it, writing the data frame as first entry
##' and adds some attributes (station information, aggregation time, data source). The only data values used are:
##' date, precipitation
##' The \code{data.frame} will have the following format:
##' | year | mon | day | hour | min | RR |
##' |------+-----+-----+------+-----+----+
##' | | | | | | |
##' @usage IDF.read(file, type)
##' @param file a \code{character string} naming the file from which the data is to be read.
##' @param type a \code{character string} defining the type of data to be read: either "stadtmessnetz" or "webwerdis", depending on if the data comes from the Stadtmessnetz Berlin
##' or WebWerdis. If type = "webwerdis", the data will be read, then sorted, formatted and missing lines added,
##' while if type = "stadtmessnetz", the data will just be read and formatted.
##' Both source types have a different layout in the original file.
##' @return Liste a \code{list} containing a \code{data.frame} of date and time information and precipitation values for each time step
##' @details This function is designed to prepare a data file for aggregation by functions \code{IDF.agg} or \code{IDF.magg}.
##' The time given in the data is the end time, so the precipitation was measured up to that time.
##' @note The first entry of the list will be a \code{data.frame} named 't1'.
##' @seealso read.table, IDF.agg, IDF.magg
##' @author Sarah Joedicke \email{sarah.joedicke@@fu-berlin.de}
#' @title Reading precipitation data
#' @description The function \code{IDF.read} reads a file in table format and creates a \code{data.frame} from it
#' and adds some attributes (station information, aggregation time, data source). The only data values used are:
#' date, precipitation
#' The \code{data.frame} will have the following format:
#' | year | mon | day | hour | min | RR |
#' |------+-----+-----+------+-----+----+
#' | | | | | | |
#' @usage IDF.read(file, type)
#' @param file a \code{character string} naming the file from which the data is to be read.
#' @param type a \code{character string} defining the type of data to be read: either "stadtmessnetz" or "webwerdis", depending on if the data comes from the Stadtmessnetz Berlin
#' or WebWerdis. If type = "webwerdis", the data will be read, then sorted, formatted and missing lines added,
#' while if type = "stadtmessnetz", the data will just be read and formatted.
#' Both source types have a different layout in the original file.
#' @return Liste a \code{data.frame} of date and time information and precipitation values for each time step
#' @details This function is designed to prepare a data file for doing an estimation on IDF parameters in function \code{IDF.fit}.
#' The time given in the data is the end time, so the precipitation was measured up to that time.
#' @seealso read.table, IDF.fit
#' @author Sarah Joedicke \email{sarah.joedicke@@fu-berlin.de}
#' @author Christoph Ritschel \email{christoph.ritschel@@met.fu-berlin.de}
IDF.read <- function(file,type){
if(type != "stadtmessnetz" && type != "webwerdis") {
......@@ -40,7 +39,7 @@ IDF.read <- function(file,type){
if (type == "stadtmessnetz") {
Tab_MN <- read.csv2(file) #STADTMESSNETZ
new_time <- strptime(Tab_MN$Zeitstempel,format="%d.%m.%Y %H:%M") #STADTMESSNETZ Datumsvektor
new_time <- strptime(Tab_MN$Zeitstempel,format="%d.%m.%Y %H:%M") #STADTMESSNETZ date vector
}
# Da die Stadtmessnetzdaten (bisher) konstistent aussehen, wird auf das Erstellen einer neuen Tabelle mit sicher allen
......@@ -51,20 +50,20 @@ IDF.read <- function(file,type){
Tab <- read.table(file,header=TRUE,sep=";") #WEBWERDIS
Tab_kurz <- Tab[,c("Date","precipitation")]
## Tabelle in gewuenschtes Format sortieren; nur WEBWERDIS, da Stadtmessnetz sortiert und keine fehlenden Zeilen
## Sort table in output format
time <- strptime(Tab_kurz$Date,format="%Y-%m-%d T %H:%M:%S")
Tab_sort <- Tab_kurz[order(as.character(time)),]
time_sort <- strptime(Tab_sort$Date,format="%Y-%m-%d T %H:%M:%S")
Tab_sort$Date <- as.character(time_sort)
## Fehlende Zeilen in Daten hinzufuegen, Niederschlagswert = NA; nur WEBWERDIS
## If dates are missing, add lines containing NA preicipitation measurments for these time steps.
h_diff <- as.numeric(difftime(format(time_sort[length(time_sort)],"%Y-%m-%d T %H:%M:%S") ,
format(time_sort[1],"%Y-%m-%d T %H:%M:%S"),units="hours")) #h_diff ist die gesamte Zeitdifferenz
format(time_sort[1],"%Y-%m-%d T %H:%M:%S"),units="hours")) #h_diff is the difference in time steps
new_time <- seq(time_sort[1], length = h_diff+1, by = "hour")
new_tab <- data.frame(Date=as.character(new_time), precipitation=NA) # Tabelle mit NAs und richtigen Zeiten vordefinieren
new_tab <- data.frame(Date=as.character(new_time), precipitation=NA) # predefine table with NAs and every time steps
Tab_na <- (merge(Tab_sort, new_tab, "Date", all.y=TRUE))[,1:2]
} #Daten sortieren und in gewuenschte Form bringen
}
new_timect <- as.POSIXct(new_time)
......@@ -77,12 +76,13 @@ IDF.read <- function(file,type){
if (type == "webwerdis") Tab_end <- data.frame(J,M,d,h,m,Tab_na$precipitation.x) #WEBWERDIS
if (type == "stadtmessnetz") Tab_end <- data.frame(J,M,d,h,m,Tab_MN[,2]) #STADTMESSNETZ
## Tabellen-Attribute benennen:
## Name table attributes:
colnames(Tab_end) <- c("year","mon","day","hour","min","RR")
attr(Tab_end,"accumulation time (min)") <- as.numeric(difftime(new_timect[2],new_timect[1], units="mins"))
Liste <- list(t1=Tab_end)
# Liste <- list(t1=Tab_end)
Liste <- Tab_end
if (type == "webwerdis"){
# WEBWERDIS:
attr(Liste,"StationName") <- as.character(Tab$Stationname[1])
......@@ -108,223 +108,147 @@ IDF.read <- function(file,type){
return(Liste)
}
# End of function IDF.read
####################################################################################################################
##### Aggregation ###
############################################
########## Moving aggragation #############
############################################
##'@title Moving aggregation of a time series
##'@description The function \code{IDF.magg} aggregates the precipitation values of a \code{data.frame} in a \code{list} by moving an aggregation
##'window over the time steps. The list should be created by funcntion \code{IDF.read} in a previous step.
##'@usage IDF.magg(x, ks=c(2,4,8,16,32,64,128))
##'@param x a \code{list} containing a \code{data.frame} named 't1', preferably generated by function \code{IDF.read}.
##'@param ks a numeric \code{vector} of aggregatation levels. Each value represents an aggregation level and will create
##'a new \code{data.frame} in the output list.
##'@details This function is designed to aggregate data which has been read by function \code{IDF.read}.
##'The time given in the data is the end time, so the precipitation was added up to that time.
##'Every line of data will be used in this process. In the end, the data frames containing aggregated precipitation
##'will be exactly as long as the original one. In the end, a new \code{list} containing several \code{data.frames} will be created.
##'The first one will be the initial \code{data.frame}, now named agg0. The other data frames will be named 'agg' + the aggregation level number from ks.
##'@return Liste a \code{list} containing \code{data.frames} of original and aggregated precipitation time series with time information and
##'precipitation values.
##'@seealso IDF.read, IDF.agg.
##' @author Sarah Joedicke \email{sarah.joedicke@@fu-berlin.de}
IDF.magg <- function(x,ks=c(2,4,8,16,32,64,128)){
NS <- x$t1$Niederschlagssumme #Niederschlagsvektor
Liste <- vector("list", length(ks)+1) # Leere Liste vordefinieren, +1 fuer Ursprungstabelle
Liste[[1]] <- x$t1
j <- 2
for (y in ks){
agg <- vector(length=nrow(x$t1))
print(paste("Aggregation ueber Aggregationsstufe", y, "laeuft"))
for (i in y:nrow(x$t1)){
agg[i] <- sum(NS[(i-(y-1)):i])
}
agg[1:(y-1)] <- NA
Tab_agg <- x$t1[,1:5]
Tab_agg[,6] <- agg
colnames(Tab_agg) <- c("Jahr","Monat","Tag","Stunde","Minute","Niederschlagssumme (mm)")
attr(Tab_agg,"Aggregationszeit (min)") <- attr(x$t1, "Aggregationszeit (min)")*y
Liste[[j]] <- Tab_agg
j <- j+1
}
attributes(Liste) <- attributes(x)
names(Liste) <- paste("agg", c(0,ks), sep="")
print("Vorgang abgeschlossen")
return(Liste)
}
############################################
########## diskret aggregation #############
############################################
##'@title Aggregation of a time series for seperated time steps
##'@description The function \code{IDF.agg} aggregates the precipitation values of a \code{data.frame} in a \code{list} in groups.
##'The list should be created by function \code{IDF.read} in a previous step.
##'@usage IDF.agg(x, ks=c(2,4,8,16,32,64,128))
##'@param x a \code{list} containing a \code{data.frame} named 't1', preferably generated by function \code{IDF.read}.
##'@param ks a numeric \code{vector} of aggregatation levels. Each value represents an aggregation level and will create
##'a new \code{data.frame} in the output list.
##'@return Liste a \code{list} containing \code{data.frames} of original and aggregated precipitation time series with time information and
##'precipitation values.
##'@details This function is designed to aggregate data which has been read by function \code{IDF.read}.
##'The time given in the data is the end time, so the precipitation was added up to that time.
##'In this process the data will be summed up in groups, depending on the chosen aggregation levels. In the end,
##'the aggregated \code{data.frames} will be much shorter than the original one.
##'In the end, a new \code{list} containing several \code{data.frames} will be created. The first one will be the
##'initial data frame, now named agg0. The other data frames will be named 'agg' + the aggregation level number from ks.
##'@seealso read.data, agg.f().
##' @author Sarah Joedicke \email{sarah.joedicke@@fu-berlin.de}
IDF.agg <- function(x,ks=c(2,4,8,16,32,64,128)){
NS <- x$t1$RR
date <- paste(as.character(x$t1$year), as.character(x$t1$mon), as.character(x$t1$day),
as.character(x$t1$hour), as.character(x$t1$min))
Liste <- vector("list", length(ks)+1) # Leere Liste vordefinieren, +1 fuer Ursprungstabelle
Liste[[1]] <- x$t1
j <- 2
for (y in ks){
agg_NS <- apply(matrix(data=NS[1:((length(NS)%/%y)*y)], nrow=(length(NS)%/%y), ncol=y, byrow=TRUE), 1, sum)
# %/% dividiert und rundet ganzzahlig ab
# Variable if NS%%y > 0 sind das die unbenutzen Resttage, die letzten "" Zeitschritte werden ignoriert
if (length(NS)%%y > 0) cat("Note: the last", length(NS)%%y,"timesteps are ignored \n")
agg_date <- matrix(data=date[1:((length(NS)%/%y)*y)], nrow=(length(NS)%/%y), ncol=y, byrow=TRUE)[,y]
new_time2 <- (strptime(agg_date, format="%Y %m %d %H %M"))
J <- as.numeric(format(new_time2,'%Y'))
M <- as.numeric(format(new_time2,'%m'))
d <- as.numeric(format(new_time2,'%d'))
h <- as.numeric(format(new_time2,'%H'))
m <- as.numeric(format(new_time2,'%M'))
Tab_agg <- data.frame(J,M,d,h,m,agg_NS)
colnames(Tab_agg) <- c("year","mon","day","hour","min","RR")
attr(Tab_agg,"accumulation time (min)") <- as.numeric(difftime(new_time2[2],new_time2[1]),units="mins")
Liste[[j]] <- Tab_agg
j <- j+1
#' \code{TS.acc} accumulates a given time series \code{x} at a given accumulation level \code{acc.val}. Minimum value
#' for acc.val is 2 [unit time]. Option for using moving sum is given.
#' @title Accumulation of a time series
#' @param x \code{vector} of a time series
#' @param acc.val \code{value} specifying the accumulation level, minimum value is 2
#' @param moving.sum \code{logical} 'TRUE' means moving sum will be applied
#' @return x.acc \code{TS.acc} returns a \code{vector} of an accumulated time series
#' @usage TS.acc(x,acc.val,moving.sum="FALSE")
#' @examples
#' TS <- rgamma(n=1000,shape=1)
#' acc.2 <- TS.acc(TS,acc.val=2)
#' \donttest{
#' acc.24 <- TS.acc(TS,acc.val=24,moving.sum=TRUE)
#' }
#' @author Christoph Ritschel \email{christoph.ritschel@@met.fu-berlin.de}
#' @author Carola Detring \email{carola.detring@@met.fu-berlin.de}
TS.acc <- function(x,acc.val=2,moving.sum="FALSE") {
## check for input value of acc.val
if(acc.val<1) cat(paste("Warning: accumulation value acc.val too small for accumulation of the time series \n"))
if(moving.sum){
x.acc <- as.numeric(filter(x,filter=rep(1,acc.val),method="convolution",sides=1))
}else{
l.new <- length(x)%/%acc.val ## calculate new length of accumulated time series
l.rest <- length(x)%%acc.val ## calculate values left over
if(l.rest==0) {
x.acc <- apply(matrix(x,nrow=l.new,byrow=T),1,sum)
}else{
x.acc <- apply(matrix(x[1:(length(x)-l.rest)],nrow=l.new,byrow=T),1,sum)
#cat(paste("Warning: ",l.rest,"time steps left and not used for accumulation \n"))
}
attributes(Liste) <- attributes(x)
names(Liste) <- paste("agg", c(0,ks), sep="")
return(Liste)
} # End of function agg
################################################################################################################
##'@title Yearly precipitation maximum
##'@description The function \code{IDF.max} generates an \code{vector} of yearly precipitation maxima from a \code{data.frame}
##'containing time information and precipitation data at a single aggregation level generated by functions \code{IDF.agg} or \code{IDF.magg}. An option
##'for selection a single month or search the whole year exists.
##'@param df a \code{data.frame} containing date and time information and the corresponding precipitation values at a single aggregation level.
##'@param month a \code{value} defining the month to use. A single number between 1 and 12 selects the month. The \code{character string} "all"
##'can be used if the maximum for all months is needed.
##'@return RR.max a \code{vector} containing the precipitation maxima for each year of a defined month or all months.
##' @author Christoph Ritschel \email{christoph.ritschel@@met.fu-berlin.de}
IDF.max <- function(df,month=1) {
years <- unique(df$year)
y.len <- length(years)
RR.max <- array(NA,dim=c(y.len))
for(y in 1:y.len) {
if(month=="all") { index <- which(df$year==years[y])
}else index <- which(df$year==years[y] & df$mon==month)
if(sum(index)>0) RR.max[y] <- max(df$RR[index],na.rm=T)
}
RR.max <- RR.max[!is.na(RR.max)]
return(RR.max)
}
## return accumulated time series
return(x.acc)
} # End of function TS.acc
#####################################################################################
#######################
## Fitting Functions ##
#######################
##'@title Density function of modified generalized extreme value distribution
##'@description The function \code{dgev.d} is a modified version of the function \code{dgev} for different durations \code{d} developed by Koutsoyiannis et al. (1998).
##'@param q Vector of quantiles
##'@param mu location value
##'@param sigma scale value
##'@param xi shape value
##'@param theta value defining the curvature of the IDF
##'@param eta value defining the slope of the IDF
##'@param d vector of durations
##'@param log \code{logical} option to use logarithmic parameter values, default=FALSE
##'@seealso \code{\link[evd]{dgev}}
##'@return dgev.d gives the density function
##' @author Christoph Ritschel \email{christoph.ritschel@@met.fu-berlin.de}
#'@title Density function of modified generalized extreme value distribution
#'@description The function \code{dgev.d} is a modified version of the function \code{dgev} for different durations \code{d} developed by Koutsoyiannis et al. (1998).
#'@param q Vector of quantiles
#'@param mu location value
#'@param sigma scale value
#'@param xi shape value
#'@param theta value defining the curvature of the IDF
#'@param eta value defining the slope of the IDF
#'@param d vector of durations
#'@param log \code{logical} option to use logarithmic parameter values, default=FALSE
#'@seealso \code{\link[evd]{dgev}}
#'@return dgev.d gives the density function
#' @author Christoph Ritschel \email{christoph.ritschel@@met.fu-berlin.de}
dgev.d <- function(q,mu=0,sigma=1,xi=0,theta=0,eta=1,d=1,log=FALSE) {
sigma.d <- sigma/(d+theta)^eta
dgev(q,loc=mu*sigma.d,scale=sigma.d,shape=xi,log=log)
##problem if sigma.d is NaN (d+theta) negative and eta smaller than 1 --> cant calculate root of negative value
sigma.d[which(is.nan(sigma.d))] <- Inf
dens <- dgev(q,loc=mu*sigma.d,scale=sigma.d,shape=xi,log=log)
dens[which(is.nan(dens))] <- NA
return(dens)
}
##'@title Quantile function of modified generalized extreme value distribution
##'@description The function \code{qgev.d} is a modified version of the function \code{qgev} for different durations \code{d} developed by Koutsoyiannis et al. (1998).
##'@param p Vector of probabilities
##'@param mu location value
##'@param sigma scale value
##'@param xi shape value
##'@param theta value defining the curvature of the IDF
##'@param eta value defining the slope of the IDF
##'@param d vector of durations
##'@param lower.tail \code{logical} if TRUE (default), probabilities are P[X <= x], otherwise, P[X > x]
##'@seealso \code{\link[evd]{qgev}}
##'@return qgev.d gives the quantile function
##' @author Christoph Ritschel \email{christoph.ritschel@@met.fu-berlin.de}
#'@title Quantile function of modified generalized extreme value distribution
#'@description The function \code{qgev.d} is a modified version of the function \code{qgev} for different durations \code{d} developed by Koutsoyiannis et al. (1998).
#'@param p Vector of probabilities
#'@param mu location value
#'@param sigma scale value
#'@param xi shape value
#'@param theta value defining the curvature of the IDF
#'@param eta value defining the slope of the IDF
#'@param d vector of durations
#'@param lower.tail \code{logical} if TRUE (default), probabilities are P[X <= x], otherwise, P[X > x]
#'@seealso \code{\link[evd]{qgev}}
#'@return qgev.d gives the quantile function
#' @author Christoph Ritschel \email{christoph.ritschel@@met.fu-berlin.de}
qgev.d <- function(p,mu=0,sigma=1,xi=0,theta=0,eta=1,d=1,lower.tail=TRUE) {
sigma.d <- sigma/(d+theta)^eta
qgev(p,loc=mu*sigma.d,scale=sigma.d,shape=xi,lower.tail=lower.tail)
##problem if sigma.d is NaN (d+theta) negative and eta smaller than 1 --> cant calculate root of negative value
sigma.d[which(is.nan(sigma.d))] <- Inf
quant <- qgev(p,loc=mu*sigma.d,scale=sigma.d,shape=xi,lower.tail=lower.tail)
quant[is.infinite(quant)] <- NA
return(quant)
}
##'@title Random generation for the modified generalized extreme value distribution
##'@description The function \code{rgev.d} is a modified version of the function \code{rgev} for different durations \code{d} developed by Koutsoyiannis et al. (1998).
##'@param n Number of observations
##'@param mu location value
##'@param sigma scale value
##'@param xi shape value
##'@param theta value defining the curvature of the IDF
##'@param eta value defining the slope of the IDF
##'@param d vector of durations
##'@seealso \code{\link[evd]{rgev}}
##'@return rgev.d generates random derivates
##' @author Christoph Ritschel \email{christoph.ritschel@@met.fu-berlin.de}
#'@title Random generation for the modified generalized extreme value distribution
#'@description The function \code{rgev.d} is a modified version of the function \code{rgev} for different durations \code{d} developed by Koutsoyiannis et al. (1998).
#'@param n Number of observations
#'@param mu location value
#'@param sigma scale value
#'@param xi shape value
#'@param theta value defining the curvature of the IDF
#'@param eta value defining the slope of the IDF
#'@param d vector of durations
#'@seealso \code{\link[evd]{rgev}}
#'@return rgev.d generates random derivates
#' @author Christoph Ritschel \email{christoph.ritschel@@met.fu-berlin.de}
rgev.d <- function(n,mu=0,sigma=1,xi=0,theta=0,eta=1,d=1) {
## gumbel
sigma.d <- sigma/(d+theta)^eta
rgev(n, loc=mu*sigma.d,scale=sigma.d,shape=xi)
}
##problem if sigma.d is NaN (d+theta) negative and eta smaller than 1 --> cant calculate root of negative value
sigma.d[which(is.nan(sigma.d))] <- Inf
x <- rgev(n, loc=mu*sigma.d,scale=sigma.d,shape=xi)
x[which(is.nan(x))] <- NA
return(x)
}
#######################################################################
##' @title Negativ log-likelihood of modified GEV
##' @description The function \code{IDF.nll} calculates the negative log-likelihood for a given set of model parameters
##' \code{mu,sigma,xi,theta,eta}, given observations \code{x} and given durations \code{d}. Options for the usage of
##' logartihmic values \code{use.log} and a debugging function \code{DEBUG} are available.
##'@param mu location value
##'@param sigma scale value
##'@param xi shape value
##'@param theta value defining the curvature of the IDF
##'@param eta value defining the slope of the IDF
##'@param x vector of observations at different durations d
##'@param d vector of durations
##'@param use.log \code{logical} value for usage of logarithmic values, default is \code{FALSE}
##'@param DEBUG \code{logical} value for usage of debugging, if \code{TRUE} the input parameters and the value of negative
##'log-likelihood are printed on console.
##'@return retruns weightes negative log-likelihood by number of observatons uesd
##' @author Christoph Ritschel \email{christoph.ritschel@@met.fu-berlin.de}
#' @title Negativ log-likelihood of modified GEV
#' @description The function \code{IDF.nll} calculates the negative log-likelihood for a given set of model parameters
#' \code{mu,sigma,xi,theta,eta}, given observations \code{x} and given durations \code{d}. Options for the usage of
#' logartihmic values \code{use.log} and a debugging function \code{DEBUG} are available.
#'@param mu location value
#'@param sigma scale value
#'@param xi shape value
#'@param theta value defining the curvature of the IDF
#'@param eta value defining the slope of the IDF
#'@param x vector of observations at different durations d
#'@param d vector of durations
#'@param use.log \code{logical} value for usage of logarithmic values, default is \code{FALSE}
#'@param DEBUG \code{logical} value for usage of debugging, if \code{TRUE} the input parameters and the value of negative
#'log-likelihood are printed on console.
#'@return retruns weightes negative log-likelihood by number of observatons uesd
#' @author Christoph Ritschel \email{christoph.ritschel@@met.fu-berlin.de}
IDF.nll <- function(mu=0,sigma=1,xi=0,theta=0,eta=1,x,d,use.log=F,DEBUG=F) {
## mu is the mu~ from Koutsoyiannis
......@@ -336,23 +260,25 @@ IDF.nll <- function(mu=0,sigma=1,xi=0,theta=0,eta=1,x,d,use.log=F,DEBUG=F) {
eta <- exp(eta)
}
sigma.d <- sigma/((d+theta)^eta)
sigma.d <- sigma/((d+theta)^eta)
if(DEBUG) debug.values <- c(mu,sigma,xi,theta,eta)
if(sum(is.nan(sigma.d))==0) {
## Weibull und Frechet
if(xi!=0){
C <- 1 + xi * (x/sigma.d - mu )
nll <- switch((sum(C<0)>0)+1,
nll <- switch((sum(C<0,na.rm=T)>0)+1,
sum(log(sigma.d),na.rm=T)+(1+1/xi)*sum(log(C),na.rm=T)+sum((C)^(-1/xi),na.rm=T),
Inf)
NA)
# + penalty*(sum(C[C<0]^2))
## Gumbel
}else if(xi==0){# & sigma<1 & eta<1)
Y <- x/sigma.d-mu
nll <- -(-sum(log(sigma.d))-sum((Y))-sum(exp(-Y)))
nll <- -(-sum(log(sigma.d),na.rm=T)-sum((Y),na.rm=T)-sum(exp(-Y),na.rm=T))
}
}else{ nll <- NA}
if(DEBUG){
cat(debug.values,nll,"\n")
options(digits.secs=6)
......@@ -360,42 +286,46 @@ IDF.nll <- function(mu=0,sigma=1,xi=0,theta=0,eta=1,x,d,use.log=F,DEBUG=F) {
## write(debug.values,file="optim.log",append=TRUE,ncolumns=length(debug.values))
## cat(debug.values,nll,sum(A<0),"\n")
}
return(nll/length(x))
} # end of function IDF.nll
######################################################################################################
##' @title Fitting function to optimize IDF model parameters
##' @description The function \code{fit.fun} fits IDF model parameters \code{mu,sigma,xi,theta,eta} to a set of given observations \code{obs},
##' typically a series of yearly maxima at different durations \code{d}. Options for using logarithmic parameter values and debugging
##' are given. Also the \code{optim} parameters \code{method} and \code{upper,lower} can be defined.
##' @param obs vector of yearly intensity maxima at different durations. Order: Y1D1, Y2D1,...,YnD1,Y1D2,...YnD2,Y1D3,...,YnDk
##' @param dur vector of durations with same length as \code{obs}. Order: n x D1, n x D2, ... n x Dk
##'@param mu location value
##'@param sigma scale value
##'@param xi shape value
##'@param theta value defining the curvature of the IDF
##'@param eta value defining the slope of the IDF
##'@param use.log \code{logical} value for usage of logarithmic values, default is \code{FALSE}
##'@param DEBUG \code{logical} value for usage of debugging, if \code{TRUE} the input parameters and the value of negative
##'log-likelihood are printed on console for each iteration during optimization.
##' @param method \code{character} defining the method to be used in \code{optim}, preferences are: "Nelder-Mead", "BFGS", "L-BFGS-B"e
##' @param lower \code{vector} specifying the lower boundary of parameters for "L-BFGS-B" method
##' @param upper \code{vector} specifying the upper boundary of parameters for "L-BFGS-B" method
##' @return $min value of negative log-likelihood at optimization minimum
##' @return $par vector of IDF parameters at optimization minimum
##' @author Christoph Ritschel \email{christoph.ritschel@@met.fu-berlin.de}
#' @title Fitting function to optimize IDF model parameters
#' @description The function \code{fit.fun} fits IDF model parameters \code{mu,sigma,xi,theta,eta} to a set of given observations \code{obs},
#' typically a series of yearly maxima at different durations \code{d}. Options for using logarithmic parameter values and debugging
#' are given. Also the \code{optim} parameters \code{method} and \code{upper,lower} can be defined.
#' @param obs vector of yearly intensity maxima at different durations. Order: Y1D1, Y2D1,...,YnD1,Y1D2,...YnD2,Y1D3,...,YnDk
#' @param dur vector of durations with same length as \code{obs}. Order: n x D1, n x D2, ... n x Dk
#' @param mu location value
#' @param sigma scale value
#' @param xi shape value
#' @param theta value defining the curvature of the IDF
#' @param eta value defining the slope of the IDF
#' @param use.log \code{logical} value for usage of logarithmic values, default is \code{FALSE}
#' @param DEBUG \code{logical} value for usage of debugging, if \code{TRUE} the input parameters and the value of negative
#' log-likelihood are printed on console for each iteration during optimization.
#' @param method \code{character} defining the method to be used in \code{optim}, preferences are: "Nelder-Mead", "BFGS", "L-BFGS-B"e
#' @param lower \code{vector} specifying the lower boundary of parameters for "L-BFGS-B" method
#' @param upper \code{vector} specifying the upper boundary of parameters for "L-BFGS-B" method
#' @return $min value of negative log-likelihood at optimization minimum
#' @return $par vector of IDF parameters at optimization minimum
#' @author Christoph Ritschel \email{christoph.ritschel@@met.fu-berlin.de}
fit.fun <- function(obs,dur,mu=1,sigma=1,xi=0.5,theta=1,eta=1,use.log=F,DEBUG=F,method="Nelder-Mead",upper=Inf,lower=-Inf) {
use.log=use.log
if(use.log) {
if(sigma<=0){sigma <- 1E-10}
if(theta<=0){theta <- 1E-10}
if(eta<=0){eta <- 1E-10}
sigma <- log(sigma)
theta <- log(theta)
eta <- log(eta)
if(method=="L-BFG-S") {
if(method=="L-BFGS-B") {
upper[2] <- log(upper[2])
upper[4] <- log(upper[4])
upper[5] <- log(upper[5])
......@@ -464,61 +394,158 @@ fit.fun <- function(obs,dur,mu=1,sigma=1,xi=0.5,theta=1,eta=1,use.log=F,DEBUG=F,
} ## end of function fit.fun
##################################################################################
##' @title Fitting IDF model parameters to observations at different durations
##' @description The function \code{IDF.fit} fits the IDF model parameters \code{mu,sigma,xi,eta,theta}
##' to a list object of observations \code{data} in a data.frame with temporal inforamtion and values of precipitation
##' at a given temporal resoultion. This precipitation time series gets aggregated at given aggregation levels
##' \code{agg.lev} and yearly maxima of intensity are caluclated for a specific month or the whole year.
##' The starting values of the IDF model parameters can be determined by the user as well as specific options to use
##' during optimization. Logartihmic transformation, debugging, the optimization method, and an option to plot the
##' IDF curves.
##'@param data a \code{list} containing a \code{data.frame} named 't1', preferably generated by function \code{IDF.read}.
##'@param agg.lev a vector of aggregation levels used to fit the IDF curves.
##'@param month a value specifying the month to be used for estimating the IDF parameters. Type "all" for all months.
##'@param mu location value
##'@param sigma scale value
##'@param xi shape value
##'@param theta value defining the curvature of the IDF
##'@param eta value defining the slope of the IDF
##'@param use.log \code{logical} value for usage of logarithmic values, default is \code{FALSE}
##'@param DEBUG \code{logical} value for usage of debugging, if \code{TRUE} the input parameters and the value of negative