gevdfit.R 22.2 KB
Newer Older
1
2
3
# This file contains the functions:
# - gev.d.fit, gev.d.init for fitting
# - gev.d.diag for diagnostic plots
4
# - gev.d.params for calculation of parameters
5
6
7
8
9
10
11
# and the documentation of the example data

#### gev.d.fit ####

#' @title Maximum-likelihood Fitting of the duration dependent GEV Distribution
#' @description Modified \code{\link[ismev]{gev.fit}} function for Maximum-likelihood fitting 
#' for the duration dependent generalized extreme 
12
#' value distribution, following Koutsoyiannis et al. (1998), including generalized linear 
13
14
15
#' modelling of each parameter.
#' @param xdat A vector containing maxima for different durations. 
#' This can be obtained from \code{\link{IDF.agg}}.
16
17
#' @param ds A vector of aggregation levels corresponding to the maxima in xdat. 
#' 1/60 corresponds to 1 minute, 1 corresponds to 1 hour.
18
19
20
21
22
23
24
#' @param ydat A matrix of covariates for generalized linear modelling of the parameters 
#' (or NULL (the default) for stationary fitting). The number of rows should be the same as the 
#' length of xdat.
#' @param  mul,sigl,shl,thetal,etal Numeric vectors of integers, giving the columns of ydat that contain
#'  covariates for generalized linear modelling of the parameters (or NULL (the default) 
#'  if the corresponding parameter is stationary).
#'  Parameters are: modified location, scale_0, shape, duration offset, duration exponent repectively.
25
26
#' @param mulink,siglink,shlink,thetalink,etalink Link functions for generalized linear 
#' modelling of the parameters, created with \code{\link{make.link}}.
27
28
#' @param init.vals list of length 5, giving initial values for all or some parameters
#' (order: mu, sigma, xi, theta, eta). If as.list(rep(NA,5)) (the default) is given, initial parameters are obtained 
29
#' internally by fitting the GEV separately for each duration and applying a linear model to obtain the 
30
#' duration dependency of the location and shape parameter.
31
#' Initial values for covariate parameters are assumed as 0 if not given.
32
33
#' @param theta_zero Logical value, indicating if theta should be estimated (FALSE, the default) or
#' should stay zero. 
34
35
36
37
38
39
#' @param show Logical; if TRUE (the default), print details of the fit.
#' @param method The optimization method used in \code{\link{optim}}.
#' @param maxit The maximum number of iterations.
#' @param ... Other control parameters for the optimization.
#' @return A list containing the following components. 
#' A subset of these components are printed after the fit. 
40
#' If \code{show} is TRUE, then assuming that successful convergence is indicated, 
41
#' the components nllh, mle and se are always printed. 
42
#' \item{nllh}{single numeric giving the negative log-likelihood value} 
43
44
#' \item{mle}{numeric vector giving the MLE's for the modified location, scale_0, shape, 
#' duration offset and duration exponent, resp.} 
45
#' \item{se}{numeric vector giving the standard errors for the MLE's (in the same order)}
46
#' \item{trans}{A logical indicator for a non-stationary fit.}
47
48
49
50
#' \item{model}{A list with components mul, sigl, shl, thetal and etal.}
#' \item{link}{A character vector giving inverse link functions.}
#' \item{conv}{The convergence code, taken from the list returned by \code{\link{optim}}. 
#' A zero indicates successful convergence.}
51
#' \item{data}{data is standardized to standard Gumbel.} 
52
#' \item{cov}{The covariance matrix.} 
53
#' \item{vals}{Parameter values for every data point.}
54
#' \item{init.vals}{Initial values that were used.}
55
#' \item{ds}{Durations for every data point.}
56
57
58
#' @seealso \code{\link{dgev.d}}, \code{\link{IDF.agg}}, \code{\link{gev.fit}}, \code{\link{optim}}
#' @export
#' @importFrom stats optim 
59
#' @importFrom stats make.link 
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
#' 
#' @examples 
#' # sampled random data from d-gev with covariates
#' # GEV parameters:
#' # mu = 4 + 0.2*cov1 +0.5*cov2
#' # sigma = 2+0.5*cov1
#' # xi = 0.5
#' # theta = 0
#' # eta = 0.5
#' 
#' data('example',package ='IDF')
#' 
#' gev.d.fit(xdat=example$dat,ds = example$d,ydat=as.matrix(example[,c('cov1','cov2')])
#' ,mul=c(1,2),sigl=1)

