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

added function to transform data to standart gumbel -> gev.d2stgumbel

parent 83791606
...@@ -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.d2stdgumbel)
export(pgev.d) export(pgev.d)
export(qgev.d) export(qgev.d)
export(rgev.d) export(rgev.d)
......
...@@ -122,7 +122,7 @@ IDF.agg <- function(data,ds,na.accept = 0, ...@@ -122,7 +122,7 @@ IDF.agg <- function(data,ds,na.accept = 0,
#' 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)
#' par <- gev.d.params(fit = fit, ydat = cbind(1,1)) #' par <- gev.d.params(fit = fit, ydat = matrix(1,1,2))
#' IDF.plot(data = example[example$cov1==1,c("dat","d")],fitparams = unlist(par), #' IDF.plot(data = example[example$cov1==1,c("dat","d")],fitparams = unlist(par),
#' calc.dur = 2^(0:5),ylim=c(1,75),st.name = 'Example') #' calc.dur = 2^(0:5),ylim=c(1,75),st.name = 'Example')
IDF.plot <- function(data,fitparams,probs=c(0.5,0.9,0.99),calc.dur=NULL, IDF.plot <- function(data,fitparams,probs=c(0.5,0.9,0.99),calc.dur=NULL,
......
...@@ -167,7 +167,7 @@ gev.d.fit<- ...@@ -167,7 +167,7 @@ gev.d.fit<-
##################################################################################### #####################################################################################
# derivations of nll after d-gev-parameters (for boosting): # derivations of nll after d-gev-parameters (for boosting):
# get parameters from covariates and a (vector containing predictors) # get parameters from covariates (and a vector containing predictors)
# mu <- mulink(mumat %*% (a[1:npmu])) # mu <- mulink(mumat %*% (a[1:npmu]))
# sigma <- siglink(sigmat %*% (a[seq(npmu + 1, length = npsc)])) # sigma <- siglink(sigmat %*% (a[seq(npmu + 1, length = npsc)]))
# xi <- shlink(shmat %*% (a[seq(npmu + npsc + 1, length = npsh)])) # xi <- shlink(shmat %*% (a[seq(npmu + npsc + 1, length = npsh)]))
...@@ -363,10 +363,10 @@ gev.d.diag <- function(fit,subset=NULL,cols=NULL,pch=NULL,which='both',mfrow=c(1 ...@@ -363,10 +363,10 @@ gev.d.diag <- function(fit,subset=NULL,cols=NULL,pch=NULL,which='both',mfrow=c(1
#' @description function to calculate mut, sigma0, xi, theta, eta #' @description function to calculate mut, sigma0, xi, theta, eta
#' (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} or \code{gev.fit}
#' @param ydat A matrix containing the covariates in the same order as used in \code{gev.d.fit}. #' @param ydat A matrix containing the covariates in the same order as used in \code{gev.d.fit}.
#' @seealso \code{\link{dgev.d}} #' @seealso \code{\link{dgev.d}}
#' @return data.frame containing mu_tilde, sigma0, xi, theta, eta #' @return data.frame containing mu_tilde, sigma0, xi, theta, eta (or mu, sigma, xi for gev.fit objects)
#' @export #' @export
#' #'
#' @examples #' @examples
...@@ -377,56 +377,45 @@ gev.d.diag <- function(fit,subset=NULL,cols=NULL,pch=NULL,which='both',mfrow=c(1 ...@@ -377,56 +377,45 @@ 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)=="gev.d.fit")stop("'fit' must be an object returned by 'gev.d.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(is.null(ncol(ydat)))stop("'ydat' must be have class matrix.") 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.")
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.")
mut <- fit$mle[1] #ydat <- rbind(0,ydat) # no error in case ncols=1
if(is.null(fit$model[[1]])){i <- 1}else{
for(i in 1:length(fit$model[[1]])){ # number of parameters
cov <- fit$model[[1]][i] npmu <- length(fit$model[[1]]) + 1
mut <- mut + fit$mle[1+i]*ydat[,cov] npsc <- length(fit$model[[2]]) + 1
} npsh <- length(fit$model[[3]]) + 1
i <- i+1 if(class(fit)=="gev.d.fit"){npth <- length(fit$model[[4]]) + 1}
} if(class(fit)=="gev.d.fit"){npet <- length(fit$model[[5]]) + 1}
sig0 <- fit$mle[i+1] # link functions
if(is.null(fit$model[[2]])){j <- 1}else{ mulink <- eval(parse(text=fit$link))[[1]]
for(j in 1:length(fit$model[[2]])){ siglink <- eval(parse(text=fit$link))[[2]]
cov <- fit$model[[2]][j] shlink <- eval(parse(text=fit$link))[[3]]
sig0 <- sig0 + fit$mle[1+i+j]*ydat[,cov] if(class(fit)=="gev.d.fit"){thetalink <- eval(parse(text=fit$link))[[4]]}
} if(class(fit)=="gev.d.fit"){etalink <- eval(parse(text=fit$link))[[5]]}
j <- j+1
} # covariates matrices
mumat <- cbind(rep(1, dim(ydat)[1]), matrix(ydat[, fit$model[[1]]],dim(ydat)[1],npmu-1))
xi <- fit$mle[i+j+1] sigmat <- cbind(rep(1, dim(ydat)[1]), matrix(ydat[, fit$model[[2]]],dim(ydat)[1],npsc-1))
if(is.null(fit$model[[3]])){k <- 1}else{ shmat <- cbind(rep(1, dim(ydat)[1]), matrix(ydat[, fit$model[[3]]],dim(ydat)[1],npsh-1))
for(k in 1:length(fit$model[[3]])){ if(class(fit)=="gev.d.fit"){thmat <- cbind(rep(1, dim(ydat)[1]), matrix(ydat[, fit$model[[4]]],dim(ydat)[1],npth-1))}
cov <- fit$model[[3]][k] if(class(fit)=="gev.d.fit"){etmat <- cbind(rep(1, dim(ydat)[1]), matrix(ydat[, fit$model[[5]]],dim(ydat)[1],npet-1))}
xi <- xi + fit$mle[1+i+j+k]*ydat[,cov]
} # calculate parameters
k <- k+1 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)]))
theta <- fit$mle[i+j+k+1] if(class(fit)=="gev.d.fit"){theta <- thetalink(thmat %*% (fit$mle[seq(npmu + npsc + npsh + 1, length = npth)]))}
if(is.null(fit$model[[4]])){l <- 1}else{ if(class(fit)=="gev.d.fit"){eta <- etalink(etmat %*% (fit$mle[seq(npmu + npsc + npsh + npth + 1, length = npet)]))}
for(l in 1:length(fit$model[[4]])){
cov <- fit$model[[4]][l] if(class(fit)=="gev.d.fit"){return(data.frame(mut=mut,sig0=sc0,xi=xi,theta=theta,eta=eta))
theta <- theta + fit$mle[1+i+j+k+l]*ydat[,cov] }else{return(data.frame(mu=mut,sig=sc0,xi=xi))}
}
l <- l+1
}
eta <- fit$mle[i+j+k+l+1]
if(!is.null(fit$model[[5]])){
for(m in 1:length(fit$model[[5]])){
cov <- fit$model[[5]][m]
eta <- eta + fit$mle[1+i+j+k+l+m]*ydat[,cov]
}
}
return(data.frame(mut=mut,sig0=sig0,xi=xi,theta=theta,eta=eta))
} }
...@@ -470,6 +459,42 @@ gev.d.rl <- function(params,p.d){ ...@@ -470,6 +459,42 @@ gev.d.rl <- function(params,p.d){
} }
#### gev.d2stdgumbel ####
#' Transform data to standart gumbel
#'
#' @param xdat A vector containing maxima for different durations.
#' @param ds A vector of aggregation levels corresponding to the maxima in xdat.
#' @param params list of parameters mu_tilde, sigma0, xi, theta, eta
#' as obtained from \code{\link{gev.d.params}}
#'
#' @return Vector containing transformed data.
#' @export
#'
#' @examples
#' data('example',package = 'IDF')
#' # fit subset
#' ind <- sample(1:10, length(example$dat), replace=TRUE)
#' fit.subs <- gev.d.fit(example$dat[ind!=1],example$d[ind!=1]
#' ,ydat = as.matrix(example[ind!=1,c("cov1","cov2")])
#' ,mul = c(1,2),sigl = 1)
#' # calculate parameters for unfitted values
#' par <- gev.d.params(fit = fit.subs
#' ,ydat = as.matrix(example[ind==1,c("cov1","cov2")]))
#' # transform unfitted values to standart gumbel
#' sg.data <- gev.d2stdgumbel(xdat = example$dat[ind==1]
#' ,ds = example$d[ind==1],params = par)
#' # check unfitted values agains standart gumbel
#' gev.d.diag(data.frame(data=sg.data,ds=example$d[ind==1]),pch=20)
gev.d2stdgumbel <- function(xdat,ds,params){
sc.d <- params$sig0/((ds+params$theta)^params$eta)
sg.data <- - log(as.vector((1 + params$xi * (xdat/sc.d-params$mut))^(-1/params$xi)))
return(sg.data)
}
#### example data #### #### example data ####
#' Sampled data for duration dependent GEV #' Sampled data for duration dependent GEV
......
...@@ -42,7 +42,7 @@ Plotting of IDF curves at a chosen station ...@@ -42,7 +42,7 @@ Plotting of IDF curves at a chosen station
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)
par <- gev.d.params(fit = fit, ydat = cbind(1,1)) par <- gev.d.params(fit = fit, ydat = matrix(1,1,2))
IDF.plot(data = example[example$cov1==1,c("dat","d")],fitparams = unlist(par), IDF.plot(data = example[example$cov1==1,c("dat","d")],fitparams = unlist(par),
calc.dur = 2^(0:5),ylim=c(1,75),st.name = 'Example') calc.dur = 2^(0:5),ylim=c(1,75),st.name = 'Example')
} }
...@@ -7,12 +7,12 @@ ...@@ -7,12 +7,12 @@
gev.d.params(fit, ydat) 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} or \code{gev.fit}}
\item{ydat}{A matrix containing the covariates in the same order as used in \code{gev.d.fit}.} \item{ydat}{A matrix containing the covariates in the same order as used in \code{gev.d.fit}.}
} }
\value{ \value{
data.frame containing mu_tilde, sigma0, xi, theta, eta data.frame containing mu_tilde, sigma0, xi, theta, eta (or mu, sigma, xi for gev.fit objects)
} }
\description{ \description{
function to calculate mut, sigma0, xi, theta, eta function to calculate mut, sigma0, xi, theta, eta
......
% Generated by roxygen2: do not edit by hand
% Please edit documentation in R/gevdfit.R
\name{gev.d2stdgumbel}
\alias{gev.d2stdgumbel}
\title{Transform data to standart gumbel}
\usage{
gev.d2stdgumbel(xdat, ds, params)
}
\arguments{
\item{xdat}{A vector containing maxima for different durations.}
\item{ds}{A vector of aggregation levels corresponding to the maxima in xdat.}
\item{params}{list of parameters mu_tilde, sigma0, xi, theta, eta
as obtained from \code{\link{gev.d.params}}}
}
\value{
Vector containing transformed data.
}
\description{
Transform data to standart gumbel
}
\examples{
data('example',package = 'IDF')
# fit subset
ind <- sample(1:10, length(example$dat), replace=TRUE)
fit.subs <- gev.d.fit(example$dat[ind!=1],example$d[ind!=1]
,ydat = as.matrix(example[ind!=1,c("cov1","cov2")])
,mul = c(1,2),sigl = 1)
# calculate parameters for unfitted values
par <- gev.d.params(fit = fit.subs
,ydat = as.matrix(example[ind==1,c("cov1","cov2")]))
# transform unfitted values to standart gumbel
sg.data <- gev.d2stdgumbel(xdat = example$dat[ind==1]
,ds = example$d[ind==1],params = par)
# check unfitted values agains standart gumbel
gev.d.diag(data.frame(data=sg.data,ds=example$d[ind==1]),pch=20)
}
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