d-gev.R 9.68 KB
Newer Older
1
2
3
4
# This file contains the functions  dgev.d, pgev.d, qgev.d, rgev.d for the duration-dependent-gev.

#### dgev.d() ####

5
#' d-GEV probability density function
6
#'
7
#' @description Probability density function of duration-dependent GEV distribution
8
#' @param q vector of quantiles
9
10
11
12
#' @param mut,sigma0,xi numeric value, giving modified location \eqn{\tilde{\mu}}, scale offset \eqn{\tilde{\sigma_0}} and 
#' shape parameter \eqn{\xi}.
#' @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)
13
#' @param tau numeric value, giving intensity offset \eqn{\tau} (defining flattening of the IDF curve)
Laura Mack's avatar
Laura Mack committed
14
#' @param d positive numeric value, giving duration
15
16
#' @param ... additional parameters passed to \code{\link[evd]{dgev}}
#' 
17
#' @details For details on the d-GEV and the parameter definitions, see \link{IDF-package}.
18
#' 
19
20
#' @return list containing vectors of density values for given quantiles.
#' The first element of the list are the density values for the first given duration etc.
21
22
23
24
25
26
27
28
#' 
#' @seealso \code{\link{pgev.d}}, \code{\link{qgev.d}}, \code{\link{rgev.d}}
#' 
#' @export
#' @importFrom evd dgev 
#'
#' @examples
#' x <- seq(4,20,0.1)
29
#' # calculate probability density for one duration
Felix Fauer's avatar
Felix Fauer committed
30
#' dgev.d(q=x,mut=4,sigma0=2,xi=0,theta=0.1,eta=0.1,d=1)
31
32
33
#' 
#' # calculate probability density for different durations
#' ds <- 1:4
Felix Fauer's avatar
Felix Fauer committed
34
#' dens <- lapply(ds,dgev.d,q=x,mut=4,sigma0=2,xi=0,theta=0.1,eta=0.1)
35
36
37
38
39
#' 
#' plot(x,dens[[1]],type='l',ylim = c(0,0.21),ylab = 'Probability Density')
#' for(i in 2:4){
#'   lines(x,dens[[i]],lty=i)
#' }
40
#' legend('topright',title = 'Duration',legend = 1:4,lty=1:4)
Felix Fauer's avatar
Felix Fauer committed
41
dgev.d <- function(q,mut,sigma0,xi,theta,eta,d,tau=0,...) {
42
43
  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, tau is a vector. ',
44
            'This is not intended and might cause an error.')}
Laura Mack's avatar
Laura Mack committed
45
  if (d<=0) {stop('The duration d has to be positive.')}
46
47
48
49
  if(any(d+theta<=0)){
    warning('Some shape parameters are negative,resulting from a negativ theta '
            ,theta, ' this will prododuce NAs.')}
  # if denominator is negative NAs will be returned
50
  if(d+theta<=0){return(rep(NA,length(q)))}else{
51
    sigma.d <-sigma0/(d+theta)^eta+ tau
52
    return(dgev(q,loc=mut*sigma.d,scale=sigma.d,shape=xi,...))}
53
54
55
56
57
}


#### pgev.d() ####

58
#' d-GEV cumulative distribution function
59
#'
60
#' @description Cumulative probability distribution function of duration-dependent GEV distribution
61
62
63
64
#' @param q vector of quantiles
#' @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 eta numeric value, giving duration exponent (defining slope of the IDF curve)
65
#' @param tau numeric value, giving intensity offset (defining flattening of the IDF curve)
Laura Mack's avatar
Laura Mack committed
66
#' @param d positive numeric value, giving duration
67
68
69
#' @param ... additional parameters passed to \code{\link[evd]{pgev}}
#' 
#' @details The duration dependent GEV distribution is defined after 
Jana Ulrich's avatar
Jana Ulrich committed
70
#' [Koutsoyiannis et al., 1998]:
Laura Mack's avatar
Laura Mack committed
71
#' \deqn{G(x)= \exp[-\left( 1+\xi(x/\sigma(d)-\mu_t) \right)^{-1/\xi}] } 
72
73
74
#' with the duration dependent scale \eqn{\sigma(d)=\sigma_0/(d+\theta)^\eta} and 
#' modified location parameter \eqn{\mu_t=\mu/\sigma(d)}.
#' 
75
76
77
78
#' @details For details on the d-GEV and the parameter definitions, see \link{IDF-package}.
#' 
#' @return list containing vectors of probability values for given quantiles. 
#' The first element of the list are the probability values for the first given duration etc.
79
80
81
82
83
84
85
86
#' 
#' @seealso \code{\link{dgev.d}}, \code{\link{qgev.d}}, \code{\link{rgev.d}}
#' 
#' @export
#' @importFrom evd pgev 
#'
#' @examples
#' x <- seq(4,20,0.1)
Felix Fauer's avatar
Felix Fauer committed
87
88
#' prob <- pgev.d(q=x,mut=4,sigma0=2,xi=0,theta=0.1,eta=0.1,d=1)
pgev.d <- function(q,mut,sigma0,xi,theta,eta,tau=0,d, ...) {
89
90
  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, tau is a vector. ',
91
            'This is not intended and might cause an error.')}