gev.d.fit<-
  function(xdat, ds, ydat = NULL, mul = NULL, sigl = NULL, shl = NULL, thetal = NULL, etal = NULL, 
77
78
           mulink = make.link("identity"), siglink = make.link("identity"), shlink = make.link("identity"),
           thetalink = make.link("identity"), etalink = make.link("identity"),  
79
           init.vals = as.list(rep(NA,5)), theta_zero = FALSE,
80
           show = TRUE, method = "Nelder-Mead", maxit = 10000, ...)
81
  {
82
83
84
    if (length(xdat) != length(ds)) {
      stop(paste0('The length of xdat is ',length(xdat),', but the length of ds is ',length(ds),'.'))
    } 
85
86
87
88
89
    z <- list()
    # number of parameters (betas) to estimate for each parameter: 
    npmu <- length(mul) + 1
    npsc <- length(sigl) + 1
    npsh <- length(shl) + 1
90
    npth <- ifelse(!theta_zero,length(thetal) + 1,0)
91
92
    npet <- length(etal) + 1
    z$trans <- FALSE  # indicates if fit is non-stationary
93
94
95
96
97
    z$model <- list(mul, sigl, shl, thetal, etal)
    z$link <- list(mulink=mulink, siglink=siglink, shlink=shlink, thetalink=thetalink, etalink=etalink)
    
    # test for NA values:
    if(any(is.na(xdat))) stop('xdat contains NA values. NA values need to be removed first.')
98
99
100
    # test for finite values:
    if(any(is.infinite(xdat))) stop('xdat contains non finite values. Inf and -Inf need to be removed first.')
    
101
102
    # test if covariates matrix is given correctly
    npar <- max(sapply(z$model,function(x){return(ifelse(is.null(x),0,max(x)))}))
103
    if(any(npar>ncol(ydat),npar>0 & is.null(ydat)))stop("Not enough columns in covariates matrix 'ydat'.")
104
    
105
    # initial values
106
107
108
    if(length(init.vals)!=5 | !is.list(init.vals)) {
      warning('Parameter init.vals is not used, because it is no list of length 5.')
      init.vals <- as.list(rep(NA,5))
109
110
    }
    if(!any(is.na(init.vals))){ #all initial values are given
111
      names(init.vals) <- c('mu','sigma','xi','theta','eta')
112
113
114
115
    }else if(any(!is.na(init.vals))) { #some initial values are given
      init.vals.user <- init.vals
      init.vals <- gev.d.init(xdat,ds,z$link) #calculate init.vals using gev.d.init
      for (i in 1:length(init.vals)){ #overwrite the calculated initial values with the ones given by the user
116
117
        if(!is.na(init.vals.user[[i]])) {
          init.vals[[i]]<-init.vals.user[[i]]
118
119
        }
      }
120
    }else{ #no initial values are given 
121
      init.vals <- gev.d.init(xdat,ds,z$link)
122
    }
123
124
    
    # generate covariates matrices: 
125
    if (is.null(mul)) { #stationary
126
      mumat <- as.matrix(rep(1, length(xdat)))
127
      muinit <- init.vals$mu
128
    }else { #non stationary
129
130
      z$trans <- TRUE
      mumat <- cbind(rep(1, length(xdat)), ydat[, mul])
131
      muinit <- c(init.vals$mu, rep(0, length(mul)))[1:npmu] #fill with 0s to length npmu
132
133
134
    }
    if (is.null(sigl)) {
      sigmat <- as.matrix(rep(1, length(xdat)))
135
      siginit <- init.vals$sigma
136
137
138
    }else {
      z$trans <- TRUE
      sigmat <- cbind(rep(1, length(xdat)), ydat[, sigl])
139
      siginit <- c(init.vals$sigma, rep(0, length(sigl)))[1:npsc]
140
141
142
    }
    if (is.null(shl)) {
      shmat <- as.matrix(rep(1, length(xdat)))
143
      shinit <- init.vals$xi 
144
145
146
    }else {
      z$trans <- TRUE
      shmat <- cbind(rep(1, length(xdat)), ydat[, shl])
147
      shinit <- c(init.vals$xi, rep(0, length(shl)))[1:npsh]
148
149
150
    }
    if (is.null(thetal)) {
      thmat <- as.matrix(rep(1, length(xdat)))
151
      thetainit <- init.vals$theta
152
153
154
    }else {
      z$trans <- TRUE
      thmat <- cbind(rep(1, length(xdat)), ydat[, thetal])
155
      thetainit <- c(init.vals$theta, rep(0, length(thetal)))[1:npth]
156
157
158
    }
    if (is.null(etal)) {
      etmat <- as.matrix(rep(1, length(xdat)))
159
      etainit <- init.vals$eta
160
161
162
    }else {
      z$trans <- TRUE
      etmat <- cbind(rep(1, length(xdat)), ydat[, etal])
163
      etainit <- c(init.vals$eta, rep(0, length(etal)))[1:npet]
164
    }
165
    
166

167
    if(!theta_zero){#When theta parameter is not included (default)
168
      init <- c(muinit, siginit, shinit, thetainit, etainit)
169
    }else{ #Do not return initial value for theta, if user does not want theta, as Hessian will fail.
170
171
      init <- c(muinit, siginit, shinit, etainit)
    }
172
173


174
175
176
    # function to calculate neg log-likelihood:
    gev.lik <- function(a) {
      # computes neg log lik of d-gev model
177
178
179
      mu <- mulink$linkinv(mumat %*% (a[1:npmu]))
      sigma <- siglink$linkinv(sigmat %*% (a[seq(npmu + 1, length = npsc)]))
      xi <- shlink$linkinv(shmat %*% (a[seq(npmu + npsc + 1, length = npsh)]))
180
      # Next line will set the theta likelihood as non-existent in case user requested it. 
181
      if(!theta_zero) {theta <- thetalink$linkinv(thmat %*% (a[seq(npmu + npsc + npsh + 1, length = npth)]))}
182
      eta <- etalink$linkinv(etmat %*% (a[seq(npmu + npsc + npsh + npth + 1, length = npet)]))
183
      
184
      ifelse(!theta_zero, ds.t <- ds+theta, ds.t <- ds) #Don't use theta if user requested not to have it.
185
186
187
188
      sigma.d <- sigma/(ds.t^eta)
      y <- xdat/sigma.d - mu
      y <- 1 + xi * y
      
189
      if(!theta_zero){ #When user wants to estimate theta parameter (default)
190
        if(any(eta <= 0) || any(theta < 0) || any(sigma.d <= 0) || any(y <= 0)) return(10^6)
191
      }else{ 
192
193
194
        if(any(eta <= 0) || any(sigma.d <= 0) || any(y <= 0)) return(10^6)
      }
      
195
196
197
198
199
200
201
202
203
204
      sum(log(sigma.d)) + sum(y^(-1/xi)) + sum(log(y) * (1/xi + 1))
    }
    
    
    # finding minimum of log-likelihood:
    x <- optim(init, gev.lik, hessian = TRUE, method = method,
               control = list(maxit = maxit, ...))
    
    # saving output parameters:
    z$conv <- x$convergence
205
206
207
    mut <- mulink$linkinv(mumat %*% (x$par[1:npmu]))
    sc0 <- siglink$linkinv(sigmat %*% (x$par[seq(npmu + 1, length = npsc)]))
    xi <- shlink$linkinv(shmat %*% (x$par[seq(npmu + npsc + 1, length = npsh)]))
208
209
210
211
212
    if(!theta_zero){ #When user does NOT set theta parameter to zero (default)
      theta <- thetalink$linkinv(thmat %*% (x$par[seq(npmu + npsc + npsh + 1, length = npth)]))
    }else{ #When user requests theta_parameter to be zero
      theta <- thetalink$linkinv(thmat %*% (0))
    }
213
    eta <- etalink$linkinv(etmat %*% (x$par[seq(npmu + npsc + npsh + npth + 1, length = npet)]))
214
    z$nllh <- x$value
215
    # normalize data to standard Gumbel:
216
217
218
    sc.d <- sc0/((ds+theta)^eta)
    z$data <-  - log(as.vector((1 + xi * (xdat/sc.d-mut))^(-1/xi))) 
    z$mle <- x$par
219
    test <- try(              # catch error 
220
    z$cov <- solve(x$hessian) # invert hessian to get estimation on var-covar-matrix
221
222
223
224
225
    ,silent = TRUE )
    if("try-error" %in% class(test)){
      warning("Hessian could not be inverted. NAs were produced.")
      z$cov <- matrix(NA,length(z$mle),length(z$mle))
        }
226
    z$se <- sqrt(diag(z$cov)) # sqrt(digonal entries) = standart error of mle's 
227
228
229
230
231
    if (!theta_zero) {#When theta parameter is returned (default)
      z$vals <- cbind(mut, sc0, xi, theta, eta)
    } else {#When theta parameter is not returned, asked by user
      z$vals <- cbind(mut, sc0, xi, eta)
    }
232
    z$init.vals <- unlist(init.vals)
233
234
235
236
237
    if(!theta_zero){ #When theta parameter is returned (default)
      colnames(z$vals) <-c('mut','sigma0','xi','theta','eta')
    } else { #When theta parameter is not returned, asked by user
      colnames(z$vals) <-c('mut','sigma0','xi','eta')
    }
238
    z$ds <- ds
239
    z$theta_zero <- theta_zero #Indicates if theta parameter was set to zero by user. 
240
    if(show) {
Laura Mack's avatar
Laura Mack committed
241
      if(z$trans) { # for nonstationary fit
242
        print(z[c(2, 4)]) # print model, link (3) , conv 
Laura Mack's avatar
Laura Mack committed
243
        # print names of link functions:
244
245
246
        cat('$link\n')
        print(c(z$link$mulink$name,z$link$siglink$name,z$link$shlink$name,z$link$thetalink$name,z$link$etalink$name))
        cat('\n')
Laura Mack's avatar
Laura Mack committed
247
248
      }else{print(z[4])} # for stationary fit print only conv
      if(!z$conv){ # if fit converged 
249
        print(z[c(5, 7, 9)]) # print nll, mle, se
Laura Mack's avatar
Laura Mack committed
250
      }
251
    }
252
    class(z) <- "gev.d.fit"
253
254
255
256
257
258
259
260
    invisible(z)
}


