Skip to content

GitLab

Commit b798d5de authored by Jana Ulrich's avatar Jana Ulrich
Browse files

initial values now work for link functions

parent 59844a24
Expand all Hide whitespace changes
Inline Side-by-side
Showing with 210 additions and 449 deletions
NAMESPACE
View file @ b798d5de
......@@ -5,10 +5,8 @@ export(IDF.plot)
export(dgev.d)
export(gev.d.diag)
export(gev.d.fit)
export(gev.d.nll)
export(gev.d.params)
export(gev.d.rl)
export(gev.d.rl2)
export(gev.d2stdgumbel)
export(pgev.d)
export(qgev.d)
export(rgev.d)
......@@ -29,5 +27,7 @@ importFrom(graphics,title)
importFrom(ismev,gev.fit)
importFrom(pbapply,pbsapply)
importFrom(stats,lm)
importFrom(stats,make.link)
importFrom(stats,optim)
importFrom(zoo,rollapply)
importFrom(zoo,as.zoo)
importFrom(zoo,rollapplyr)
R/IDF.R
View file @ b798d5de
......@@ -34,7 +34,8 @@
#'
#' @export
#' @importFrom pbapply pbsapply
#' @importFrom zoo rollapply
#' @importFrom zoo rollapplyr
#' @importFrom zoo as.zoo
#'
#' @examples
#' dates <- as.Date("2019-01-01")+0:729
......@@ -68,7 +69,7 @@ IDF.agg <- function(data,ds,na.accept = 0,
# function 1: aggregate over single durations and find annual maxima:
agg.ts <- function(ds){
runmean <- rollapply(data.s[,names[2]],ds/dtime,FUN=sum,fill =NA,align='right')
runmean <- rollapplyr(as.zoo(data.s[,names[2]]),ds/dtime,FUN=sum,fill =NA,align='right')
runmean <- runmean/ds #intensity per hour
subset <- is.element(as.POSIXlt(data.s[,names[1]])$mon,which.mon)
max <- tapply(runmean[subset],(as.POSIXlt(data.s[,names[1]])$year+1900)[subset],
......@@ -133,8 +134,9 @@ IDF.plot <- function(durations,fitparams,probs=c(0.5,0.9,0.99),
if(length(cols)!=length(probs))cols <- rep_len(cols,length.out=length(probs))
## calculate IDF values for given probability and durations
idf.array <- simplify2array(qgev.d(probs,mut=fitparams[1],sigma0=fitparams[2],xi=fitparams[3],
theta=fitparams[4],eta=fitparams[5],d=durations)) # array[probs,durs]
qs <- lapply(durations,qgev.d,p=probs,mut=fitparams[1],sigma0=fitparams[2],xi=fitparams[3],
theta=fitparams[4],eta=fitparams[5])
idf.array <- simplify2array(qs) # array[probs,durs]
if(!add){ #new plot
## initialize plot window
# check if limits were passed
......@@ -166,4 +168,4 @@ IDF.plot <- function(durations,fitparams,probs=c(0.5,0.9,0.99),
legend(x="topright",title = 'p-quantile',legend=probs,
col=cols,lty=lty,lwd=lwd)
}
}
\ No newline at end of file
}
R/d-gev.R
View file @ b798d5de
......@@ -28,7 +28,12 @@
#'
#' @examples
#' x <- seq(4,20,0.1)
#' dens <- dgev.d(q=x,mut=4,sigma0=2,xi=0,theta=0.1,eta=0.1,d=1:4)
#' # calculate probability density for one duration
#' dgev.d(q=x,mut=4,sigma0=2,xi=0,theta=0.1,eta=0.1,d=1)
#'
#' # calculate probability density for different durations
#' ds <- 1:4
#' dens <- lapply(ds,dgev.d,q=x,mut=4,sigma0=2,xi=0,theta=0.1,eta=0.1)
#'
#' plot(x,dens[[1]],type='l',ylim = c(0,0.21),ylab = 'Probability Density')
#' for(i in 2:4){
......@@ -43,15 +48,9 @@ dgev.d <- function(q,mut,sigma0,xi,theta,eta,d,...) {
warning('Some shape parameters are negative,resulting from a negativ theta '
,theta, ' this will prododuce NAs.')}
# if denominator is negative NAs will be returned
densfun <- function(d){
if(d+theta<=0){return(rep(NA,length(q)))}else{
sigma.d <-sigma0/(d+theta)^eta
return(dgev(q,loc=mut*sigma.d,scale=sigma.d,shape=xi,...))}
}
# calculate for each duration in d
durations <- as.list(d)
names(durations) <- d
return(sapply(durations,densfun,simplify = FALSE))
if(d+theta<=0){return(rep(NA,length(q)))}else{
sigma.d <-sigma0/(d+theta)^eta
return(dgev(q,loc=mut*sigma.d,scale=sigma.d,shape=xi,...))}
}
......@@ -83,7 +82,7 @@ dgev.d <- function(q,mut,sigma0,xi,theta,eta,d,...) {
#'
#' @examples
#' x <- seq(4,20,0.1)
#' prob <- pgev.d(q=x,mu=4,sigma=2,xi=0,theta=0.1,eta=0.1,d=c(1,2))
#' prob <- pgev.d(q=x,mu=4,sigma0=2,xi=0,theta=0.1,eta=0.1,d=1)
pgev.d <- function(q,mut,sigma0,xi,theta,eta,d,...) {
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. ',
......@@ -92,15 +91,9 @@ pgev.d <- function(q,mut,sigma0,xi,theta,eta,d,...) {
warning('Some shape parameters are negative,resulting from a negativ theta '
,theta, ' this will prododuce NAs.')}
# if denominator is negative NAs will be returned
pfun <- function(d){
if(d+theta<=0){return(rep(NA,length(q)))}else{
sigma.d <-sigma0/(d+theta)^eta
return(pgev(q,loc=mut*sigma.d,scale=sigma.d,shape=xi,...))}
}
# calculate for each duration in d
durations <- as.list(d)
names(durations) <- d
return(sapply(durations,pfun,simplify = FALSE))
if(d+theta<=0){return(rep(NA,length(q)))}else{
sigma.d <-sigma0/(d+theta)^eta
return(pgev(q,loc=mut*sigma.d,scale=sigma.d,shape=xi,...))}
}
......@@ -132,15 +125,19 @@ pgev.d <- function(q,mut,sigma0,xi,theta,eta,d,...) {
#'
#' @examples
#' p <- c(0.5,0.9,0.99)
#' d <- 2^seq(0,4,0.1)
#' qs <- qgev.d(p=p,mu=4,sigma=2,xi=0,theta=0.1,eta=0.3,d=d)
#' # calulate quantiles for one duration
#' qgev.d(p=p,mu=4,sigma0=2,xi=0,theta=0.1,eta=0.3,d=1)
#'
#' # calculate quantiles for sequence of durations
#' ds <- 2^seq(0,4,0.1)
#' qs <- lapply(ds,qgev.d,p=p,mu=4,sigma0=2,xi=0,theta=0.1,eta=0.3)
#' qs <- simplify2array(qs)
#'
#' plot(d,qs[1,],ylim=c(3,20),type='l',log = 'xy',ylab='Intensity',xlab = 'Duration')
#' plot(ds,qs[1,],ylim=c(3,20),type='l',log = 'xy',ylab='Intensity',xlab = 'Duration')
#' for(i in 2:3){
#' lines(d,qs[i,],lty=i)
#' lines(ds,qs[i,],lty=i)
#' }
#' legend(8,19,title = 'annual frequency \n of exceedance 1/T',
#' legend('topright',title = 'Annual frequency of exceedance',
#' legend = 1-p,lty=1:3,bty = 'n')
qgev.d <- function(p,mut,sigma0,xi,theta,eta,d,...) {
if(any(c(length(mut),length(sigma0),length(xi),length(theta),length(eta))>1)){
......@@ -150,16 +147,10 @@ qgev.d <- function(p,mut,sigma0,xi,theta,eta,d,...) {
warning('Some shape parameters are negative,resulting from a negativ theta '
,theta, ' this will prododuce NAs.')}
# if denominator is negative NAs will be returned
qfun <- function(d){
if(d+theta<=0){return(rep(NA,length(p)))}else{
sigma.d <-sigma0/(d+theta)^eta
return(qgev(p,loc=as.numeric(mut*sigma.d)
,scale=as.numeric(sigma.d),shape=as.numeric(xi),...))}
}
# calculate for each duration in d
durations <- as.list(d)
names(durations) <- d
return(sapply(durations,qfun,simplify = FALSE))
if(d+theta<=0){return(rep(NA,length(p)))}else{
sigma.d <-sigma0/(d+theta)^eta
return(qgev(p,loc=as.numeric(mut*sigma.d)
,scale=as.numeric(sigma.d),shape=as.numeric(xi),...))}
}
......@@ -189,8 +180,13 @@ qgev.d <- function(p,mut,sigma0,xi,theta,eta,d,...) {
#' @importFrom evd rgev
#'
#' @examples
#' samp <- rgev.d(100,4,2,0,0,0.5,c(1,2))
#'
#' # random sample for one duration
#' rgev.d(n=100,mut=4,sigma=2,xi=0,theta=0.1,eta=0.3,d=1)
#'
#' # compare randomn samples for different durations
#' ds <- c(1,4)
#' samp <- lapply(ds,rgev.d,n=100,mut=4,sigma=2,xi=0,theta=0.1,eta=0.3)
#'
#' hist(samp[[1]],breaks = 10,col=rgb(1,0,0,0.5),freq = FALSE
#' ,ylim=c(0,0.3),xlab='x',main = 'd-GEV samples for two different durations')
#' hist(samp[[2]],breaks = 10,add=TRUE,col=rgb(0,0,1,0.5),freq = FALSE)
......@@ -202,15 +198,9 @@ rgev.d <- function(n,mut,sigma0,xi,theta,eta,d) {
warning('Some shape parameters are negative,resulting from a negativ theta '
,theta, ' this will prododuce NAs.')}
# if denominator is negative NAs will be returned
rfun <- function(d){
if(d+theta<=0){return(rep(NA,n))}else{
sigma.d <-sigma0/(d+theta)^eta
return(rgev(n,loc=mut*sigma.d,scale=sigma.d,shape=xi))}
}
# calculate for each duration in d
durations <- as.list(d)
names(durations) <- d
return(sapply(durations,rfun,simplify = FALSE))
if(d+theta<=0){return(rep(NA,n))}else{
sigma.d <-sigma0/(d+theta)^eta
return(rgev(n,loc=mut*sigma.d,scale=sigma.d,shape=xi))}
}
R/gevdfit.R
View file @ b798d5de
This diff is collapsed.
man/dgev.d.Rd
View file @ b798d5de
......@@ -35,7 +35,12 @@ modified location parameter \eqn{\mu_t=\mu/\sigma(d)}.
}
\examples{
x <- seq(4,20,0.1)
dens <- dgev.d(q=x,mut=4,sigma0=2,xi=0,theta=0.1,eta=0.1,d=1:4)
# calculate probability density for one duration
dgev.d(q=x,mut=4,sigma0=2,xi=0,theta=0.1,eta=0.1,d=1)
# calculate probability density for different durations
ds <- 1:4
dens <- lapply(ds,dgev.d,q=x,mut=4,sigma0=2,xi=0,theta=0.1,eta=0.1)
plot(x,dens[[1]],type='l',ylim = c(0,0.21),ylab = 'Probability Density')
for(i in 2:4){
......
man/gev.d.fit.Rd
View file @ b798d5de
......@@ -5,10 +5,11 @@
\title{Maximum-likelihood Fitting of the duration dependent GEV Distribution}
\usage{
gev.d.fit(xdat, ds, 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,
shl = NULL, thetal = NULL, etal = NULL,
mulink = make.link("identity"), siglink = make.link("identity"),
shlink = make.link("identity"), thetalink = make.link("identity"),
etalink = make.link("identity"), muinit = NULL, siginit = NULL,
shinit = NULL, thetainit = NULL, etainit = NULL, show = TRUE,
method = "Nelder-Mead", maxit = 10000, init.vals = NULL, ...)
}
\arguments{
......@@ -26,14 +27,12 @@ covariates for generalized linear modelling of the parameters (or NULL (the defa
if the corresponding parameter is stationary).
Parameters are: modified location, scale_0, shape, duration offset, duration exponent repectively.}
\item{mulink, siglink, shlink, thetalink, etalink}{Inverse link functions for generalized linear
modelling of the parameters.}
\item{mulink, siglink, shlink, thetalink, etalink}{Link functions for generalized linear
modelling of the parameters, created with \code{\link{make.link}}.}
\item{muinit, siginit, shinit, thetainit, etainit}{initial values as numeric of length
equal to total number of parameters
used to model the parameters. If NULL (the default) is given, initial parameters are obtained
internally by fitting the GEV seperately for each duration and applying a linear model to optain the
duration dependency of the location and shape parameter.}
\item{muinit, siginit, shinit, thetainit, etainit}{Initial values as numeric of length
equal to total number of parameters. Alternatively initial values for only the parameter intercepts
can be passed to \code{init.vals}.}
\item{show}{Logical; if TRUE (the default), print details of the fit.}
......@@ -41,7 +40,10 @@ duration dependency of the location and shape parameter.}
\item{maxit}{The maximum number of iterations.}
\item{init.vals}{vector of length 5, giving initial values for parameter intercepts}
\item{init.vals}{vector of length 5, giving initial values for parameter intercepts
used to model the parameters. If NULL (the default) is given, initial parameters are obtained
internally by fitting the GEV seperately for each duration and applying a linear model to optain the
duration dependency of the location and shape parameter.}
\item{...}{Other control parameters for the optimization.}
}
......@@ -50,10 +52,10 @@ A list containing the following components.
A subset of these components are printed after the fit.
If show is TRUE, then assuming that successful convergence is indicated,
the components nllh, mle and se are always printed.
\item{nllh}{single numeric giving the negative log-likelihood value.}
\item{nllh}{single numeric giving the negative log-likelihood value}
\item{mle}{numeric vector giving the MLE's for the modified location, scale_0, shape,
duration offset and duration exponent, resp.}
\item{se}{numeric vector giving the standard errors for the MLE's (in the same order).}
\item{se}{numeric vector giving the standard errors for the MLE's (in the same order)}
\item{trans}{An logical indicator for a non-stationary fit.}
\item{model}{A list with components mul, sigl, shl, thetal and etal.}
\item{link}{A character vector giving inverse link functions.}
......
man/gev.d.init.Rd
View file @ b798d5de
......@@ -4,13 +4,15 @@
\alias{gev.d.init}
\title{get initial values for gev.d.fit}
\usage{
gev.d.init(xdat, ds, thetainit)
gev.d.init(xdat, ds, link)
}
\arguments{
\item{xdat}{vector of maxima for differnt durations}
\item{ds}{vector of durations belonging to maxima in xdat}
\item{link}{list of 5, link functions for parameters, created with \code{\link{make.link}}}
\item{thetainit}{initial parameter for theta}
}
\value{
......
man/gev.d.nll.Rd 0 → 100644
View file @ b798d5de
% Generated by roxygen2: do not edit by hand
% Please edit documentation in R/gevdfit.R
\name{gev.d.nll}
\alias{gev.d.nll}
\title{computes negative log-likelihood of d-gev model}
\usage{
gev.d.nll(xdat, ds, mut, sig0, xi, theta, eta)
}
\arguments{
\item{xdat}{numeric vector containing observations}
\item{ds}{numeric vector containing coresponding durations}
\item{mut, sig0, xi, theta, eta}{numeric vectors containing corresponding mles for each of the parameters}
}
\value{
single value containing negative log likelihood
}
\description{
computes negative log-likelihood of d-gev model
}
\examples{
# compute nll of values not included in fit
train.set <- example[example$d!=2,]
test.set <- example[example$d==2,]
fit <- gev.d.fit(train.set$dat,train.set$d,mul = c(1,2),sigl = 1
,ydat = as.matrix(train.set[c('cov1','cov2')]))
params <- gev.d.params(fit,ydat = as.matrix(test.set[c('cov1','cov2')]))
gev.d.nll(xdat = test.set$dat,ds = test.set$d,mut = params[,1],sig0 = params[,2],xi = params[,3]
,theta = params[,4],eta = params[,5])
}
man/gev.d.rl.Rd deleted 100644 → 0
View file @ 59844a24
% Generated by roxygen2: do not edit by hand
% Please edit documentation in R/gevdfit.R
\name{gev.d.rl}
\alias{gev.d.rl}
\title{Calculate Returnlevel}
\usage{
gev.d.rl(params, p.d)
}
\arguments{
\item{params}{list of parameters mu_tilde, sigma0, xi, theta, eta (single values and/or compatible matrices)
as obtained from \code{\link{gev.d.fit}} or \code{\link{gev.d.params}}}
\item{p.d}{numeric vector of length 2 containing one value the for exeedance probability p = 1-1/RP
and one value for the duration at which to calculate the return level}
}
\value{
one return level value or matrix with return levels (depending on input to params)
unit: e.g. mm/h
}
\description{
calculate Returnlevel for chosen duration and return period
from \code{\link{gev.d.fit}} parameters
}
\examples{
data('example',package = 'IDF')
fit <- gev.d.fit(example$dat,example$d,ydat = as.matrix(example[,c("cov1","cov2")])
,mul = c(1,2),sigl = 1)
### calculate rl on grid:
# covariates values
cov1 <- rep(seq(-1,1,0.1),11)
cov2 <- rep(seq(0,1,0.1),each=21)
# calculate parameters of d-gev on grid, based on output of gev.d.fit
par <- gev.d.params(fit = fit,ydat = cbind(cov1,cov2))
# calculate 100 year (p=0.99) rl for a duration of d=24 hours
rl <- gev.d.rl(params = par,p.d = c(0.99,24))
dim(rl) <- c(21,11)
# rl map:
image(x=seq(-1,1,0.1),y=seq(0,1,0.1),z=rl,xlab = 'cov1', ylab = 'cov2')
}
man/gev.d.rl2.Rd deleted 100644 → 0
View file @ 59844a24
% Generated by roxygen2: do not edit by hand
% Please edit documentation in R/gevdfit.R
\name{gev.d.rl2}
\alias{gev.d.rl2}
\title{Calculate Returnlevel with standart deviation (through delta method)}
\usage{
gev.d.rl2(p, fit, ds, ydat, sd = FALSE)
}
\arguments{
\item{p}{single numeric value containing the non-exceedance probability p=1-1/T}
\item{fit}{\code{\link{gev.d.fit}} object}
\item{ds}{numeric vector of durations for which to calculate the return levels}
\item{ydat}{matrix containing the covariates in the same order as used in \code{\link{gev.d.fit}}}
\item{sd}{logical, if TRUE, standart deviation of the returnlevels obtained with the delta method
are also returned}
}
\value{
either vector or matrix containing the p-quantile or list containing both the p-quantile and its standart deviation
for durations in ds (and different rows for the covariates in ydat)
}
\description{
calculate returnlevel for chosen duration and return period
for \code{\link{gev.d.fit}} objects
}
\examples{
data('example',package = 'IDF')
fit <- gev.d.fit(example$dat,example$d,ydat = as.matrix(example[,c("cov1","cov2")])
,mul = c(1,2),sigl = 1)
# calculate 100 year (p=0.99) rl for a duration of d=24 hours
rl <- gev.d.rl2(p = 0.99,fit = fit,ds = 1:30,ydat = matrix(0.7,1,2),sd=TRUE)
plot(1:30,rl$q,log='xy',type = 'l',xlab = 'Duration',ylab = 'Intensity',ylim = c(5,100))
# add 0.95-confidence intervals
lines(1:30,rl$q-1.96*rl$sd,lty=2)
lines(1:30,rl$q+1.96*rl$sd,lty=2)
}
man/gev.d2stdgumbel.Rd deleted 100644 → 0
View file @ 59844a24
% Generated by roxygen2: do not edit by hand
% Please edit documentation in R/gevdfit.R
\name{gev.d2stdgumbel}
\alias{gev.d2stdgumbel}
\title{Transform data to standart gumbel}
\usage{
gev.d2stdgumbel(xdat, ds, params)
}
\arguments{
\item{xdat}{A vector containing maxima for different durations.}
\item{ds}{A vector of aggregation levels corresponding to the maxima in xdat.}
\item{params}{list of parameters mu_tilde, sigma0, xi, theta, eta
as obtained from \code{\link{gev.d.params}}}
}
\value{
Vector containing transformed data.
}
\description{
Transform data to standart gumbel
}
\examples{
data('example',package = 'IDF')
# fit subset
ind <- sample(1:10, length(example$dat), replace=TRUE)
fit.subs <- gev.d.fit(example$dat[ind!=1],example$d[ind!=1]
,ydat = as.matrix(example[ind!=1,c("cov1","cov2")])
,mul = c(1,2),sigl = 1)
# calculate parameters for unfitted values
par <- gev.d.params(fit = fit.subs
,ydat = as.matrix(example[ind==1,c("cov1","cov2")]))
# transform unfitted values to standart gumbel
sg.data <- gev.d2stdgumbel(xdat = example$dat[ind==1]
,ds = example$d[ind==1],params = par)
# check unfitted values agains standart gumbel
gev.d.diag(data.frame(data=sg.data,ds=example$d[ind==1]),pch=20)
}
man/pgev.d.Rd
View file @ b798d5de
......@@ -35,7 +35,7 @@ modified location parameter \eqn{\mu_t=\mu/\sigma(d)}.
}
\examples{
x <- seq(4,20,0.1)
prob <- pgev.d(q=x,mu=4,sigma=2,xi=0,theta=0.1,eta=0.1,d=c(1,2))
prob <- pgev.d(q=x,mu=4,sigma0=2,xi=0,theta=0.1,eta=0.1,d=1)
}
\references{
Koutsoyannis et al., 1998, doi:10.1016/S0022-1694(98)00097-3
......
man/qgev.d.Rd
View file @ b798d5de
......@@ -35,15 +35,19 @@ modified location parameter \eqn{\mu_t=\mu/\sigma(d)}.
}
\examples{
p <- c(0.5,0.9,0.99)
d <- 2^seq(0,4,0.1)
qs <- qgev.d(p=p,mu=4,sigma=2,xi=0,theta=0.1,eta=0.3,d=d)
# calulate quantiles for one duration
qgev.d(p=p,mu=4,sigma0=2,xi=0,theta=0.1,eta=0.3,d=1)
# calculate quantiles for sequence of durations
ds <- 2^seq(0,4,0.1)
qs <- lapply(ds,qgev.d,p=p,mu=4,sigma0=2,xi=0,theta=0.1,eta=0.3)
qs <- simplify2array(qs)
plot(d,qs[1,],ylim=c(3,20),type='l',log = 'xy',ylab='Intensity',xlab = 'Duration')
plot(ds,qs[1,],ylim=c(3,20),type='l',log = 'xy',ylab='Intensity',xlab = 'Duration')
for(i in 2:3){
lines(d,qs[i,],lty=i)
lines(ds,qs[i,],lty=i)
}
legend(8,19,title = 'annual frequency \\n of exceedance 1/T',
legend('topright',title = 'Annual frequency of exceedance',
legend = 1-p,lty=1:3,bty = 'n')
}
\references{
......
man/rgev.d.Rd
View file @ b798d5de
......@@ -32,7 +32,12 @@ with the duration dependent scale \eqn{\sigma(d)=\sigma_0/(d+\theta)^\eta} and
modified location parameter \eqn{\mu_t=\mu/\sigma(d)}.
}
\examples{
samp <- rgev.d(100,4,2,0,0,0.5,c(1,2))
# random sample for one duration
rgev.d(n=100,mut=4,sigma=2,xi=0,theta=0.1,eta=0.3,d=1)
# compare randomn samples for different durations
ds <- c(1,4)
samp <- lapply(ds,rgev.d,n=100,mut=4,sigma=2,xi=0,theta=0.1,eta=0.3)
hist(samp[[1]],breaks = 10,col=rgb(1,0,0,0.5),freq = FALSE
,ylim=c(0,0.3),xlab='x',main = 'd-GEV samples for two different durations')
......
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