Commit 6c1b09f1 authored by Jana Ulrich's avatar Jana Ulrich
Browse files

Didnt have the newest Version of other Functions. Was not sure if I could...

Didnt have the newest Version of other Functions. Was not sure if I could merge branches, so instead just copied everything but gev.d.fit from master branch.
parent 2b117368
##############################################################
##################################################
## IDF package
## Authors: Sarah Joedicke, Carola Detring, Christoph Ritschel
## Update: 15.09.2017
## revise for integration of covariates Sep. 2018
##############################################################
## Update: 15.09.2017
###################################################
##############################################################
## Read data
##############################################################
###############################################
############# Read Data function ##############
###############################################
#' @title Reading precipitation data
#' @description The function \code{IDF.read} reads a file in table format and creates a \code{data.frame} from it
......@@ -109,11 +108,10 @@ IDF.read <- function(file,type){
return(Liste)
}
## End of function IDF.read
# End of function IDF.read
####################################################################################################################
##############################################################
## Accumulation
##############################################################
##### Aggregation ###
#' \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.
......@@ -156,15 +154,16 @@ TS.acc <- function(x,acc.val=2,moving.sum="FALSE") {
## return accumulated time series
return(x.acc)
}
## End of function TS.acc
} # End of function TS.acc
#####################################################################################
##############################################################
## Define duration dependent GEV, d/p/q/rgev.d
##############################################################
#######################
## 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).
#'@description The function \code{dgev.d} is a modified version of the function \code{\link[evd]{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
......@@ -188,7 +187,7 @@ dgev.d <- function(q,mu=0,sigma=1,xi=0,theta=0,eta=1,d=1,log=FALSE) {
#'@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).
#'@description The function \code{qgev.d} is a modified version of the function \code{\link[evd]{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
......@@ -211,7 +210,7 @@ qgev.d <- function(p,mu=0,sigma=1,xi=0,theta=0,eta=1,d=1,lower.tail=TRUE) {
}
#'@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).
#'@description The function \code{rgev.d} is a modified version of the function \code{\link[evd]{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
......@@ -233,10 +232,6 @@ rgev.d <- function(n,mu=0,sigma=1,xi=0,theta=0,eta=1,d=1) {
}
##############################################################
## Define negative log-likelihood
##############################################################
#######################################################################
#' @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
......@@ -255,9 +250,16 @@ rgev.d <- function(n,mu=0,sigma=1,xi=0,theta=0,eta=1,d=1) {
#'@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,DEBUG=F) {
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
if(use.log){
## ensure that critical parameters are positive
sigma <- exp(sigma)
theta <- exp(theta)
eta <- exp(eta)
}
sigma.d <- sigma/((d+theta)^eta)
if(DEBUG) debug.values <- c(mu,sigma,xi,theta,eta)
......@@ -287,8 +289,7 @@ IDF.nll <- function(mu=0,sigma=1,xi=0,theta=0,eta=1,x,d,DEBUG=F) {
return(nll/length(x))
}
## end of function IDF.nll
} # end of function IDF.nll
##################################################################################
### copied gev.fit from ismev to be adapted to gev.d.fit
......@@ -461,12 +462,13 @@ IDF.nll <- function(mu=0,sigma=1,xi=0,theta=0,eta=1,x,d,DEBUG=F) {
class( z) <- "gev.d.fit"
invisible(z)
}
#########################################################################################
######################################################################################################
#' @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.
#' are given. Also the \code{\link[stats]{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
......@@ -477,14 +479,12 @@ IDF.nll <- function(mu=0,sigma=1,xi=0,theta=0,eta=1,x,d,DEBUG=F) {
#' @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
#' @param ... Further arguments to pass to \code{\link[stats]{optim}}.
#' @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,...) {
fit.fun <- function(obs,dur,mu=1,sigma=1,xi=0.5,theta=1,eta=1,use.log=F,DEBUG=F,method="Nelder-Mead",...) {
use.log=use.log
......@@ -510,7 +510,7 @@ fit.fun <- function(obs,dur,mu=1,sigma=1,xi=0.5,theta=1,eta=1,use.log=F,DEBUG=F,
## check initial value of negative log-Likelihood function
nll <- IDF.nll(mu,sigma,xi,theta,eta,x=obs,d=dur,use.log=use.log,DEBUG=DEBUG)
## if initial value is acceptable...
if(!is.infinite(nll)&!is.na(nll)) {
......@@ -519,18 +519,16 @@ fit.fun <- function(obs,dur,mu=1,sigma=1,xi=0.5,theta=1,eta=1,use.log=F,DEBUG=F,
## problem: optimization algrorithm often has difficulities concerning infinite or NA-difference values betweeen iterations
## solution: ignore this error message using functon tryCatch and return NULL if there was an error during optimization
fit <- tryCatch(mle(IDF.nll,start=list(mu=mu,sigma=sigma,xi=xi,theta=theta,eta=eta),fixed=list(x=obs,d=dur,use.log=use.log,DEBUG=DEBUG),
control=list(...),
method=method,upper=upper,lower=lower), error=function(e) NULL)#,
fit <- tryCatch(mle(IDF.nll,start=list(mu=mu,sigma=sigma,xi=xi,theta=theta,eta=eta),
fixed=list(x=obs,d=dur,use.log=use.log,DEBUG=DEBUG),...), error=function(e) NULL)#,
#upper=upper,lower=lower)
}else{
## problem: optimization algrorithm often has difficulities concerning infinite or NA-difference values betweeen iterations
## solution: ignore this error message using functon tryCatch and return NULL if there was an error during optimization
fit <- tryCatch(mle(IDF.nll,start=list(mu=mu,sigma=sigma,xi=xi,theta=theta,eta=eta),fixed=list(x=obs,d=dur,use.log=use.log,DEBUG=DEBUG),
control=list(...),
method=method), error=function(e) NULL)#,
fit <- tryCatch(mle(IDF.nll,start=list(mu=mu,sigma=sigma,xi=xi,theta=theta,eta=eta),
fixed=list(x=obs,d=dur,use.log=use.log,DEBUG=DEBUG),...), error=function(e) NULL)#,
#upper=upper,lower=lower)
......@@ -562,15 +560,12 @@ fit.fun <- function(obs,dur,mu=1,sigma=1,xi=0.5,theta=1,eta=1,use.log=F,DEBUG=F,
return(list("min"=fit.min,"par"=fit.par))
}
## end of function fit.fun
} ## end of function fit.fun
##################################################################################
#' @title Data aggregation for IDF parameter estimation
#' @description The function \code{IDF.agg} aggregates a data.frame of observations \code{data} with temporal inforamtion (at least years) and values of precipitation
#' at a given temporal resoultion at given aggregation levels \code{agg.lev} and yearly maxima of intensity are caluclated for a specific month or the whole year/dataset.
#' @param data a \code{data,frame}, preferably generated by function \code{IDF.read}. It should at least contain a \code{$RR} and \code{$year} element for the
#' @param data a \code{data,frame}, preferably generated by function \code{\link{IDF.read}}. It should at least contain a \code{$RR} and \code{$year} element for the
#' function tow work properly. Also an option to use \code{moving.sum} is given. The function returns a vector of intensities and durations as well as the number of years of data.
#' @param agg.lev a vector of aggregation levels used to fit the IDF curves.
#' @param month \code{integer} value specifying the month to be used for estimating the IDF parameters. Type "all" for all months or if
......@@ -586,7 +581,7 @@ fit.fun <- function(obs,dur,mu=1,sigma=1,xi=0.5,theta=1,eta=1,use.log=F,DEBUG=F,
#' data.agg <- IDF.agg(data)
#' @author Christoph Ritschel \email{christoph.ritschel@@met.fu-berlin.de}
IDF.agg <- function(data,agg.lev=c(2,3,6,12,24,48,72,96),month="all",moving.sum="FALSE",DEBUG=FALSE) {
IDF.agg <- function(data,agg.lev=c(2,3,6,12,24,48,72,96),month="all",moving.sum=FALSE,DEBUG=FALSE) {
RR <- data$RR ## get precipitation time series from data.frame
years <- unique(data$year) # get years from data.frame
......@@ -628,54 +623,77 @@ IDF.agg <- function(data,agg.lev=c(2,3,6,12,24,48,72,96),month="all",moving.sum=
return(list(int.vec=int.vec,durs=durs,n.y=n.y))
} #
###############################################################################
#' @title Estimation of initial values for IDF fitting.
#' @description The function \code{IDF.init} estimates inital values for \code{mu,sigma,xi and eta} assuming \code{theta}
#' equals zero. A generalized extreme value distribution is fitted individually for each year and then the inital values
#' for the duration dependent gev fit are estimated from those by applying a linear regression to the scale parameters of each year.
#' @param xdat a \code{vector} of yearly maxima of intensity sorted by year and aggregatin level
#' @param ds a \code{vector} of durations used to fit the model. Has to have same length and order as \code{int.vec}
#' @param int.vec a \code{vector} of yearly maxima of intensity sorted by year and aggregatin level
#' @param durs a \code{vector} of durations used to fit the model. Has to have same length and order as \code{int.vec}
#' @param n.y \code{integer} value specifying the number of years of data.
#' @param method \code{character} defining the method to be used in \code{\link[stats]{optim}}, preferences are: "Nelder-Mead", "BFGS", "L-BFGS-B"
#' @param ... Other contral parameters for the optimization. These are passed to components of the control argument of \code{\link[stats]{optim}}.
#' @return $mu initial estimation of location parameter
#' @return $sigma initial estimation of scale parameter
#' @return $xi inital estimation of shape parameter
#' @return $eta intial estimation of slope parameter for sigma-power law.
#' @examples
#' RR <- rgamma(10*30*24,shape=1)
#' year <- sort(rep(1:(10),30*24))
#' data <- data.frame(RR,year)
#' data.agg <- IDF.agg(data,agg.lev=c(2,6,12,24))
#' pars.init <- IDF.init(data.agg$int.vec,data.agg$durs,data.agg$n.y)
#' @author Christoph Ritschel \email{christoph.ritschel@@met.fu-berlin.de}
#' @author Henning Rust \email{henning.rust@fu-berlin.de}
IDF.init <- function(xdat,ds) {
IDF.init <- function(int.vec,durs,n.y,method="Nelder-Mead",...) {
## Fit a generalized extreme value distribution to the maximum intensities of each year for a single
## aggregation level and write the estimated parameters in an array for further analyisis.
pars <- simplify2array(tapply(xdat,ds,function(xdat) gev.fit(xdat,show=FALSE)$mle,simplify=TRUE))
if(anyNA(pars)){
cat("Warning: optimization did not converge in some cases and no parameters were estimated.\n")
mu <- sigma <- xi <- eta <- NA
}else{
#############################################################
### Derive starting parameters for duration-dependent GEV ###
#############################################################
## Fit a linear model to the individual sigmas for individual aggregation times in a log-log environment
## The slope coefficient is an estimate for the slope in the duration-dependent GEV, namely parameter eta
## The intersection is an estimation of the starting parameter sigma
## Parameter mu is estimated as mean value of individual mus divided by indiviudal sigmas
## The initial value for xi will be the mean of all individual xi, since it is approximately independent of duration
formel <- lm(log(pars[2,]) ~ log(as.numeric(dimnames(pars)[[2]])))
sigma <- as.numeric(exp(formel$coefficients[1]))
mu <- mean(pars[1,]/pars[2,])
eta <- as.numeric(-formel$coefficients[2])
xi <- max(0,mean(pars[3,],na.rm=T))
}
## Fit a generalized extreme value distribution to the maximum intensities of each year for a single
## aggregation level and write the estimated parameters in an array for further analyisis.
d.all <- unique(durs)
ints.all <- matrix(int.vec,nrow=n.y) ## sort intensities in a matrix, rows are years, columns are aggregation levels
pars <- array(NA,dim=c(3,length(d.all)))
## In case of NA values the optimization fails, therefore years with NA values need to be removed.
ints.all <- matrix(ints.all[rowSums(!is.na(ints.all)) == length(d.all)],ncol=length(d.all))
if(nrow(ints.all)<3) {
cat("Warning: optimization did not converge and no parameters were estimated. Time Series contains less than 3 years of valid data. \n")
mu=NA
sigma=NA
xi=NA
eta=NA
}else{
## loop over all aggregation levels
for(d in 1:length(d.all)) {
#fit <- fit.fun.emp(obs=ints.all[,d],mu=mu,sigma=sigma,xi=xi,use.log=use.log,
# DEBUG=DEBUG,method=method,upper=upper,lower=lower)
fit <- gev.fit(xdat=ints.all[,d],method=method,show=F,...)
pars[,d] <- fit$mle
} ## end loop over aggregation levels
#############################################################
### Derive starting parameters for duration-dependent GEV ###
#############################################################
## Fit a linear model to the individual sigmas for individual aggregation times in a log-log environment
## The slope coefficient is an estimate for the slope in the duration-dependent GEV, namely parameter eta
## The intersection is an estimation of the starting parameter sigma
## Parameter mu is estimated as mean value of individual mus divided by indiviudal sigmas
## The initial value for xi will be the mean of all individual xi, since it is approximately independent of duration
formel <- lm(log(pars[2,]) ~ log(d.all))
sigma <- as.numeric(exp(formel$coefficients[1]))
mu <- mean(pars[1,]/pars[2,])
eta <- as.numeric(-formel$coefficients[2])
xi <- max(0,mean(pars[3,],na.rm=T))
}
return(list("mu"=mu,"sigma"=sigma,"xi"=xi,"eta"=eta))
}
} # EOF
#################################################################################
......@@ -690,6 +708,7 @@ IDF.init <- function(xdat,ds) {
#' IDF curves.
#' @param data a \code{data,frame}, preferably generated by function \code{IDF.read}. It should at least contain a \code{$RR} and \code{$year} element for the
#' function tow work properly.
#' @param ... Arguments to be passed to function \code{\link[graphics]{plot}}, such as \code{graphical parameters} (see \code{\link[graphics]{par}}).
#' @param agg.lev a vector of aggregation levels used to fit the IDF curves.
#' @param month \code{integer} value specifying the month to be used for estimating the IDF parameters. Type "all" for all months or if
#' the whole time series should be fitted.
......@@ -702,7 +721,7 @@ IDF.init <- function(xdat,ds) {
#' @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 method \code{character} defining the method to be used in \code{\link[stats]{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
#' @param plot \code{logical} option of creating a plot of IDF curves with estimated parameters.
......@@ -722,9 +741,9 @@ IDF.init <- function(xdat,ds) {
#' pars <- fit$par
#' @author Christoph Ritschel \email{christoph.ritschel@@met.fu-berlin.de}
IDF <- function(data,agg.lev=c(2,3,6,12,24,48,72,96),month="all",moving.sum="FALSE",mu.init=NA,sigma.init=NA,xi.init=NA,theta.init=0,eta.init=NA,
IDF <- function(data,...,agg.lev=c(2,3,6,12,24,48,72,96),month="all",moving.sum=FALSE,mu.init=NA,sigma.init=NA,xi.init=NA,theta.init=0,eta.init=NA,
use.log=FALSE,DEBUG=FALSE,method="Nelder-Mead",upper=Inf,lower=-Inf,plot=FALSE,
probs=c(0.5,0.9,0.99),cols=c(rgb(1,0,0,1),rgb(0,1,0,1),rgb(0,0,1,1)),station.name="Berlin",data.name="obs",...) {
probs=c(0.5,0.9,0.99),cols=rainbow(length(probs)),station.name="Berlin",data.name="obs") {
#########################################################################
### Calculate extreme values for each year and each aggregation level ###
......@@ -752,7 +771,7 @@ IDF <- function(data,agg.lev=c(2,3,6,12,24,48,72,96),month="all",moving.sum="FAL
######################################################
if(!is.na(mu.init) | !is.na(sigma.init) | !is.na(xi.init) | !is.na(eta.init)) {
fit <- fit.fun(obs=int.vec,dur=durs,mu=mu.init,sigma=sigma.init,xi=xi.init,theta=theta.init,eta=eta.init,use.log=use.log,
DEBUG=DEBUG,method=method,upper=upper,lower=lower,...)
DEBUG=DEBUG,method=method,upper=upper,lower=lower)
}else {
cat("Warning: Optimization not carried out due to invalid initial values. \n")
fit.min <- NA
......@@ -762,8 +781,9 @@ IDF <- function(data,agg.lev=c(2,3,6,12,24,48,72,96),month="all",moving.sum="FAL
######################################################
if(plot&& !is.na(fit$min)) {
d.all <- unique(durs)
ds <- sort(rep(d.all,length(int.vec)/length(d.all)))
IDF.plot(pars=fit$par,probs,st.name=station.name,dt.name=data.name,ints=int.vec,ds=durs)
IDF.plot(fit$par,...,probs=probs,st.name=station.name,dt.name=data.name,ints=int.vec,ds=durs)
}
......@@ -777,9 +797,10 @@ IDF <- function(data,agg.lev=c(2,3,6,12,24,48,72,96),month="all",moving.sum="FAL
return(list("ints"=int.vec,"durs"=durs,"min"=fit$min,"par"=fit$par))
}
## End of function IDF.fit
} ## End of function IDF.fit
######################################################################################################################
####
#' @title Fitting IDF model parameters to annual maximum intensity time series
#' @description The function \code{IDF.short} fits the IDF model parameters \code{mu,sigma,xi,eta,theta}
......@@ -790,6 +811,7 @@ IDF <- function(data,agg.lev=c(2,3,6,12,24,48,72,96),month="all",moving.sum="FAL
#' @param ints.vec a \code{vector} of yearly maxima of intensity sorted by year and aggregatin level
#' @param durs a vector of aggregation levels used to fit the IDF curves. One value for each year. Has to have same lenght as \code{int.vec}
#' @param n.y \code{integer} value specifying the number of years of data
#' @param ... Arguments to be passed to function \code{\link[graphics]{plot}}, such as \code{graphical parameters} (see \code{\link[graphics]{par}}).
#' @param mu.init initial estimation of location parameter, default is NA. Initial value estimated by fitting individual gev parameters
#' @param sigma.init initial estimation of scale parameter,default is NA. Initial value estimated by fitting individual gev parameters
#' @param xi.init inital estimation of shape parameter, default is NA. Initial value estimated by fitting individual gev parameters
......@@ -798,7 +820,7 @@ IDF <- function(data,agg.lev=c(2,3,6,12,24,48,72,96),month="all",moving.sum="FAL
#' @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 method \code{character} defining the method to be used in \code{\link[stats]{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
#' @param plot \code{logical} option of creating a plot of IDF curves with estimated parameters.
......@@ -824,7 +846,7 @@ IDF <- function(data,agg.lev=c(2,3,6,12,24,48,72,96),month="all",moving.sum="FAL
IDF.short <- function(int.vec,durs,n.y,mu.init=NA,sigma.init=NA,xi.init=NA,theta.init=0,eta.init=NA,
use.log=FALSE,DEBUG=FALSE,method="Nelder-Mead",upper=Inf,lower=-Inf,plot=FALSE,
probs=c(0.5,0.9,0.99),cols=c(rgb(1,0,0,1),rgb(0,1,0,1),rgb(0,0,1,1)),
probs=c(0.5,0.9,0.99),cols=rainbow(length(probs)),
station.name="Station",data.name="obs",...) {
###################################################################################
......@@ -842,9 +864,9 @@ IDF.short <- function(int.vec,durs,n.y,mu.init=NA,sigma.init=NA,xi.init=NA,theta
######################################################
### Estimate parameters for duration-dependent GEV ###
######################################################
if(!is.na(mu.init) | !is.na(sigma.init) | !is.na(xi.init) | !is.na(eta.init)) {
if(!is.na(mu.init) | !is.na(sigma.init) | !is.na(xi.init) | !is.na(eta.init)) {
fit <- fit.fun(obs=int.vec,dur=durs,mu=mu.init,sigma=sigma.init,xi=xi.init,theta=theta.init,eta=eta.init,use.log=use.log,
DEBUG=DEBUG,method=method,upper=upper,lower=lower,...)
DEBUG=DEBUG,method=method,upper=upper,lower=lower)
}else {
cat("Warning: Optimization not carried out due to invalid initial values. \n")
fit.min <- NA
......@@ -856,7 +878,7 @@ IDF.short <- function(int.vec,durs,n.y,mu.init=NA,sigma.init=NA,xi.init=NA,theta
if(plot&& !is.na(fit$min)) {
d.all <- unique(durs)
ds <- sort(rep(d.all,length(int.vec)/length(d.all)))
IDF.plot(pars=fit$par,probs,st.name=station.name,dt.name=data.name,ints=int.vec,ds=durs)
IDF.plot(fit$par,...,probs=probs,st.name=station.name,dt.name=data.name,ints=int.vec,ds=durs)
}
......@@ -880,6 +902,7 @@ IDF.short <- function(int.vec,durs,n.y,mu.init=NA,sigma.init=NA,xi.init=NA,theta
#' several probability levels \code{probs} at given durations \code{dur}. The colors of the curves can be defined with
#' parameter \code{cols} (need to have same length as \code{probs}). The \code{station.name} will be printed in the legend.
#' @param pars a vector of IDF model parameters mu,sigma,xi,eta,theta
#' @param ... Arguments to be passed to methods, such as \code{graphical parameters} (see \code{\link[graphics]{par}}).
#' @param probs a vector of probabilities for which the IDF curves are calculated
#' @param dur a vector of durations at which the IDF curves are calculated
#' @param cols a vector of colors for the seperate IDF curves, needs same length as \code{probs}
......@@ -896,10 +919,10 @@ IDF.short <- function(int.vec,durs,n.y,mu.init=NA,sigma.init=NA,xi.init=NA,theta
#' IDF.plot(pars=param,st.name="example",dt.name="rgamma")
#' @author Christoph Ritschel \email{christoph.ritschel@@met.fu-berlin.de}
IDF.plot <- function(pars,probs=c(0.5,0.9,0.99),
IDF.plot <- function(pars,...,probs=c(0.5,0.9,0.99),
dur=c(0.5,1,2,3,6,12,24,48,72,96),
cols=rainbow(length(probs)),lty=1,
st.name="Station",dt.name="obs",ints=NA,ds=NA,ylim=NA,add=FALSE,...) {
st.name="Station",dt.name="obs",ints=NA,ds=NA,ylim=c(NA,NA),add=FALSE) {
## initialize array for IDF values at different durations and for different probabilities
idf.array <- array(NA,dim=c(length(dur),length(probs)))
......@@ -911,46 +934,25 @@ IDF.plot <- function(pars,probs=c(0.5,0.9,0.99),
idf.array[,i] <- qgev.d(probs[i],mu=pars[1],sigma=pars[2],xi=pars[3],theta=pars[4],eta=pars[5],d=dur)
} ## end of loop over probs
if(!add){
## initiialize plot window with limits of IDF values
y.range <- ifelse(is.na(ylim), c(min(idf.array[,1],na.rm=T),max(idf.array[,3],na.rm=T)),ylim)
plot(NA,axes=F,xlim=c(min(dur,na.rm=T),max(dur,na.rm=T)),ylim=y.range,xlab="duration [h]",ylab="intensity [mm/h]",log="xy",...)
axis(1,at=dur,labels=dur)
axis(2)
points(ds,ints,pch=16,col=rgb(0,0,0,0.5))
## loop over probabilities
## plot IDF curve
legend.text.2 <- "quantile"
## plot legend
legend(x="topright",legend=c(st.name,dt.name,paste(probs,legend.text.2,sep=" ")),
col=c(1,rgb(0,0,0,0.5),cols),lty=c(NA,NA,rep(1,length(cols))),pch=c(NA,16,rep(NA,length(cols))))
}
for(i in 1:length(probs))
lines(dur,idf.array[,i],col=cols[i],lwd=1.5,lty=lty)
if(!add){
## initiialize plot window with limits of IDF values
y.range <- ifelse(is.na(ylim), c(min(idf.array[,1],na.rm=T),max(idf.array[,length(probs)],na.rm=T)),ylim)
plot(NA,...,axes=F,xlim=c(min(dur,na.rm=T),max(dur,na.rm=T)),ylim=y.range,xlab="duration [h]",ylab="intensity [mm/h]",log="xy")
axis(1,at=dur,labels=dur)
axis(2)
points(ds,ints,pch=16,col=rgb(0,0,0,0.5))
## loop over probabilities
## plot IDF curve
legend.text.2 <- "quantile"
## plot legend
legend(x="topright",legend=c(st.name,dt.name,paste(probs,legend.text.2,sep=" ")),
col=c(1,rgb(0,0,0,0.5),cols),lty=c(NA,NA,rep(1,length(cols))),pch=c(NA,16,rep(NA,length(cols))))
}
for(i in 1:length(probs))
lines(dur,idf.array[,i],col=cols[i],lwd=1.5,lty=lty)
} ## end of function IDF.plot
###################################################################################
ydat = NULL
mul = NULL
sigl = NULL
shl = NULL
thetal = NULL
etal = NULL
mulink = identity
siglink = identity
shlink = identity
thetalink = identity
etalink = identity
muinit = NULL
siginit = NULL
shinit = NULL
thetainit = NULL
etainit = NULL
show = TRUE
method = "Nelder-Mead"
maxit = 10000
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