#### gev.d.init ####

# function to get initial values for gev.d.fit:
# obtain initial values 
261
# by fitting every duration separately
262
263

# possible ways to improve:
264
# take given initial values into account, if there are any
265
266
267
268
269
270
271
272
# xi -> mean vs. median ... how do we improve that?
# mu_tilde -> is not very good for small sample sizes yet
# improved inital value for eta, by fitting both mu~d and sigma~d in log-log scale

#' @title get initial values for gev.d.fit
#' @description obtain initial values by fitting every duration seperately
#' @param xdat vector of maxima for differnt durations
#' @param ds vector of durations belonging to maxima in xdat
273
#' @param link list of 5, link functions for parameters, created with \code{\link{make.link}}
274
275
#' @return list of initail values for mu_tilde, sigma_0, xi, eta
#' @importFrom stats lm 
276
#' @importFrom stats median 
277
278
279
#' @importFrom ismev gev.fit
#' @keywords internal 

280
gev.d.init <- function(xdat,ds,link){
281
282
283
  durs <- unique(ds)
  mles <- matrix(NA, nrow=length(durs), ncol= 3)
  for(i in 1:length(durs)){
284
285
    test <- try(fit <- gev.fit(xdat[ds==durs[i]],show = FALSE),silent = TRUE)
    if("try-error" %in% class(test) | fit$conv!=0){mles[i,] <- rep(NA,3)}else{mles[i,] <- fit$mle}
286
  }
287
  if(all(is.na(mles))){stop('Initial values could not be computed for this dataset.')}
288
  # get values for sig0 and eta (also mu_0) from linear model in log-log scale
289
290
  lmsig <- lm(log(mles[,2])~log(durs))
  lmmu <- lm(log(mles[,1])~log(durs))
291
292
  
  # sig0 <- exp Intercept
293
  siginit <- link$siglink$linkfun(exp(lmsig$coefficients[[1]]))
294
  # eta <- mean of negativ slopes 
295
  etainit <- link$etalink$linkfun(mean(c(-lmsig$coefficients[[2]],-lmmu$coefficients[[2]])))
296
297
  # mean of mu_d/sig_d 
  # could try:
298
  # mu0/sig0 = exp(lmmu$coefficients[[1]])/exp(lmsig$coefficients[[1]])
299
  muinit <- link$mulink$linkfun(median(c(mles[,1]/mles[,2]),na.rm = TRUE))
300
  # mean of shape parameters 
301
  shinit <- link$shlink$linkfun(median(mles[,3],na.rm = TRUE))
302
  thetainit <- link$thetalink$linkfun(0)
303
  
Jana Ulrich's avatar
Jana Ulrich committed
304
  return(list(mu=muinit,sigma=siginit,xi=shinit,theta=thetainit,eta=etainit))
305
306
}

