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

- changed gev.d.params input to covariates matrix instread of list

- catch solve(x$hessian) error in gev.d.fit
parent 5212b6cd
...@@ -202,7 +202,13 @@ gev.d.fit<- ...@@ -202,7 +202,13 @@ gev.d.fit<-
sc.d <- sc0/((ds+theta)^eta) sc.d <- sc0/((ds+theta)^eta)
z$data <- - log(as.vector((1 + xi * (xdat/sc.d-mut))^(-1/xi))) z$data <- - log(as.vector((1 + xi * (xdat/sc.d-mut))^(-1/xi)))
z$mle <- x$par z$mle <- x$par
test <- try( # catch error
z$cov <- solve(x$hessian) # invert hessian to get estimation on var-covar-matrix z$cov <- solve(x$hessian) # invert hessian to get estimation on var-covar-matrix
,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))
}
z$se <- sqrt(diag(z$cov)) # sqrt(digonal entries) = standart error of mle's z$se <- sqrt(diag(z$cov)) # sqrt(digonal entries) = standart error of mle's
z$vals <- cbind(mut, sc0, xi, theta, eta) z$vals <- cbind(mut, sc0, xi, theta, eta)
colnames(z$vals) <- c('mut','sigma0','xi','theta','eta') colnames(z$vals) <- c('mut','sigma0','xi','theta','eta')
...@@ -358,24 +364,29 @@ gev.d.diag <- function(fit,subset=NULL,cols=NULL,pch=NULL,which='both',mfrow=c(1 ...@@ -358,24 +364,29 @@ gev.d.diag <- function(fit,subset=NULL,cols=NULL,pch=NULL,which='both',mfrow=c(1
#' (modified location, scale, shape, duration offset, duration exponent) #' (modified location, scale, shape, duration offset, duration exponent)
#' from results of \code{\link{gev.d.fit}} with covariates #' from results of \code{\link{gev.d.fit}} with covariates
#' @param fit fit object returned by \code{gev.d.fit} #' @param fit fit object returned by \code{gev.d.fit}
#' @param cov.list list of covariates. Either single values - to calculate #' @param ydat A matrix containing the covariates in the same order as used in \code{gev.d.fit}.
#' parameters at a single station or compatible vectors or matrices - to calculate
#' parameters on a grid
#' @seealso \code{\link{dgev.d}} #' @seealso \code{\link{dgev.d}}
#' @return list containing mu_tilde, sigma0, xi, theta, eta #' @return data.frame containing mu_tilde, sigma0, xi, theta, eta
#' @export #' @export
#' #'
#' @examples #' @examples
#' data('example',package = 'IDF') #' data('example',package = 'IDF')
#' fit <- gev.d.fit(example$dat,example$d,ydat = as.matrix(example[,c("cov1","cov2")]) #' fit <- gev.d.fit(example$dat,example$d,ydat = as.matrix(example[,c("cov1","cov2")])
#' ,mul = c(1,2),sigl = 1) #' ,mul = c(1,2),sigl = 1)
#' gev.d.params(fit = fit,cov.list = list(0.9,0.5)) #' gev.d.params(fit = fit,ydat = cbind(c(0.9,1),c(0.5,1)))
gev.d.params <- function(fit,cov.list){
gev.d.params <- function(fit,ydat){
if(!class(fit)=="gev.d.fit")stop("'fit' must be an object returned by 'gev.d.fit'.")
if(is.null(ncol(ydat)))stop("'ydat' must be have class matrix.")
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.")
mut <- fit$mle[1] mut <- fit$mle[1]
if(is.null(fit$model[[1]])){i <- 1}else{ if(is.null(fit$model[[1]])){i <- 1}else{
for(i in 1:length(fit$model[[1]])){ for(i in 1:length(fit$model[[1]])){
cov <- fit$model[[1]][i] cov <- fit$model[[1]][i]
mut <- mut + fit$mle[1+i]*cov.list[[cov]] mut <- mut + fit$mle[1+i]*ydat[,cov]
} }
i <- i+1 i <- i+1
} }
...@@ -384,7 +395,7 @@ gev.d.params <- function(fit,cov.list){ ...@@ -384,7 +395,7 @@ gev.d.params <- function(fit,cov.list){
if(is.null(fit$model[[2]])){j <- 1}else{ if(is.null(fit$model[[2]])){j <- 1}else{
for(j in 1:length(fit$model[[2]])){ for(j in 1:length(fit$model[[2]])){
cov <- fit$model[[2]][j] cov <- fit$model[[2]][j]
sig0 <- sig0 + fit$mle[1+i+j]*cov.list[[cov]] sig0 <- sig0 + fit$mle[1+i+j]*ydat[,cov]
} }
j <- j+1 j <- j+1
} }
...@@ -393,7 +404,7 @@ gev.d.params <- function(fit,cov.list){ ...@@ -393,7 +404,7 @@ gev.d.params <- function(fit,cov.list){
if(is.null(fit$model[[3]])){k <- 1}else{ if(is.null(fit$model[[3]])){k <- 1}else{
for(k in 1:length(fit$model[[3]])){ for(k in 1:length(fit$model[[3]])){
cov <- fit$model[[3]][k] cov <- fit$model[[3]][k]
xi <- xi + fit$mle[1+i+j+k]*cov.list[[cov]] xi <- xi + fit$mle[1+i+j+k]*ydat[,cov]
} }
k <- k+1 k <- k+1
} }
...@@ -402,7 +413,7 @@ gev.d.params <- function(fit,cov.list){ ...@@ -402,7 +413,7 @@ gev.d.params <- function(fit,cov.list){
if(is.null(fit$model[[4]])){l <- 1}else{ if(is.null(fit$model[[4]])){l <- 1}else{
for(l in 1:length(fit$model[[4]])){ for(l in 1:length(fit$model[[4]])){
cov <- fit$model[[4]][l] cov <- fit$model[[4]][l]
theta <- theta + fit$mle[1+i+j+k+l]*cov.list[[cov]] theta <- theta + fit$mle[1+i+j+k+l]*ydat[,cov]
} }
l <- l+1 l <- l+1
} }
...@@ -411,14 +422,15 @@ gev.d.params <- function(fit,cov.list){ ...@@ -411,14 +422,15 @@ gev.d.params <- function(fit,cov.list){
if(!is.null(fit$model[[5]])){ if(!is.null(fit$model[[5]])){
for(m in 1:length(fit$model[[5]])){ for(m in 1:length(fit$model[[5]])){
cov <- fit$model[[5]][m] cov <- fit$model[[5]][m]
eta <- eta + fit$mle[1+i+j+k+l+m]*cov.list[[cov]] eta <- eta + fit$mle[1+i+j+k+l+m]*ydat[,cov]
} }
} }
return(list(mut=mut,sig0=sig0,xi=xi,theta=theta,eta=eta)) return(data.frame(mut=mut,sig0=sig0,xi=xi,theta=theta,eta=eta))
} }
#### gev.d.rl #### #### gev.d.rl ####
#' Calculate (spatial) Returnlevel #' Calculate (spatial) Returnlevel
...@@ -431,7 +443,7 @@ gev.d.params <- function(fit,cov.list){ ...@@ -431,7 +443,7 @@ gev.d.params <- function(fit,cov.list){
#' and one value for the duration at which to calculate the return level #' and one value for the duration at which to calculate the return level
#' #'
#' @return one return level value or matrix with return levels (depending on input to params) #' @return one return level value or matrix with return levels (depending on input to params)
#' unit: e.g. mm/(given duration) #' unit: e.g. mm/h
#' @export #' @export
#' #'
#' @examples #' @examples
...@@ -439,21 +451,22 @@ gev.d.params <- function(fit,cov.list){ ...@@ -439,21 +451,22 @@ gev.d.params <- function(fit,cov.list){
#' fit <- gev.d.fit(example$dat,example$d,ydat = as.matrix(example[,c("cov1","cov2")]) #' fit <- gev.d.fit(example$dat,example$d,ydat = as.matrix(example[,c("cov1","cov2")])
#' ,mul = c(1,2),sigl = 1) #' ,mul = c(1,2),sigl = 1)
#' ### calculate rl on grid: #' ### calculate rl on grid:
#' # create matrixes for the covariates values #' # covariates values
#' cov1 <- matrix(seq(-1,1,0.1),ncol=11,nrow=21) #' cov1 <- rep(seq(-1,1,0.1),11)
#' cov2 <- matrix(seq(0,1,0.1),ncol=11,nrow=21,byrow = TRUE) #' cov2 <- rep(seq(0,1,0.1),each=21)
#' # calculate parameters of d-gev on grid, based on output of gev.d.fit #' # calculate parameters of d-gev on grid, based on output of gev.d.fit
#' par <- gev.d.params(fit = fit,cov.list = list(cov1,cov2)) #' par <- gev.d.params(fit = fit,ydat = cbind(cov1,cov2))
#' # calculate 100 year (p=0.99) rl for a duration of d=24 hours #' # calculate 100 year (p=0.99) rl for a duration of d=24 hours
#' rl <- gev.d.rl(params = par,p.d = c(0.99,24)) #' rl <- gev.d.rl(params = par,p.d = c(0.99,24))
#' # plot of 'spatial rl': #' dim(rl) <- c(21,11)
#' # rl map:
#' image(x=seq(-1,1,0.1),y=seq(0,1,0.1),z=rl,xlab = 'cov1', ylab = 'cov2') #' image(x=seq(-1,1,0.1),y=seq(0,1,0.1),z=rl,xlab = 'cov1', ylab = 'cov2')
gev.d.rl <- function(params,p.d){ gev.d.rl <- function(params,p.d){
sigma.d <- params[[2]]/((p.d[2]+params[[4]])^params[[5]]) sigma.d <- params[[2]]/((p.d[2]+params[[4]])^params[[5]])
mu <- params[[1]]*sigma.d mu <- params[[1]]*sigma.d
yt <- -1/log(p.d[1]) yt <- -1/log(p.d[1])
rl <- mu+sigma.d/params[[3]]*(yt^params[[3]]-1) rl <- mu+sigma.d/params[[3]]*(yt^params[[3]]-1)
return(rl*p.d[2]) return(rl)
} }
......
...@@ -4,17 +4,15 @@ ...@@ -4,17 +4,15 @@
\alias{gev.d.params} \alias{gev.d.params}
\title{Calculate gev(d) parameters from \code{gev.d.fit} output} \title{Calculate gev(d) parameters from \code{gev.d.fit} output}
\usage{ \usage{
gev.d.params(fit, cov.list) gev.d.params(fit, ydat)
} }
\arguments{ \arguments{
\item{fit}{fit object returned by \code{gev.d.fit}} \item{fit}{fit object returned by \code{gev.d.fit}}
\item{cov.list}{list of covariates. Either single values - to calculate \item{ydat}{A matrix containing the covariates in the same order as used in \code{gev.d.fit}.}
parameters at a single station or compatible vectors or matrices - to calculate
parameters on a grid}
} }
\value{ \value{
list containing mu_tilde, sigma0, xi, theta, eta data.frame containing mu_tilde, sigma0, xi, theta, eta
} }
\description{ \description{
function to calculate mut, sigma0, xi, theta, eta function to calculate mut, sigma0, xi, theta, eta
...@@ -25,7 +23,7 @@ from results of \code{\link{gev.d.fit}} with covariates ...@@ -25,7 +23,7 @@ from results of \code{\link{gev.d.fit}} with covariates
data('example',package = 'IDF') data('example',package = 'IDF')
fit <- gev.d.fit(example$dat,example$d,ydat = as.matrix(example[,c("cov1","cov2")]) fit <- gev.d.fit(example$dat,example$d,ydat = as.matrix(example[,c("cov1","cov2")])
,mul = c(1,2),sigl = 1) ,mul = c(1,2),sigl = 1)
gev.d.params(fit = fit,cov.list = list(0.9,0.5)) gev.d.params(fit = fit,ydat = cbind(c(0.9,1),c(0.5,1)))
} }
\seealso{ \seealso{
\code{\link{dgev.d}} \code{\link{dgev.d}}
......
...@@ -15,7 +15,7 @@ and one value for the duration at which to calculate the return level} ...@@ -15,7 +15,7 @@ and one value for the duration at which to calculate the return level}
} }
\value{ \value{
one return level value or matrix with return levels (depending on input to params) one return level value or matrix with return levels (depending on input to params)
unit: e.g. mm/(given duration) unit: e.g. mm/h
} }
\description{ \description{
calculate (spatial) Returnlevel for chosen duration and return period calculate (spatial) Returnlevel for chosen duration and return period
...@@ -26,13 +26,14 @@ data('example',package = 'IDF') ...@@ -26,13 +26,14 @@ data('example',package = 'IDF')
fit <- gev.d.fit(example$dat,example$d,ydat = as.matrix(example[,c("cov1","cov2")]) fit <- gev.d.fit(example$dat,example$d,ydat = as.matrix(example[,c("cov1","cov2")])
,mul = c(1,2),sigl = 1) ,mul = c(1,2),sigl = 1)
### calculate rl on grid: ### calculate rl on grid:
# create matrixes for the covariates values # covariates values
cov1 <- matrix(seq(-1,1,0.1),ncol=11,nrow=21) cov1 <- rep(seq(-1,1,0.1),11)
cov2 <- matrix(seq(0,1,0.1),ncol=11,nrow=21,byrow = TRUE) cov2 <- rep(seq(0,1,0.1),each=21)
# calculate parameters of d-gev on grid, based on output of gev.d.fit # calculate parameters of d-gev on grid, based on output of gev.d.fit
par <- gev.d.params(fit = fit,cov.list = list(cov1,cov2)) par <- gev.d.params(fit = fit,ydat = cbind(cov1,cov2))
# calculate 100 year (p=0.99) rl for a duration of d=24 hours # calculate 100 year (p=0.99) rl for a duration of d=24 hours
rl <- gev.d.rl(params = par,p.d = c(0.99,24)) rl <- gev.d.rl(params = par,p.d = c(0.99,24))
# plot of 'spatial rl': dim(rl) <- c(21,11)
# rl map:
image(x=seq(-1,1,0.1),y=seq(0,1,0.1),z=rl,xlab = 'cov1', ylab = 'cov2') image(x=seq(-1,1,0.1),y=seq(0,1,0.1),z=rl,xlab = 'cov1', ylab = 'cov2')
} }
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