Commit 12f162d0 authored by jana ulrich's avatar jana ulrich
Browse files

added confidence intervals to diagnostic plotting function

parent 9ed4faa1
Package: IDF Package: IDF
Type: Package Type: Package
Title: Estimation and Plotting of IDF Curves Title: Estimation and Plotting of IDF Curves
Version: 2.1.0 Version: 2.1.0.9000
Date: 2021-04-07 Date: 2021-05-12
Authors@R: c(person("Jana", "Ulrich", email = "jana.ulrich@met.fu-berlin.de", role = c("aut", "cre")), Authors@R: c(person("Jana", "Ulrich", email = "jana.ulrich@met.fu-berlin.de", role = c("aut", "cre")),
person("Laura","Mack", email= "laura.mack@fu-berlin.de",role="ctb"), person("Laura","Mack", email= "laura.mack@fu-berlin.de",role="ctb"),
person("Oscar E.","Jurado", email= "jurado@zedat.fu-berlin.de",role="ctb"), person("Oscar E.","Jurado", email= "jurado@zedat.fu-berlin.de",role="ctb"),
......
# IDF 2.1.0.9000
- added confidence intervals for diagnostic plots (gev.d.diag)
# IDF 2.1.0 # IDF 2.1.0
## New features: ## New features:
......
...@@ -479,11 +479,13 @@ gev.d.lik <- function(xdat,ds,mut,sigma0,xi,theta,eta,log=FALSE,tau=0,eta2=NULL) ...@@ -479,11 +479,13 @@ gev.d.lik <- function(xdat,ds,mut,sigma0,xi,theta,eta,log=FALSE,tau=0,eta2=NULL)
#' @param legend logical indicating if legends should be plotted #' @param legend logical indicating if legends should be plotted
#' @param title character vector of length 2, giving the titles for the pp- and the qq-plot #' @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 emp.lab,mod.lab character string containing names for empirical and model axis
#' @param ci logical indicating whether 0.95 confidence intervals should be plotted
#' @param ... additional parameters passed on to the plotting function #' @param ... additional parameters passed on to the plotting function
#' #'
#' @export #' @export
#' @importFrom graphics plot abline par title #' @importFrom graphics plot abline par title
#' @importFrom grDevices rainbow #' @importFrom grDevices rainbow
#' @importFrom evd rgev
#' #'
#' @examples #' @examples
#' data('example',package ='IDF') #' data('example',package ='IDF')
...@@ -496,10 +498,10 @@ gev.d.lik <- function(xdat,ds,mut,sigma0,xi,theta,eta,log=FALSE,tau=0,eta2=NULL) ...@@ -496,10 +498,10 @@ gev.d.lik <- function(xdat,ds,mut,sigma0,xi,theta,eta,log=FALSE,tau=0,eta2=NULL)
#' gev.d.diag(fit,subset = example$cov1==1,pch=1) #' gev.d.diag(fit,subset = example$cov1==1,pch=1)
gev.d.diag <- function(fit,subset=NULL,cols=NULL,pch=NULL,which='both',mfrow=c(1,2),legend=TRUE, 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'), title=c('Residual Probability Plot','Residual Quantile Plot'),
emp.lab='Empirical',mod.lab='Model',...){ emp.lab='Empirical',mod.lab='Model',ci=FALSE,...){
# check parameter: # check parameter:
if(!is.element(which,c('both','pp','qq'))) stop("Parameter 'which'= ",which, if(!is.element(which,c('both','pp','qq'))) stop("Parameter 'which'= ",which,
" but only 'both','pp' or 'qq' are allowed.") " but only 'both','pp' or 'qq' are allowed.")
# subset data # subset data
df <- data.frame(data=fit$data,ds=fit$ds) df <- data.frame(data=fit$data,ds=fit$ds)
if(!is.null(subset)){ if(!is.null(subset)){
...@@ -511,34 +513,50 @@ gev.d.diag <- function(fit,subset=NULL,cols=NULL,pch=NULL,which='both',mfrow=c(1 ...@@ -511,34 +513,50 @@ gev.d.diag <- function(fit,subset=NULL,cols=NULL,pch=NULL,which='both',mfrow=c(1
durs <- sort(unique(df$ds)) durs <- sort(unique(df$ds))
# rescale durations to assign colors # rescale durations to assign colors
df$cval <- sapply(df$ds,function(d){which(durs==d)}) df$cval <- sapply(df$ds,function(d){which(durs==d)})
# sort data # sort data
df <- df[order(df$data),] df <- df[order(df$data),]
# plotting position # plotting position
n <- length(df$data) n <- length(df$data)
px <- (1:n)/(n + 1) px <- (1:n)/(n + 1)
# get 95% confidence intervals
if(ci){
samp <- rgev(n * 99, loc = 0, scale = 1, shape = 0)
samp <- matrix(samp, n, 99)
samp <- apply(samp, 2, sort)
samp <- apply(samp, 1, sort)
ci.vals <- t(samp[c(3, 97), ])
}else{ci.vals <- NULL}
# create plots: # create plots:
if(which=='both') par(mfrow=mfrow) # 2 subplots if(which=='both') par(mfrow=mfrow) # 2 subplots
# colors and symbols # colors and symbols
if(is.null(cols))cols <- rainbow(length(durs)) if(is.null(cols))cols <- rainbow(length(durs))
if(is.null(pch))pch <- df$cval if(is.null(pch))pch <- df$cval
if(which=='both'|which=='pp'){ if(which=='both'|which=='pp'){
# pp # pp
plot(px, exp( - exp( - df$data)), xlab = plot(px, exp( - exp( - df$data)), xlab = emp.lab, ylab = mod.lab,col=cols[df$cval]
emp.lab, ylab = mod.lab,col=cols[df$cval],pch=pch,...) ,pch=pch,ylim=range(exp( - exp(-c(ci.vals,df$data))),na.rm = TRUE),...)
abline(0, 1, col = 1,lwd=1) abline(0, 1, col = 1,lwd=1)
if(ci){
lines(px,exp( - exp( - ci.vals[,1])),lty=2)
lines(px,exp( - exp( - ci.vals[,2])),lty=2)
}
title(title[1]) title(title[1])
if(legend){legend('bottomright',legend = round(durs,digits = 2),pch=pch, if(legend){legend('bottomright',legend = round(durs,digits = 2),pch=pch,
col = cols[1:length(durs)],title = 'Duration [h]',ncol = 2)} col = cols[1:length(durs)],title = 'Duration [h]',ncol = 2)}
} }
if(which=='both'|which=='qq'){ if(which=='both'|which=='qq'){
# qq # qq
plot( - log( - log(px)), df$data, ylab = plot( - log( - log(px)), df$data, ylab = emp.lab, xlab = mod.lab,col=cols[df$cval]
emp.lab, xlab = mod.lab,col=cols[df$cval],pch=pch,...) ,pch=pch,ylim=range(c(ci.vals,df$data),na.rm = TRUE),...)
abline(0, 1, col = 1,lwd=1) abline(0, 1, col = 1,lwd=1)
if(ci){
lines(-log(-log(px)),ci.vals[,1],lty=2)
lines(-log(-log(px)),ci.vals[,2],lty=2)
}
title(title[2]) title(title[2])
if(legend){legend('bottomright',legend = round(durs,digits = 2),pch=pch, if(legend){legend('bottomright',legend = round(durs,digits = 2),pch=pch,
col = cols[1:length(durs)],title = 'Duration [h]',ncol = 2)} col = cols[1:length(durs)],title = 'Duration [h]',ncol = 2)}
...@@ -546,6 +564,7 @@ gev.d.diag <- function(fit,subset=NULL,cols=NULL,pch=NULL,which='both',mfrow=c(1 ...@@ -546,6 +564,7 @@ gev.d.diag <- function(fit,subset=NULL,cols=NULL,pch=NULL,which='both',mfrow=c(1
if(which=='both') par(mfrow=c(1,1)) # reset par if(which=='both') par(mfrow=c(1,1)) # reset par
} }
#### gev.d.params #### #### gev.d.params ####
#' Calculate gev(d) parameters from \code{gev.d.fit} output #' Calculate gev(d) parameters from \code{gev.d.fit} output
......
...@@ -15,6 +15,7 @@ gev.d.diag( ...@@ -15,6 +15,7 @@ gev.d.diag(
title = c("Residual Probability Plot", "Residual Quantile Plot"), title = c("Residual Probability Plot", "Residual Quantile Plot"),
emp.lab = "Empirical", emp.lab = "Empirical",
mod.lab = "Model", mod.lab = "Model",
ci = FALSE,
... ...
) )
} }
...@@ -41,6 +42,8 @@ set to \code{c(1,1)}.} ...@@ -41,6 +42,8 @@ set to \code{c(1,1)}.}
\item{emp.lab, mod.lab}{character string containing names for empirical and model axis} \item{emp.lab, mod.lab}{character string containing names for empirical and model axis}
\item{ci}{logical indicating whether 0.95 confidence intervals should be plotted}
\item{...}{additional parameters passed on to the plotting function} \item{...}{additional parameters passed on to the plotting function}
} }
\description{ \description{
......
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