d-gev.R 10 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)
Laura Mack's avatar
Laura Mack committed
13
#' @param d positive numeric value, giving duration
14
15
#' @param ... additional parameters passed to \code{\link[evd]{dgev}}
#' 
16
#' @details For details on the d-GEV and the parameter definitions, see \link{IDF-package}.
17
#' 
18
19
#' @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.
20
21
22
23
24
25
26
27
#' 
#' @seealso \code{\link{pgev.d}}, \code{\link{qgev.d}}, \code{\link{rgev.d}}
#' 
#' @export
#' @importFrom evd dgev 
#'
#' @examples
#' x <- seq(4,20,0.1)
28
29
30
31
32
33
#' # 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)
34
35
36
37
38
#' 
#' 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)
#' }
39
#' legend('topright',title = 'Duration',legend = 1:4,lty=1:4)
40
41
42
43
dgev.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. ',
            'This is not intended and might cause an error.')}
Laura Mack's avatar
Laura Mack committed
44
  if (d<=0) {stop('The duration d has to be positive.')}
45
46
47
48
  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
49
50
51
  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,...))}
52
53
54
55
56
}


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

57
#' d-GEV cumulative distribution function
58
#'
59
#' @description Cumulative probability distribution function of duration-dependent GEV distribution
60
61
62
63
#' @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)
Laura Mack's avatar
Laura Mack committed
64
#' @param d positive numeric value, giving duration
65
66
67
#' @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
68
#' [Koutsoyiannis et al., 1998]:
Laura Mack's avatar
Laura Mack committed
69
#' \deqn{G(x)= \exp[-\left( 1+\xi(x/\sigma(d)-\mu_t) \right)^{-1/\xi}] } 
70
71
72
#' with the duration dependent scale \eqn{\sigma(d)=\sigma_0/(d+\theta)^\eta} and 
#' modified location parameter \eqn{\mu_t=\mu/\sigma(d)}.
#' 
73
74
75
76
#' @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.
77
78
79
80
81
82
83
84
#' 
#' @seealso \code{\link{dgev.d}}, \code{\link{qgev.d}}, \code{\link{rgev.d}}
#' 
#' @export
#' @importFrom evd pgev 
#'
#' @examples
#' x <- seq(4,20,0.1)
85
#' prob <- pgev.d(q=x,mut=4,sigma0=2,xi=0,theta=0.1,eta=0.1,d=1)
86
87
88
89
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. ',
            'This is not intended and might cause an error.')}
Laura Mack's avatar
Laura Mack committed
90
  if (d<=0) {stop('The duration d has to be positive.')}
91
92
93
94
  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
95
96
97
  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,...))}
98
99
100
101
102
}


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

103
#' d-GEV quantile function 
104
#'
105
#' @description Quantile function of duration-dependent GEV distribution (inverse of the cumulative probability distribution function)
106
107
108
109
#' @param p vector of probabilities
#' @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)
Laura Mack's avatar
Laura Mack committed
110
#' @param d positive numeric value, giving duration
111
112
113
#' @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
114
#' [Koutsoyiannis et al., 1998]:
Laura Mack's avatar
typo    
Laura Mack committed
115
#' \deqn{ G(x)= \exp[-\left( 1+\xi(x/\sigma(d)-\mu_t) \right)^{-1/\xi}] } 
116
117
118
#' with the duration dependent scale \eqn{\sigma(d)=\sigma_0/(d+\theta)^\eta} and 
#' modified location parameter \eqn{\mu_t=\mu/\sigma(d)}.
#' 
119
120
121
122
#' @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}.
123
124
125
126
127
128
129
130
#' 
#' @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)
131
#' # calulate quantiles for one duration
132
#' qgev.d(p=p,mut=4,sigma0=2,xi=0,theta=0.1,eta=0.3,d=1)
133
134
135
#' 
#' # calculate quantiles for sequence of durations
#' ds <- 2^seq(0,4,0.1)
136
#' qs <- lapply(ds,qgev.d,p=p,mut=4,sigma0=2,xi=0,theta=0.1,eta=0.3)
137
138
#' qs <- simplify2array(qs)
#' 
139
#' plot(ds,qs[1,],ylim=c(3,20),type='l',log = 'xy',ylab='Intensity',xlab = 'Duration')
140
#' for(i in 2:3){
141
#'   lines(ds,qs[i,],lty=i)
142
#' }
143
144
#' legend('topright',title = 'p-quantile',
#'        legend = p,lty=1:3,bty = 'n')
145
qgev.d.old <- function(p,mut,sigma0,xi,theta,eta,d,...) { # cannot deal with tau
146
147
148
  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. ',
            'This is not intended and might cause an error.')}