307
308
309
310
#### gev.d.nll ####
#' computes negative log-likelihood of d-gev model
#'
#' @param xdat numeric vector containing observations
311
#' @param ds numeric vector containing coresponding durations (1/60 corresponds to 1 minute, 1 corresponds to 1 hour)
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
#' @param mut,sig0,xi,theta,eta numeric vectors containing corresponding mles for each of the parameters
#'
#' @return single value containing negative log likelihood 
#' @export
#'
#' @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])
gev.d.nll <- function(xdat,ds,mut,sig0,xi,theta,eta) {
  # computes neg log lik of d-gev model
328
329
  if(any(! c(length(ds),length(mut),length(sig0),length(xi),length(theta),length(eta)) %in% 
         c(1,length(xdat)))){
330
331
332
333
334
335
336
337
338
339
    warning('Input vectors differ in length, but must have the same length.')
  }
  
  ds.t <- ds+theta
  sigma.d <- sig0/(ds.t^eta)
  y <- xdat/sigma.d - mut
  y <- 1 + xi * y
  
  sum(log(sigma.d)) + sum(y^(-1/xi)) + sum(log(y) * (1/xi + 1))
}
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358

#### gev.d.diag ####

#' Diagnostic Plots for d-gev Models
#'
#' @description  Produces diagnostic plots for d-gev models using 
#' the output of the function \code{\link{gev.d.fit}}. Values for different durations can be plotted in 
#' different colors of with different symbols.
#' @param fit object returned by \code{\link{gev.d.fit}}
#' @param subset an optional vector specifying a subset of observations to be used in the plot
#' @param cols optional either one value or vector of same length as \code{unique(durations)} to
#' specify the colors of plotting points. 
#' The default uses the \code{rainbow} function.
#' @param pch optional either one value or vector of same length as \code{unique(durations)} containing
#' integers or symbols to specify the plotting points.
#' @param which string containing 'both', 'pp' or 'qq' to specify, which plots should be produced.
#' @param mfrow vector specifying layout of plots. If both plots should be produced seperately,
#' set to \code{c(1,1)}.
#' @param legend logical indicating if legends should be plotted
Jana Ulrich's avatar
Jana Ulrich committed
359
360
361
#' @param title character vector of length 2, giving the titles for the pp- and the qq-plot
#' @param emp.lab,mod.lab character string containing names for empirical and model axis
#' @param ... additional parameters passed on to the plotting function 
362
363
364
365
366
367
368
369
370
371
372
373
374
375
#'
#' @export
#' @importFrom graphics plot abline par title
#' @importFrom grDevices rainbow
#'
#' @examples
#' data('example',package ='IDF')
#' 
#' fit <- gev.d.fit(xdat=example$dat,ds = example$d,ydat=as.matrix(example[,c('cov1','cov2')])
#'                  ,mul=c(1,2),sigl=1)
#' # diagnostic plots for complete data                
#' gev.d.diag(fit)    
#' # diagnostic plots for subset of data (e.g. one station)            
#' gev.d.diag(fit,subset = example$cov1==1)
Jana Ulrich's avatar
Jana Ulrich committed
376
377
378
gev.d.diag <- function(fit,subset=NULL,cols=NULL,pch=NULL,which='both',mfrow=c(1,2),legend=TRUE,
                       title=c('Residual Probability Plot','Residual Quantile Plot'),
                       emp.lab='Empirical',mod.lab='Model',...){
379
  # check parameter:
Jana Ulrich's avatar
Jana Ulrich committed
380
  if(!is.element(which,c('both','pp','qq'))) stop("Parameter 'which'= ",which,
381
382
383
384
                                                 " but only 'both','pp' or 'qq' are allowed.")
  # subset data
  df <- data.frame(data=fit$data,ds=fit$ds)
  if(!is.null(subset))df <- df[subset,]
Jana Ulrich's avatar
Jana Ulrich committed
385
386
387
388
389
  # get single durations
  durs <- sort(unique(df$ds))
  # rescale durations to assign colors
  df$cval <- sapply(df$ds,function(d){which(durs==d)})

390
391
392
393
394
395
396
397
398
399
  # sort data 
  df <- df[order(df$data),]
  
  # plotting position
  n <- length(df$data)
  px <- (1:n)/(n + 1)

  # create plots:
  if(which=='both') par(mfrow=mfrow) # 2 subplots
  # colors and symbols
Jana Ulrich's avatar
Jana Ulrich committed
400
  if(is.null(cols))cols <- rainbow(length(durs))
401
402
403
404
405
  if(is.null(pch))pch <- df$cval
  
  if(which=='both'|which=='pp'){
    # pp
    plot(px, exp( - exp( - df$data)), xlab =
Jana Ulrich's avatar
Jana Ulrich committed
406
407
408
           emp.lab, ylab = mod.lab,col=cols[df$cval],pch=pch,...)
    abline(0, 1, col = 1,lwd=1)
    title(title[1])
Jana Ulrich's avatar
Jana Ulrich committed
409
    if(legend){legend('bottomright',legend = round(durs,digits = 2),pch=pch,
410
                      col = cols[1:length(durs)],title = 'Durations[h]',ncol = 2)}
411
412
413
414
  }
  if(which=='both'|which=='qq'){
    # qq
    plot( - log( - log(px)), df$data, ylab =
Jana Ulrich's avatar
Jana Ulrich committed
415
416
417
            emp.lab, xlab = mod.lab,col=cols[df$cval],pch=pch,...)
    abline(0, 1, col = 1,lwd=1)
    title(title[2])
Jana Ulrich's avatar
Jana Ulrich committed
418
    if(legend){legend('bottomright',legend = round(durs,digits = 2),pch=pch,
419
                      col = cols[1:length(durs)],title = 'Durations [h]',ncol = 2)}
420
421
422
423
424
425
426
427
428
429
430
  }
  if(which=='both') par(mfrow=c(1,1)) # reset par
}