Laura Mack's avatar
Laura Mack committed
92
  if (d<=0) {stop('The duration d has to be positive.')}
93
94
95
96
  if(any(d+theta<=0)){
    warning('Some shape parameters are negative,resulting from a negativ theta '
            ,theta, ' this will prododuce NAs.')}
  # if denominator is negative NAs will be returned
97
  if(d+theta<=0){return(rep(NA,length(q)))}else{
98
    sigma.d <-sigma0/(d+theta)^eta+tau
99
    return(pgev(q,loc=mut*sigma.d,scale=sigma.d,shape=xi,...))}
100
101
102
103
104
}


#### qgev.d() ####

105
#' d-GEV quantile function 
106
#'
107
#' @description Quantile function of duration-dependent GEV distribution (inverse of the cumulative probability distribution function)
108
109
#' @param p vector of probabilities
#' @param mut,sigma0,xi numeric value, giving modified location, modified scale and shape parameter
Felix Fauer's avatar
Felix Fauer committed
110
#' @param theta numeric value, giving duration offset (defining curvature of the IDF curve for short durations)
111
#' @param eta numeric value, giving duration exponent (defining slope of the IDF curve)
Felix Fauer's avatar
Felix Fauer committed
112
#' @param tau numeric value, giving intensity offset (defining flattening of the IDF curve for long durations)
Laura Mack's avatar
Laura Mack committed
113
#' @param d positive numeric value, giving duration
114
115
116
#' @param ... additional parameters passed to \code{\link[evd]{qgev}}
#' 
#' @details The duration dependent GEV distribution is defined after 
Jana Ulrich's avatar
Jana Ulrich committed
117
#' [Koutsoyiannis et al., 1998]:
Laura Mack's avatar
typo    
Laura Mack committed
118
#' \deqn{ G(x)= \exp[-\left( 1+\xi(x/\sigma(d)-\mu_t) \right)^{-1/\xi}] } 
119
120
121
#' with the duration dependent scale \eqn{\sigma(d)=\sigma_0/(d+\theta)^\eta} and 
#' modified location parameter \eqn{\mu_t=\mu/\sigma(d)}.
#' 
122
123
124
125
#' @return list containing vectors of quantile values for given probabilities. 
#' The first element of the list are the q. values for the first given duration etc.
#' 
#' @details For details on the d-GEV and the parameter definitions, see \link{IDF-package}.
126
127
128
129
130
131
132
133
#' 
#' @seealso \code{\link{pgev.d}}, \code{\link{dgev.d}}, \code{\link{rgev.d}}
#' 
#' @export
#' @importFrom evd qgev 
#'
#' @examples
#' p <- c(0.5,0.9,0.99)
134
#' # calulate quantiles for one duration
Felix Fauer's avatar
Felix Fauer committed
135
#' qgev.d(p=p,mut=4,sigma0=2,xi=0,theta=0.1,eta=0.3, d=1)
136
137
138
#' 
#' # calculate quantiles for sequence of durations
#' ds <- 2^seq(0,4,0.1)
Felix Fauer's avatar
Felix Fauer committed
139
#' qs <- lapply(ds,qgev.d,p=p,mut=4,sigma0=2,xi=0,theta=0.1,eta=0.3)
140
141
#' qs <- simplify2array(qs)
#' 
142
#' plot(ds,qs[1,],ylim=c(3,20),type='l',log = 'xy',ylab='Intensity',xlab = 'Duration')
143
#' for(i in 2:3){
144
#'   lines(ds,qs[i,],lty=i)
145
#' }
146
147
#' legend('topright',title = 'p-quantile',
#'        legend = p,lty=1:3,bty = 'n')
Felix Fauer's avatar
Felix Fauer committed
148
qgev.d <- function(p,mut,sigma0,xi,theta,eta,d,tau=0, ...) {
149
150
151
152
153
154
155
156
157
158
  if(any(c(length(mut),length(sigma0),length(xi),length(theta),length(eta), length(tau))>1)){
    message('One of the parameters mut, sigma0, xi, theta, eta is a vector. ',
            'This is not intended and might cause an error.')}
  if (d<=0) {stop('The duration d has to be positive.')}
  if(any(d+theta<=0)){
    warning('Some shape parameters are negative, resulting from a negativ theta '
            ,theta, ' this will prododuce NAs.')}
  # if denominator is negative NAs will be returned
  if(d+theta<=0){return(rep(NA,length(p)))}else{
    #sigma.d <-sigma0/(d+theta)^eta
159
    sigma.d <-sigma0/(d+theta)^eta+tau
160
161
    return(qgev(p,loc=as.numeric(mut*sigma.d)
                ,scale=as.numeric(sigma.d),shape=as.numeric(xi),...))}
162
163
164
165
166
}


