d-gev.R 8.92 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
# This file contains the functions  dgev.d, pgev.d, qgev.d, rgev.d for the duration-dependent-gev.

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

#' Density function of duration dependent GEV distribution
#'
#' @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)
#' @param d numeric value, giving duration
#' @param ... additional parameters passed to \code{\link[evd]{dgev}}
#' 
#' @details The duration dependent GEV distribution is defined after 
#' [Koutsoyannis et al., 1998]:
#' \deqn{ G(x)= exp[-{ 1+\xi(x/\sigma(d)-\mu_t) }^(-1/\xi)] } 
#' with the duration dependent scale \eqn{\sigma(d)=\sigma_0/(d+\theta)^\eta} and 
#' modified location parameter \eqn{\mu_t=\mu/\sigma(d)}.
#' 
#' @return list containing vectors of density values for given quantiles 
#' the first element of the list are the dens. values for the first given duration and so on
#' 
#' @seealso \code{\link{pgev.d}}, \code{\link{qgev.d}}, \code{\link{rgev.d}}
#' @references Koutsoyannis et al., 1998, doi:10.1016/S0022-1694(98)00097-3
#' 
#' @export
#' @importFrom evd dgev 
#'
#' @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)
#' 
#' 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)
#' }
#' legend('topright',title = 'duration',legend = 1:4,lty=1:4)
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.')}
  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
  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
52
53
54
  durations <- as.list(d)
  names(durations) <- d
  return(sapply(durations,densfun,simplify = FALSE))
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
}


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

#' Distribution function of duration dependent GEV distribution
#'
#' @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)
#' @param d numeric value, giving duration
#' @param ... additional parameters passed to \code{\link[evd]{pgev}}
#' 
#' @details The duration dependent GEV distribution is defined after 
#' [Koutsoyannis et al., 1998]:
#' \deqn{ G(x)= exp[-{ 1+\xi(x/\sigma(d)-\mu_t) }^(-1/\xi)] } 
#' with the duration dependent scale \eqn{\sigma(d)=\sigma_0/(d+\theta)^\eta} and 
#' modified location parameter \eqn{\mu_t=\mu/\sigma(d)}.
#' 
#' @return list containing vectors of probability values for given quantiles 
#' the first element of the list are the prob. values for the first given duration and so on
#' 
#' @seealso \code{\link{dgev.d}}, \code{\link{qgev.d}}, \code{\link{rgev.d}}
#' @references Koutsoyannis et al., 1998, doi:10.1016/S0022-1694(98)00097-3
#' 
#' @export
#' @importFrom evd pgev 
#'
#' @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))
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.')}
  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
  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
101
102
103
  durations <- as.list(d)
  names(durations) <- d
  return(sapply(durations,pfun,simplify = FALSE))
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
}


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

#' Quantile function of duration dependent GEV distribution
#'
#' @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)
#' @param d numeric value, giving duration
#' @param ... additional parameters passed to \code{\link[evd]{qgev}}
#' 
#' @details The duration dependent GEV distribution is defined after 
#' [Koutsoyannis et al., 1998]:
#' \deqn{ G(x)= exp[-{ 1+\xi(x/\sigma(d)-\mu_t) }^(-1/\xi)] } 
#' with the duration dependent scale \eqn{\sigma(d)=\sigma_0/(d+\theta)^\eta} and 
#' modified location parameter \eqn{\mu_t=\mu/\sigma(d)}.
#' 
#' @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 and so on
#' 
#' @seealso \code{\link{pgev.d}}, \code{\link{dgev.d}}, \code{\link{rgev.d}}
#' @references Koutsoyannis et al., 1998, doi:10.1016/S0022-1694(98)00097-3
#' 
#' @export
#' @importFrom evd qgev 
#'
#' @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)
#' qs <- simplify2array(qs)
#' 
#' plot(d,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)
#' }
#' legend(8,19,title = 'annual frequency \n of exceedance 1/T',
#'        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)){
    message('One of the parameters mut, sigma0, xi, theta, eta is a vector. ',
            'This is not intended and might cause an error.')}
  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
  qfun <- function(d){
    if(d+theta<=0){return(rep(NA,length(p)))}else{
      sigma.d <-sigma0/(d+theta)^eta
156
157
      return(qgev(p,loc=as.numeric(mut*sigma.d)
                  ,scale=as.numeric(sigma.d),shape=as.numeric(xi),...))}
158
159
  }
  # calculate for each duration in d
160
161
162
  durations <- as.list(d)
  names(durations) <- d
  return(sapply(durations,qfun,simplify = FALSE))
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
}


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

#' random generation of variables from duration dependent GEV distribution
#'
#' @param n number of observations.
#' @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)
#' @param d numeric value, giving duration
#' 
#' @details The duration dependent GEV distribution is defined after 
#' [Koutsoyannis et al., 1998]:
#' \deqn{ G(x)= exp[-{ 1+\xi(x/\sigma(d)-\mu_t) }^(-1/\xi)] } 
#' with the duration dependent scale \eqn{\sigma(d)=\sigma_0/(d+\theta)^\eta} and 
#' modified location parameter \eqn{\mu_t=\mu/\sigma(d)}.
#' 
#' @return list containing vectors of random values  
#' the first element of the list are the random values for the first given duration and so on
#' 
#' @seealso \code{\link{pgev.d}}, \code{\link{qgev.d}}, \code{\link{dgev.d}}
#' @references Koutsoyannis et al., 1998, doi:10.1016/S0022-1694(98)00097-3
#' 
#' @export
#' @importFrom evd rgev 
#'
#' @examples
#' samp <- rgev.d(100,4,2,0,0,0.5,c(1,2))
#'
#' 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)
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.')}
  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
  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
211
212
213
  durations <- as.list(d)
  names(durations) <- d
  return(sapply(durations,rfun,simplify = FALSE))
214
215
216
}