#### gev.d.params ####

#' Calculate gev(d) parameters from \code{gev.d.fit} output
#'
#' @description function to calculate mut, sigma0, xi, theta, eta 
#' (modified location, scale, shape, duration offset, duration exponent) 
#' from results of \code{\link{gev.d.fit}} with covariates
431
#' @param fit fit object returned by \code{gev.d.fit} or \code{gev.fit}
432
#' @param ydat A matrix containing the covariates in the same order as used in \code{gev.d.fit}.
433
#' @seealso \code{\link{dgev.d}}
434
#' @return data.frame containing mu_tilde, sigma0, xi, theta, eta (or mu, sigma, xi for gev.fit objects)
435
436
437
438
439
440
#' @export
#' 
#' @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)
441
442
443
444
#' gev.d.params(fit = fit,ydat = cbind(c(0.9,1),c(0.5,1)))


gev.d.params <- function(fit,ydat){
445
446
  if(!class(fit)%in%c("gev.d.fit","gev.fit"))stop("'fit' must be an object returned by 'gev.d.fit' or 'gev.fit'.")
  if(!is.matrix(ydat))stop("'ydat' must be of class matrix.")
447
448
449
  n.par <- max(sapply(fit$model,function(x){return(ifelse(is.null(x),0,max(x)))}))
  if(n.par>ncol(ydat))stop("Covariates-Matrix 'ydat' has ",ncol(ydat), " columns, but ", n.par," are required.")
  
450
451
452
453
  # number of parameters
  npmu <- length(fit$model[[1]]) + 1
  npsc <- length(fit$model[[2]]) + 1
  npsh <- length(fit$model[[3]]) + 1
454
455
456
457
458
459
  if(class(fit)=="gev.d.fit"){
    if(!fit$theta_zero){npth <- length(fit$model[[4]]) + 1 #Including theta parameter (default)]
    }else{npth <- 0}#With no theta parameter, asked by user
    npet <- length(fit$model[[5]]) + 1
  }
 
460
  # inverse link functions
461
462
463
464
  if(class(fit)=="gev.d.fit"){
    mulink <- fit$link$mulink$linkinv
    siglink <- fit$link$siglink$linkinv
    shlink <- fit$link$shlink$linkinv
465
    if(!fit$theta_zero) thetalink <- fit$link$thetalink$linkinv
466
467
468
469
470
471
    etalink <- fit$link$etalink$linkinv
  }else{
    mulink <- eval(parse(text=fit$link))[[1]]
    siglink <- eval(parse(text=fit$link))[[2]]
    shlink <- eval(parse(text=fit$link))[[3]]
  }
472
  
473
474
475
476
  # covariates matrices
  mumat <- cbind(rep(1, dim(ydat)[1]), matrix(ydat[, fit$model[[1]]],dim(ydat)[1],npmu-1))
  sigmat <- cbind(rep(1, dim(ydat)[1]), matrix(ydat[, fit$model[[2]]],dim(ydat)[1],npsc-1))
  shmat <- cbind(rep(1, dim(ydat)[1]), matrix(ydat[, fit$model[[3]]],dim(ydat)[1],npsh-1))
477
478
479
480
  if(class(fit)=="gev.d.fit"){
    if(!fit$theta_zero){thmat <- cbind(rep(1, dim(ydat)[1]), matrix(ydat[, fit$model[[4]]],dim(ydat)[1],npth-1))}
    etmat <- cbind(rep(1, dim(ydat)[1]), matrix(ydat[, fit$model[[5]]],dim(ydat)[1],npet-1))
  }
481
  
482
483
484
485
  # calculate parameters
  mut <- mulink(mumat %*% (fit$mle[1:npmu]))
  sc0 <- siglink(sigmat %*% (fit$mle[seq(npmu + 1, length = npsc)]))
  xi <- shlink(shmat %*% (fit$mle[seq(npmu + npsc + 1, length = npsh)]))
486
487
488
  if(class(fit)=="gev.d.fit" ){
    if(!fit$theta_zero){theta <- thetalink(thmat %*% (fit$mle[seq(npmu + npsc + npsh + 1, length = npth)]))
    }else{theta <- rep(0,dim(ydat)[1])}}
489
  if(class(fit)=="gev.d.fit"){eta <- etalink(etmat %*% (fit$mle[seq(npmu + npsc + npsh + npth + 1, length = npet)]))}
490
  
491
  if(class(fit)=="gev.d.fit"){
492
    return(data.frame(mut=mut,sig0=sc0,xi=xi,theta=theta,eta=eta))
493
  }else{return(data.frame(mu=mut,sig=sc0,xi=xi))}
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
}


#### example data ####

#' Sampled data for duration dependent GEV
#'
#' A dataset containing:
#' \itemize{
#'   \item \code{$xdat}: 'annual' maxima values
#'   \item \code{$ds}: corresponding durations
#'   \item \code{$cov1}, \code{$cov2}: covariates}
#' GEV parameters:
#' \itemize{
#'   \item mu = 4 + 0.2*cov1 +0.5*cov2
#'   \item sigma = 2+0.5*cov1
#'   \item xi = 0.5
#'   \item theta = 0
#'   \item eta = 0.5}
#'
#' @docType data
#' @keywords datasets
#' @name example
#' @usage data('example',package ='IDF')
NULL