#### rgev.d() ####

167
#' Generation of random variables from d-GEV
168
#'
169
170
#' @description Generation of random variables following duration-dependent GEV.
#' @param n number of random variables per duration
171
172
173
#' @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 eta numeric value, giving duration exponent (defining slope of the IDF curve)
174
#' @param tau numeric value, giving intensity offset (defining flattening of the IDF curve)
Laura Mack's avatar
Laura Mack committed
175
#' @param d positive numeric value, giving duration
176
#' 
177
#' @details For details on the d-GEV and the parameter definitions, see \link{IDF-package}
178
#' 
179
180
181
#' @return list containing vectors of random variables.  
#' The first element of the list are the random values for the first given duration etc.
#' Note that the random variables for different durations are nor ordered (contrary to precipitation maxima of different durations).
182
183
184
185
186
187
188
#' 
#' @seealso \code{\link{pgev.d}}, \code{\link{qgev.d}}, \code{\link{dgev.d}}
#' 
#' @export
#' @importFrom evd rgev 
#'
#' @examples
189
#' # random sample for one duration
Felix Fauer's avatar
Felix Fauer committed
190
#' rgev.d(n=100,mut=4,sigma0=2,xi=0,theta=0.1,eta=0.3,d=1)
191
192
193
#' 
#' # compare randomn samples for different durations
#' ds <- c(1,4)
Felix Fauer's avatar
Felix Fauer committed
194
#' samp <- lapply(ds,rgev.d,n=100,mut=4,sigma0=2,xi=0,theta=0.1,eta=0.3)
195
#' 
196
#' hist(samp[[1]],breaks = 10,col=rgb(1,0,0,0.5),freq = FALSE
197
#'      ,ylim=c(0,0.3),xlim=c(3,20),xlab='x',main = 'Random d-GEV samples')
198
#' hist(samp[[2]],breaks = 10,add=TRUE,col=rgb(0,0,1,0.5),freq = FALSE)
199
200
#' legend('topright',fill = c(rgb(1,0,0,0.5),rgb(0,0,1,0.5)),
#' legend = paste('d=',1:2,'h'),title = 'Duration')
Felix Fauer's avatar
Felix Fauer committed
201
rgev.d <- function(n,mut,sigma0,xi,theta,eta,d,tau=0) {
202
203
  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, tau is a vector. ',
204
            'This is not intended and might cause an error.')}
Laura Mack's avatar
Laura Mack committed
205
  if (d<=0) {stop('The duration d has to be positive.')}
206
  if(any(d+theta<=0)){
Laura Mack's avatar
Laura Mack committed
207
    warning('Some shape parameters are negative, resulting from a negative theta '
208
209
            ,theta, ' this will prododuce NAs.')}
  # if denominator is negative NAs will be returned
210
  if(d+theta<=0){return(rep(NA,n))}else{
211
    sigma.d <-sigma0/(d+theta)^eta+tau
212
    return(rgev(n,loc=mut*sigma.d,scale=sigma.d,shape=xi))}
213
214
215
}