Nothing
#' Function for heuristically choosing the domain where extreme value theory can be applied
#'
#' @param data a numeric vector containing the observation used for verification
#' @param thresh vector of thresholds to try
#' @param guess starting values for GPD's sigma and xi (0<xi<1)
#' @param plots which parameter plots do you want
#' @param R number of bootstrap estimates for confidence intervals
#' @param ncpus if you want to make bootstrap on several cores
#' @return three plots summarizing the stability of the parameters to threshold. The starting threshold admits kappa=1 and its confidence interval ;
#' according Papastathopoulos & Tawn (2013)
#' @return a list with thresholds used, GP parameters and CIs, optimal threshold and xi.
#' @export
#'
#' @importFrom graphics abline arrows legend lines par plot points polygon title
choosethres=function(data,thresh,guess=c(1,0.1),plots=1:3,R=200,ncpus=1){
thresh=sort(thresh)
pe=matrix(0,ncol=3,nrow=length(thresh))
se=array(0,dim=c(length(thresh),3,2))
fi=suppressWarnings(.gp.pw.fit(data[data>thresh[1]]-thresh[1],init=c(1,guess),CI=T,R=R,ncpus=ncpus))
pe[1,]=fi$fit$PWM
se[1,,]=t(fi$CI$PWM)
for (i in 2:length(thresh)){
print(i)
fi=suppressWarnings(.gp.pw.fit(data[data>thresh[i]]-thresh[i],init=c(1,guess),CI=T,R=R,ncpus=ncpus))
pe[i,]=fi$fit$PWM
se[i,,]=t(fi$CI$PWM)
}
iopt=which((se[,1,1]<1)&(se[,1,2]>1))[1]
if (is.na(iopt)) {
warning("A pareto approximation of data is not available here")
}
if (!(length(plots) == 1 && is.na(plots))) {
plots <- sort(unique(plots))
if (!isTRUE(all(plots %in% 1:3))) {
stop("Invalid plot selection. Must be a vector of integers containing indices 1, 2 or 3.")
}
old.par <- par(no.readonly = TRUE)
par(mfrow = c(length(plots), 1), mar = c(4.5, 4.5, 3.1,
0.1))
on.exit(par(old.par))
for (i in plots) {
ylims = c(min(se[, i,],pe[,i]), max(se[, i,],pe[,i]))
plot(x = thresh, y = pe[, i], pch = 20, xlab = "Threshold",
bty = "l", ylab = switch(i, expression(kappa),
expression(sigma), expression(xi)),
ylim = ylims, type = "n")
polygon(c(thresh, rev(thresh)), c( se[, i,1], rev(se[, i,2])),
col = "gray57", border = FALSE)
abline(v=thresh[iopt],lty=2,col="seagreen")
if (i == min(plots)) {
title(paste0("Parameter stability plots for EGP"))
}
if (i == 1) {
abline(h = 1, lwd = 0.5, col = "gray20", lty = 2)
}
arrows(x0 = thresh, y0 = se[, i,1], y1 = se[, i,2],
length = 0.05, angle = 90, code = 3)
points(x = thresh, y = pe[, i], type = "b", pch = 20)
}
}
return(invisible(list(par = pe, CI = se, thresh = thresh, opt=c(thresh[iopt],pe[iopt,3]))))
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.