Laura Mack's avatar
Laura Mack committed
149
  if (d<=0) {stop('The duration d has to be positive.')}
150
  if(any(d+theta<=0)){
151
    warning('Some shape parameters are negative, resulting from a negativ theta '
152
153
            ,theta, ' this will prododuce NAs.')}
  # if denominator is negative NAs will be returned
154
155
156
157
  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),...))}
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
} # can deal with tau!
qgev.d <- function(p,mut,sigma0,xi,theta,eta,d, tau=NULL, ...) {
  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
    ifelse(!is.null(tau), sigma.d <-sigma0/(d+theta)^eta+tau, sigma.d <-sigma0/(d+theta)^eta)
    return(qgev(p,loc=as.numeric(mut*sigma.d)
                ,scale=as.numeric(sigma.d),shape=as.numeric(xi),...))}
173
174
175
176
177
}


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

178
#' Generation of random variables from d-GEV
179
#'
180
181
#' @description Generation of random variables following duration-dependent GEV.
#' @param n number of random variables per duration
182
183
184
#' @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)
Laura Mack's avatar
Laura Mack committed
185
#' @param d positive numeric value, giving duration
186
#' 
187
#' @details For details on the d-GEV and the parameter definitions, see \link{IDF-package}
188
#' 
189
190
191
#' @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).
192
193
194
195
196
197
198
#' 
#' @seealso \code{\link{pgev.d}}, \code{\link{qgev.d}}, \code{\link{dgev.d}}
#' 
#' @export
#' @importFrom evd rgev 
#'
#' @examples
199
#' # random sample for one duration
200
#' rgev.d(n=100,mut=4,sigma0=2,xi=0,theta=0.1,eta=0.3,d=1)
201
202
203
#' 
#' # compare randomn samples for different durations
#' ds <- c(1,4)
204
#' samp <- lapply(ds,rgev.d,n=100,mut=4,sigma0=2,xi=0,theta=0.1,eta=0.3)
205
#' 
206
#' hist(samp[[1]],breaks = 10,col=rgb(1,0,0,0.5),freq = FALSE
207
#'      ,ylim=c(0,0.3),xlim=c(3,20),xlab='x',main = 'Random d-GEV samples')
208
#' hist(samp[[2]],breaks = 10,add=TRUE,col=rgb(0,0,1,0.5),freq = FALSE)
209
210
#' legend('topright',fill = c(rgb(1,0,0,0.5),rgb(0,0,1,0.5)),
#' legend = paste('d=',1:2,'h'),title = 'Duration')
211
212
213
214
rgev.d <- function(n,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. ',
            'This is not intended and might cause an error.')}
Laura Mack's avatar
Laura Mack committed
215
  if (d<=0) {stop('The duration d has to be positive.')}
216
  if(any(d+theta<=0)){
Laura Mack's avatar
Laura Mack committed
217
    warning('Some shape parameters are negative, resulting from a negative theta '
218
219
            ,theta, ' this will prododuce NAs.')}
  # if denominator is negative NAs will be returned
220
221
222
  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))}
223
224
225
}