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

small changes

parent ee93972a
...@@ -7,6 +7,7 @@ export(gev.d.diag) ...@@ -7,6 +7,7 @@ export(gev.d.diag)
export(gev.d.fit) export(gev.d.fit)
export(gev.d.params) export(gev.d.params)
export(gev.d.rl) export(gev.d.rl)
export(gev.d.rl2)
export(gev.d2stdgumbel) export(gev.d2stdgumbel)
export(pgev.d) export(pgev.d)
export(qgev.d) export(qgev.d)
......
...@@ -89,6 +89,9 @@ gev.d.fit<- ...@@ -89,6 +89,9 @@ gev.d.fit<-
# calculate initial values for mu.d, sigma_0, xi, eta using IDF.init: (thetainit=0) # calculate initial values for mu.d, sigma_0, xi, eta using IDF.init: (thetainit=0)
init.vals <- gev.d.init(xdat,ds,ifelse(is.null(thetainit),0,thetainit[1])) init.vals <- gev.d.init(xdat,ds,ifelse(is.null(thetainit),0,thetainit[1]))
# TODO: transform initial values with link function
# use make.link()
# generate covariates matrices: # generate covariates matrices:
if (is.null(mul)) { if (is.null(mul)) {
...@@ -144,6 +147,7 @@ gev.d.fit<- ...@@ -144,6 +147,7 @@ gev.d.fit<-
z$model <- list(mul, sigl, shl, thetal, etal) z$model <- list(mul, sigl, shl, thetal, etal)
z$link <- deparse(substitute(c(mulink, siglink, shlink, thetalink, etalink))) z$link <- deparse(substitute(c(mulink, siglink, shlink, thetalink, etalink)))
# TODO: store as functions not as strings!
init <- c(muinit, siginit, shinit, thetainit, etainit) init <- c(muinit, siginit, shinit, thetainit, etainit)
# function to calculate neg log-likelihood: # function to calculate neg log-likelihood:
...@@ -252,7 +256,8 @@ gev.d.init <- function(xdat,ds,thetainit){ ...@@ -252,7 +256,8 @@ gev.d.init <- function(xdat,ds,thetainit){
durs <- unique(ds) durs <- unique(ds)
mles <- matrix(NA, nrow=length(durs), ncol= 3) mles <- matrix(NA, nrow=length(durs), ncol= 3)
for(i in 1:length(durs)){ for(i in 1:length(durs)){
mles[i,] <- gev.fit(xdat[ds==durs[i]],show = FALSE)$mle test <- try(mles[i,] <- gev.fit(xdat[ds==durs[i]],show = FALSE)$mle,silent = TRUE)
if("try-error" %in% class(test)){mles[i,] <- rep(NA,3)}
} }
# get values for sig0 and eta (also mu_0) from linear model in log-log scale # get values for sig0 and eta (also mu_0) from linear model in log-log scale
lmsig <- lm(log(mles[,2])~log(durs+thetainit)) lmsig <- lm(log(mles[,2])~log(durs+thetainit))
...@@ -265,9 +270,9 @@ gev.d.init <- function(xdat,ds,thetainit){ ...@@ -265,9 +270,9 @@ gev.d.init <- function(xdat,ds,thetainit){
# mean of mu_d/sig_d # mean of mu_d/sig_d
# could try: # could try:
# mu0/sig0 is also an estimate but needs to be weighted in mean # mu0/sig0 is also an estimate but needs to be weighted in mean
muinit <- mean(c(mles[,1]/mles[,2])) #exp(lmmu$coefficients[[1]])/exp(lmsig$coefficients[[1]]) muinit <- mean(c(mles[,1]/mles[,2]),na.rm = TRUE) #exp(lmmu$coefficients[[1]])/exp(lmsig$coefficients[[1]])
# mean of shape parameters # mean of shape parameters
shinit <- mean(mles[,3]) shinit <- mean(mles[,3],na.rm = TRUE)
return(list(mu=muinit,sigma=siginit,xi=shinit,eta=etainit)) return(list(mu=muinit,sigma=siginit,xi=shinit,eta=etainit))
} }
...@@ -378,9 +383,10 @@ gev.d.diag <- function(fit,subset=NULL,cols=NULL,pch=NULL,which='both',mfrow=c(1 ...@@ -378,9 +383,10 @@ gev.d.diag <- function(fit,subset=NULL,cols=NULL,pch=NULL,which='both',mfrow=c(1
gev.d.params <- function(fit,ydat){ gev.d.params <- function(fit,ydat){
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(!class(fit)%in%c("gev.d.fit","gev.fit"))stop("'fit' must be an object returned by 'gev.d.fit' or 'gev.fit'.")
if(!fit$trans){warning('No glm for parameters. Max. likelihood estimates are returned.')
return(fit$mle)}
if(!is.matrix(ydat))stop("'ydat' must be of class matrix.") if(!is.matrix(ydat))stop("'ydat' must be of class matrix.")
if(!fit$trans){warning('No vglm for parameters. Max. likelihood estimates are returned.')
return(data.frame(mut=rep(fit$mle[1],dim(ydat)[1]),sig0=fit$mle[2],xi=fit$mle[3]
,theta=fit$mle[4],eta=fit$mle[5]))}
n.par <- max(sapply(fit$model,function(x){return(ifelse(is.null(x),0,max(x)))})) 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.") if(n.par>ncol(ydat))stop("Covariates-Matrix 'ydat' has ",ncol(ydat), " columns, but ", n.par," are required.")
...@@ -422,9 +428,9 @@ gev.d.params <- function(fit,ydat){ ...@@ -422,9 +428,9 @@ gev.d.params <- function(fit,ydat){
#### gev.d.rl #### #### gev.d.rl ####
#' Calculate (spatial) Returnlevel #' Calculate Returnlevel
#' #'
#' @description calculate (spatial) Returnlevel for chosen duration and return period #' @description calculate Returnlevel for chosen duration and return period
#' from \code{\link{gev.d.fit}} parameters #' from \code{\link{gev.d.fit}} parameters
#' @param params list of parameters mu_tilde, sigma0, xi, theta, eta (single values and/or compatible matrices) #' @param 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}} #' as obtained from \code{\link{gev.d.fit}} or \code{\link{gev.d.params}}
...@@ -458,6 +464,108 @@ gev.d.rl <- function(params,p.d){ ...@@ -458,6 +464,108 @@ gev.d.rl <- function(params,p.d){
return(rl) return(rl)
} }
### alternativ function that can also calculate delta-method-sd:
#' Calculate Returnlevel with standart deviation (through delta method)
#'
#'@description calculate returnlevel for chosen duration and return period
#' for \code{\link{gev.d.fit}} objects
#' @param p single numeric value containing the non-exceedance probability p=1-1/T
#' @param fit \code{\link{gev.d.fit}} object
#' @param ds numeric vector of durations for which to calculate the return levels
#' @param ydat matrix containing the covariates in the same order as used in \code{\link{gev.d.fit}}
#' @param sd logical, if TRUE, standart deviation of the returnlevels obtained with the delta method
#' are also returned
#'
#' @return 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)
#' @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)
#' # 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)
gev.d.rl2<-function(p,fit,ds,ydat,sd=FALSE){
if(length(p)>1){
warning("'p' can only be a single value. Only first element of 'p' is used.")
p <- p[1]
}
# parameter estimates
if(!fit$trans){ # model without covariates
mut <- fit$mle[1]
sig0 <- fit$mle[2]
xi <- fit$mle[3]
theta <- fit$mle[4]
eta <- fit$mle[5]
}else{ # model with covariates
n.cov <- max(sapply(fit$model,function(x){return(ifelse(is.null(x),0,max(x)))}))
if(is.null(ydat)){ # if no ydat, 0's will be used as values for covariates
warning("No covariates-matrix was given. Zeros are used as values of covariates.")
ydat <- matrix(0,1,n.cov)
}
if(!is.matrix(ydat)){
warning("'ydat' was converted to a matrix.")
ydat <- matrix(ydat,1,n.cov)
}
yrows <- nrow(ydat)
ds <- matrix(ds,nrow(ydat),length(ds),byrow = TRUE)
params <- gev.d.params(fit,ydat) # calculate parameters
ydat <- cbind(rep(1,yrows),ydat) # add first column of ones
mut <- params$mut
sig0 <- params$sig0
xi <- params$xi
theta <- params$theta
eta <- params$eta
}
# quantile
y.xi <- (-1/log(p))^xi
dtheta <- ds+theta
sigma <- sig0*dtheta^(-eta)
q <- mut*sigma+sigma/xi*(y.xi-1)
if(fit$trans){colnames(q) <- ds[1,]}
if(!sd){return(q)} # only qantiles
# partial derivatives
dq.dmut <- sigma
dq.dsig0 <- mut*dtheta^(-eta)-dtheta^(-eta)/xi*(y.xi-1)
dq.dxi <- -sigma/xi^2*(y.xi-1)-sigma/xi*y.xi*log(-log(p))
dq.dtheta <- mut*sig0*dtheta^(-eta-1)-sig0*(-eta*dtheta^(-eta-1))/xi*(y.xi-1)
dq.deta <- -mut*sigma*log(dtheta)+sigma/xi*log(dtheta)*(y.xi-1)
if(!fit$trans){
dq.dpar <- rbind(dq.dmut,dq.dsig0,dq.dxi,dq.dtheta,dq.deta)
}else{
dq.dpar <- lapply(1:yrows,function(i){
rbind(dq.dmut[i,],dq.dsig0[i,],dq.dxi[i,],dq.dtheta[i,],dq.deta[i,])})
}
# calculate sd
if(!fit$trans){ # model without covariates
sd<-apply(dq.dpar,2,function(d.vec) sqrt(t(d.vec) %*% fit$cov %*% d.vec))
}else{ # model with covariates
n.par <- sapply(fit$model,length)+1
cov.ind <- sapply(fit$model,function(vec){c(1,vec+1)})
cov.ind <- do.call(c,cov.ind) # index of ydat
calc.sd <- function(j){
inner <- ydat[j,cov.ind] # inner derivative
outer <- dq.dpar[[j]]
outer <- sapply(1:5,function(i){ # outer derivative
matrix(rep(outer[i,],n.par[i]),n.par[i],length(ds[1,]),byrow = TRUE)})
outer <- do.call(rbind,outer)
sd<-apply(outer,2,function(d.vec) sqrt(t(d.vec*inner) %*% fit$cov %*% (d.vec*inner)))
}
sd <- lapply(1:yrows,calc.sd) # calculate sd for every row in ydat
sd <- do.call(rbind,sd)
colnames(sd) <- ds[1,]
}
return(list(q=q,sd=sd)) # qantiles and standart deviation
}
#### gev.d2stdgumbel #### #### gev.d2stdgumbel ####
......
...@@ -2,7 +2,7 @@ ...@@ -2,7 +2,7 @@
% Please edit documentation in R/gevdfit.R % Please edit documentation in R/gevdfit.R
\name{gev.d.rl} \name{gev.d.rl}
\alias{gev.d.rl} \alias{gev.d.rl}
\title{Calculate (spatial) Returnlevel} \title{Calculate Returnlevel}
\usage{ \usage{
gev.d.rl(params, p.d) gev.d.rl(params, p.d)
} }
...@@ -18,7 +18,7 @@ one return level value or matrix with return levels (depending on input to param ...@@ -18,7 +18,7 @@ one return level value or matrix with return levels (depending on input to param
unit: e.g. mm/h unit: e.g. mm/h
} }
\description{ \description{
calculate (spatial) Returnlevel for chosen duration and return period calculate Returnlevel for chosen duration and return period
from \code{\link{gev.d.fit}} parameters from \code{\link{gev.d.fit}} parameters
} }
\examples{ \examples{
......
% 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)
}
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