[Robast-commits] r1016 - in branches/robast-1.1/pkg/RobExtremes: . R inst/scripts man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sun Jul 22 13:33:48 CEST 2018
Author: ruckdeschel
Date: 2018-07-22 13:33:47 +0200 (Sun, 22 Jul 2018)
New Revision: 1016
Added:
branches/robast-1.1/pkg/RobExtremes/R/checkIC.R
branches/robast-1.1/pkg/RobExtremes/R/makeIC.R
branches/robast-1.1/pkg/RobExtremes/inst/scripts/RobFitsAtRealData.R
branches/robast-1.1/pkg/RobExtremes/man/checkmakeIC-methods.Rd
Modified:
branches/robast-1.1/pkg/RobExtremes/NAMESPACE
branches/robast-1.1/pkg/RobExtremes/R/GEVFamily.R
branches/robast-1.1/pkg/RobExtremes/R/GEVFamilyMuUnknown.R
branches/robast-1.1/pkg/RobExtremes/R/GParetoFamily.R
branches/robast-1.1/pkg/RobExtremes/R/getStartIC.R
branches/robast-1.1/pkg/RobExtremes/R/getStartICPareto.R
branches/robast-1.1/pkg/RobExtremes/R/gevgpddiag.R
branches/robast-1.1/pkg/RobExtremes/R/internal-getpsi.R
branches/robast-1.1/pkg/RobExtremes/R/startEstGEV.R
branches/robast-1.1/pkg/RobExtremes/R/startEstGPD.R
branches/robast-1.1/pkg/RobExtremes/man/GEVFamily.Rd
branches/robast-1.1/pkg/RobExtremes/man/GEVFamilyMuUnknown.Rd
branches/robast-1.1/pkg/RobExtremes/man/GParetoFamily.Rd
branches/robast-1.1/pkg/RobExtremes/man/getStartIC-methods.Rd
branches/robast-1.1/pkg/RobExtremes/man/ismevgpdgevdiag-methods.Rd
Log:
[RobExtremes] branch 1.1
+ new checkIC and makeIC methods for GEV[U], GPD, and Pareto
+ new script RobFitsAtRealData
+ GEVFamily, GParetoFamily and GEVFamilyMuUnknown gain argument withMDE (by default TRUE) which controls usage of MDEs at finding startPars
+ due to speed gains moved back diagnostic plot examples out of \donttest
+ added MakeIC calls in getStartIC to enhance accuracy
+ in order to avoid time-costly double calls to MakeIC, the respective modifyIC-slots gain attribute "hasMakeICin.modifyIC" once MakeIC is in the slot
+ added different route to get hand on L2Fam in gevgpddiag.R to make it applicable to return values of roptest
+ speeded up startEstGEV and startEstGPD
Modified: branches/robast-1.1/pkg/RobExtremes/NAMESPACE
===================================================================
--- branches/robast-1.1/pkg/RobExtremes/NAMESPACE 2018-07-22 11:19:26 UTC (rev 1015)
+++ branches/robast-1.1/pkg/RobExtremes/NAMESPACE 2018-07-22 11:33:47 UTC (rev 1016)
@@ -47,6 +47,7 @@
exportMethods("modifyModel", "getStartIC")
exportMethods("moveL2Fam2RefParam",
"moveICBackFromRefParam")
+exportMethods("checkIC", "makeIC")
export("EULERMASCHERONICONSTANT","APERYCONSTANT")
export("getCVaR", "getVaR", "getEL")
export("Gumbel", "Pareto", "GPareto", "GEV")
Modified: branches/robast-1.1/pkg/RobExtremes/R/GEVFamily.R
===================================================================
--- branches/robast-1.1/pkg/RobExtremes/R/GEVFamily.R 2018-07-22 11:19:26 UTC (rev 1015)
+++ branches/robast-1.1/pkg/RobExtremes/R/GEVFamily.R 2018-07-22 11:33:47 UTC (rev 1016)
@@ -181,6 +181,7 @@
secLevel = 0.7,
withCentL2 = FALSE,
withL2derivDistr = FALSE,
+ withMDE = FALSE,
..ignoreTrafo = FALSE,
..withWarningGEV = TRUE){
theta <- c(loc, scale, shape)
@@ -281,7 +282,7 @@
# PF <- GEVFamily(loc = theta[1], scale = theta[2], shape = theta[3])
# e1 <- PickandsEstimator(x,ParamFamily=PF)
# e0 <- estimate(e1)
- e0 <- .getBetaXiGEV(x=x, mu=mu, xiGrid=.getXiGrid(), withPos=withPos)
+ e0 <- .getBetaXiGEV(x=x, mu=mu, xiGrid=.getXiGrid(), withPos=withPos, withMDE=withMDE)
}else{
if(is(start0Est,"function")){
e1 <- start0Est(x, ...)
Modified: branches/robast-1.1/pkg/RobExtremes/R/GEVFamilyMuUnknown.R
===================================================================
--- branches/robast-1.1/pkg/RobExtremes/R/GEVFamilyMuUnknown.R 2018-07-22 11:19:26 UTC (rev 1015)
+++ branches/robast-1.1/pkg/RobExtremes/R/GEVFamilyMuUnknown.R 2018-07-22 11:33:47 UTC (rev 1016)
@@ -145,6 +145,7 @@
secLevel = 0.7,
withCentL2 = FALSE,
withL2derivDistr = FALSE,
+ withMDE = FALSE,
..ignoreTrafo = FALSE,
..withWarningGEV = TRUE,
..name =""){
@@ -247,7 +248,7 @@
# PF <- GEVFamily(loc = theta[1], scale = theta[2], shape = theta[3])
# e1 <- PickandsEstimator(x,ParamFamily=PF)
# e0 <- estimate(e1)
- e0 <- .getMuBetaXiGEV(x=x, xiGrid=.getXiGrid(), withPos=withPos)
+ e0 <- .getMuBetaXiGEV(x=x, xiGrid=.getXiGrid(), withPos=withPos, withMDE=withMDE)
}else{
if(is(start0Est,"function")){
e1 <- start0Est(x, ...)
Modified: branches/robast-1.1/pkg/RobExtremes/R/GParetoFamily.R
===================================================================
--- branches/robast-1.1/pkg/RobExtremes/R/GParetoFamily.R 2018-07-22 11:19:26 UTC (rev 1015)
+++ branches/robast-1.1/pkg/RobExtremes/R/GParetoFamily.R 2018-07-22 11:33:47 UTC (rev 1016)
@@ -44,6 +44,7 @@
secLevel = 0.7,
withCentL2 = FALSE,
withL2derivDistr = FALSE,
+ withMDE = FALSE,
..ignoreTrafo = FALSE){
theta <- c(loc, scale, shape)
@@ -139,7 +140,7 @@
medkMADhybr(c(x), k=10, ParamFamily = PF,
q.lo = 1e-3, q.up = 15), silent =TRUE)
if(is(e1,"try-error")){ e0 <- .getBetaXiGPD(x=x, mu=tr,
- xiGrid=.getXiGrid(), withPos=withPos)
+ xiGrid=.getXiGrid(), withPos=withPos, withMDE=withMDE)
}else e0 <- estimate(e1)
}else{
if(is(start0Est,"function")){
Added: branches/robast-1.1/pkg/RobExtremes/R/checkIC.R
===================================================================
--- branches/robast-1.1/pkg/RobExtremes/R/checkIC.R (rev 0)
+++ branches/robast-1.1/pkg/RobExtremes/R/checkIC.R 2018-07-22 11:33:47 UTC (rev 1016)
@@ -0,0 +1,24 @@
+..checkIC.qtl <- function (IC, L2Fam, out = TRUE, ...){
+ mc <- match.call(expand.dots = TRUE)
+ mcF <- match.call(expand.dots = FALSE)
+ dots <- mcF$"..."
+ mcl <- as.list(mc)[-1]
+ mcl$out <- out
+ mcl$IC <- IC
+ D1 <- L2Fam at distribution
+ D1 <- as(D1,"DistributionsIntegratingByQuantiles")
+ L2Fam at distribution <- D1
+ L2Fam <- as(L2Fam, "L2ParamFamily")
+ mcl$L2Fam <- L2Fam
+ if(!is.null(dots)) mcl <- c(mcl,dots)
+ do.call("checkIC", mcl)
+ }
+setMethod("checkIC", signature(IC = "IC", L2Fam = "GParetoFamily"),
+ ..checkIC.qtl)
+
+setMethod("checkIC", signature(IC = "IC", L2Fam = "GEVFamilyMuUnknown"),
+ ..checkIC.qtl)
+setMethod("checkIC", signature(IC = "IC", L2Fam = "GEVFamily"),
+ ..checkIC.qtl)
+setMethod("checkIC", signature(IC = "IC", L2Fam = "ParetoFamily"),
+ ..checkIC.qtl)
Modified: branches/robast-1.1/pkg/RobExtremes/R/getStartIC.R
===================================================================
--- branches/robast-1.1/pkg/RobExtremes/R/getStartIC.R 2018-07-22 11:19:26 UTC (rev 1015)
+++ branches/robast-1.1/pkg/RobExtremes/R/getStartIC.R 2018-07-22 11:33:47 UTC (rev 1016)
@@ -8,6 +8,7 @@
gridn <- gsub("\\.","",type(risk))
nam <- paste(".",gsub("[F,f]amily","",gsub(" ","",name(model))),sep="")
+
if(nam==".GeneralizedPareto") nam <- ".GPareto"
param1 <- param(model)
@@ -28,15 +29,22 @@
para <- param(L2Fam)
if(!.is.na.Psi(para, interpolfct, shnam))
return(.getPsi(para, interpolfct, L2Fam, type(risk)))
- else
- return(do.call(getStartIC, as.list(mc[-1]),
- envir=parent.frame(2)))
+ else{
+ IC0 <- do.call(getStartIC, as.list(mc[-1]),
+ envir=parent.frame(2))
+ IC0 <- makeIC(IC0, L2Fam)
+ return(IC0)
+ }
}
+ attr(.modifyIC0,"hasMakeICin.modifyIC") <- TRUE
+
.modifyIC <- function(L2Fam,IC){
psi.0 <- .modifyIC0(L2Fam,IC)
psi.0 at modifyIC <- .modifyIC
return(psi.0)
}
+ attr(.modifyIC,"hasMakeICin.modifyIC") <- TRUE
+
if(!.is.na.Psi(param1, interpolfct, shnam)){
IC0 <- .getPsi(param1, interpolfct, model, type(risk))
IC0 at modifyIC <- .modifyIC
@@ -44,7 +52,9 @@
}
}
}
- return(do.call(getStartIC, as.list(mc[-1]), envir=parent.frame(2)))
+ IC <- do.call(getStartIC, as.list(mc[-1]), envir=parent.frame(2))
+ IC <- makeIC(IC,model)
+ return(IC)
})
setMethod("getStartIC",signature(model = "L2LocScaleShapeUnion", risk = "interpolRisk"),
@@ -74,15 +84,21 @@
para <- param(L2Fam)
if(!.is.na.Psi(para, interpolfct, shnam))
return(.getPsi.wL(para, interpolfct, L2Fam, type(risk)))
- else
- return(do.call(getStartIC, as.list(mc[-1]),
- envir=parent.frame(2)))
+ else{
+ IC0 <- do.call(getStartIC, as.list(mc[-1]),
+ envir=parent.frame(2))
+ IC0 <- makeIC(IC0, L2Fam)
+ return(IC0)
+ }
}
+ attr(.modifyIC0,"hasMakeICin.modifyIC") <- TRUE
.modifyIC <- function(L2Fam,IC){
psi.0 <- .modifyIC0(L2Fam,IC)
psi.0 at modifyIC <- .modifyIC
return(psi.0)
}
+ attr(.modifyIC,"hasMakeICin.modifyIC") <- TRUE
+
if(!.is.na.Psi(param1, interpolfct, shnam)){
IC0 <- .getPsi.wL(param1, interpolfct, model, type(risk))
IC0 at modifyIC <- .modifyIC
@@ -90,6 +106,8 @@
}
}
}
- return(do.call(getStartIC, as.list(mc[-1]), envir=parent.frame(2)))
+ IC <- do.call(getStartIC, as.list(mc[-1]), envir=parent.frame(2))
+ IC <- makeIC(IC,model)
+ return(IC)
})
Modified: branches/robast-1.1/pkg/RobExtremes/R/getStartICPareto.R
===================================================================
--- branches/robast-1.1/pkg/RobExtremes/R/getStartICPareto.R 2018-07-22 11:19:26 UTC (rev 1015)
+++ branches/robast-1.1/pkg/RobExtremes/R/getStartICPareto.R 2018-07-22 11:33:47 UTC (rev 1016)
@@ -7,11 +7,13 @@
xi0 <- main(param(L2Fam))
return(.getPsi.P(xi0, L2Fam, type(risk)))
}
+ attr(.modifyIC0,"hasMakeICin.modifyIC") <- TRUE
.modifyIC <- function(L2Fam,IC){
psi.0 <- .modifyIC0(L2Fam,IC)
psi.0 at modifyIC <- .modifyIC
return(psi.0)
}
+ attr(.modifyIC,"hasMakeICin.modifyIC") <- TRUE
IC0 <- .getPsi.P(xi, model, type(risk))
IC0 at modifyIC <- .modifyIC
return(IC0)
@@ -68,5 +70,6 @@
IC <- generateIC(nb, L2Fam, res)
+ IC <- makeIC(IC,L2Fam)
return(IC)
}
Modified: branches/robast-1.1/pkg/RobExtremes/R/gevgpddiag.R
===================================================================
--- branches/robast-1.1/pkg/RobExtremes/R/gevgpddiag.R 2018-07-22 11:19:26 UTC (rev 1015)
+++ branches/robast-1.1/pkg/RobExtremes/R/gevgpddiag.R 2018-07-22 11:33:47 UTC (rev 1016)
@@ -73,12 +73,16 @@
else param <- ParamFamParameter(main=utfe, nuisance=nuisance(z))
es.call <- z at estimate.call
nm.call <- names(es.call)
- PFam <- NULL
- if("ParamFamily" %in% nm.call)
- PFam <- eval(as.list(es.call)[["ParamFamily"]])
- if(is.null(PFam))
- stop("There is no object of class 'ProbFamily' in the call of 'z'")
- PFam0 <- modifyModel(PFam, param)
+ if("pIC" %in% names(getSlots(class(z)))){
+ PFam0 <- eval(z at pIC@CallL2Fam)
+ }else{
+ PFam <- NULL
+ if("ParamFamily" %in% nm.call)
+ PFam <- eval(as.list(es.call)[["ParamFamily"]])
+ if(is.null(PFam))
+ stop("There is no object of class 'ProbFamily' in the call of 'z'")
+ PFam0 <- modifyModel(PFam, param)
+ }
thresh <- if(GPD) fixed(param) else NULL
x <- eval(es.call$x)
n0 <- length(x)
Modified: branches/robast-1.1/pkg/RobExtremes/R/internal-getpsi.R
===================================================================
--- branches/robast-1.1/pkg/RobExtremes/R/internal-getpsi.R 2018-07-22 11:19:26 UTC (rev 1015)
+++ branches/robast-1.1/pkg/RobExtremes/R/internal-getpsi.R 2018-07-22 11:33:47 UTC (rev 1016)
@@ -52,6 +52,7 @@
IC <- generateIC(nb, L2Fam, res)
+ IC <- makeIC(IC,L2Fam)
return(IC)
}
@@ -108,6 +109,7 @@
IC <- generateIC(nb, L2Fam, res)
+ IC <- makeIC(IC,L2Fam)
return(IC)
}
Added: branches/robast-1.1/pkg/RobExtremes/R/makeIC.R
===================================================================
--- branches/robast-1.1/pkg/RobExtremes/R/makeIC.R (rev 0)
+++ branches/robast-1.1/pkg/RobExtremes/R/makeIC.R 2018-07-22 11:33:47 UTC (rev 1016)
@@ -0,0 +1,20 @@
+..makeIC.qtl <- function (IC, L2Fam){
+ mc <- match.call()
+ mcl <- as.list(mc)[-1]
+ mcl$IC <- IC
+ D1 <- L2Fam at distribution
+ D1 <- as(D1,"DistributionsIntegratingByQuantiles")
+ L2Fam at distribution <- D1
+ L2Fam <- as(L2Fam, "L2ParamFamily")
+ mcl$L2Fam <- L2Fam
+ do.call("makeIC", mcl)
+ }
+setMethod("makeIC", signature(IC = "IC", L2Fam = "GParetoFamily"),
+ ..makeIC.qtl)
+
+setMethod("makeIC", signature(IC = "IC", L2Fam = "GEVFamilyMuUnknown"),
+ ..makeIC.qtl)
+setMethod("makeIC", signature(IC = "IC", L2Fam = "GEVFamily"),
+ ..makeIC.qtl)
+setMethod("makeIC", signature(IC = "IC", L2Fam = "ParetoFamily"),
+ ..makeIC.qtl)
Modified: branches/robast-1.1/pkg/RobExtremes/R/startEstGEV.R
===================================================================
--- branches/robast-1.1/pkg/RobExtremes/R/startEstGEV.R 2018-07-22 11:19:26 UTC (rev 1015)
+++ branches/robast-1.1/pkg/RobExtremes/R/startEstGEV.R 2018-07-22 11:33:47 UTC (rev 1016)
@@ -1,8 +1,8 @@
-.getXiGrid <- function(){c(-0.48,-0.1,0,0.1,0.4,1,3,6)}
+.getXiGrid <- function(){c(1, -0.48,0,3,-0.1,0.4 ,0.1,6)}
.getBetaXiGEV <- function(x, mu, xiGrid = .getXiGrid(), withPos=TRUE, secLevel = 0.7,
- .issueIntermediateParams = FALSE){
+ .issueIntermediateParams = FALSE, withMDE = FALSE){
n <- length(x)
epsn <- min(floor(secLevel*sqrt(n))+1,n)
@@ -11,23 +11,28 @@
s0 <- max(x0)-min(x0)
crit0 <- Inf
- fu <- function(x,...) .getBetaXiGEV(x,mu,xiGrid = xiGrid,withPos=withPos)
+ fu <- function(x,...) .getBetaXiGEV(x,0,xiGrid = xiGrid,withPos=withPos)
e0 <- NULL
es <- c(NA,NA)
### first try (to ensure global consistency): PickandsEstimator
try({mygev <- GEVFamily(loc=0,scale=1,shape=0.1, withPos=withPos,
..withWarningGEV=FALSE)
- e1 <- PickandsEstimator(x,ParamFamily=mygev)
+ e1 <- PickandsEstimator(x0,ParamFamily=mygev)
if(.issueIntermediateParams){
cat("Pickands:\n");print(e1) }
e0 <- estimate(e1)}, silent=TRUE)
+ validi <- 0
+ es0 <- c(NA,NA)
if(!is.null(e0)) if(!is(e0,"try-error")){
+ if(!withMDE) {
+ names(e0) <- c("scale","shape")
+ return(e0)
+ }
mygev <- GEVFamily(loc=0,scale=e0[1],shape=e0[2], withPos=withPos,
start0Est = fu, ..withWarningGEV=FALSE)
mde0 <- try(MDEstimator(x0, mygev, distance=CvMDist, startPar=c("scale"=e0[1],"shape"=e0[2])),silent=TRUE)
- es0 <- c(NA,NA)
if(!is(mde0,"try-error")){
es <- estimate(mde0)
crit1 <- criterion(mde0)
@@ -35,12 +40,15 @@
cat("1st candidate:\n", round(es,6), " crit:", round(crit1,6), , " ")
}
if(quantile(1+es[2]*x0/es[1], epsn/n)>0){
+ validi <- 1
mdeb <- mde0
crit0 <- crit1
es0 <- es
if(.issueIntermediateParams){
cat("side condition '1+sc/sh (x-mu) >0' fulfilled;\n")
}
+ names(es) <- c("scale","shape")
+ return(es)
}else{
if(.issueIntermediateParams){
cat("side condition '1+sc/sh (x-mu) >0' violated;\n")
@@ -50,13 +58,18 @@
}
i <- 0
- for(xi in xiGrid){
+ sd0 <- c(Inf,Inf)
+ esS <- matrix(NA,length(xiGrid)+validi,2)
+ if(validi>0) esS[1,] <- es0
+
+ while((i<length(xiGrid))&&max(sd0)>1e-3){
i <- i + 1
+ xi <- xiGrid[i]
funl <- function(sig){
mygev1 <- GEV(loc=0,scale=sig,shape=xi)
CvMDist(x0,mygev1)
}
- intlo <- quantile(-xi*(x-mu),1-epsn/n)
+ intlo <- quantile(-xi*x0,1-epsn/n)
intv <- c(max(1e-5,intlo), s0)
sigCvMMD1 <- optimize(funl, interval=intv)$minimum
mygev <- GEVFamily(loc=0,scale=sigCvMMD1,shape=xi, withPos=withPos,
@@ -70,6 +83,9 @@
cat("candidate no",i+1, ":\n", round(es,6), " crit:", round(crit1,6), " ")
}
if(quantile(1+es[2]*x0/es[1], epsn/n)>0){
+ validi <- validi+1
+ esS[validi,] <- es
+ if(validi>2) sd0 <- apply(esS[1:validi,],2,sd)
if(.issueIntermediateParams){
cat("side condition '1+sc/sh (x-mu) >0' fulfilled;\n")
}
@@ -85,6 +101,9 @@
es[1] <- intlo+1e-5
mygev2 <- GEV(loc=0,scale=es[1],shape=es[2])
crit1 <- CvMDist(x0,mygev2)
+ validi <- validi+1
+ esS[validi,] <- es
+ if(validi>2) sd0 <- apply(esS[1:validi,],2,sd)
if(.issueIntermediateParams){
cat("candidate no",i+1, "(b):\n", round(es,6), " crit:", round(crit1,6), " ")
}
@@ -100,11 +119,12 @@
}
.getMuBetaXiGEV <- function(x, xiGrid = .getXiGrid(), withPos=TRUE, secLevel = 0.7,
- .issueIntermediateParams = FALSE){
+ .issueIntermediateParams = FALSE, withMDE = FALSE){
mu <- quantile(x,exp(-1))
es <- .getBetaXiGEV(x=x, mu=mu, xiGrid=xiGrid, withPos=withPos,
secLevel = secLevel,
- .issueIntermediateParams = .issueIntermediateParams)
+ .issueIntermediateParams = .issueIntermediateParams,
+ withMDE = withMDE)
es0 <- c(mu,es)
names(es0) <- c("loc","scale","shape")
return(es0)
Modified: branches/robast-1.1/pkg/RobExtremes/R/startEstGPD.R
===================================================================
--- branches/robast-1.1/pkg/RobExtremes/R/startEstGPD.R 2018-07-22 11:19:26 UTC (rev 1015)
+++ branches/robast-1.1/pkg/RobExtremes/R/startEstGPD.R 2018-07-22 11:33:47 UTC (rev 1016)
@@ -1,5 +1,5 @@
.getBetaXiGPD <- function(x, mu, xiGrid = .getXiGrid(), withPos=TRUE, secLevel = 0.7,
- .issueIntermediateParams = FALSE){
+ .issueIntermediateParams = FALSE, withMDE = FALSE){
n <- length(x)
epsn <- min(floor(secLevel*sqrt(n))+1,n)
@@ -8,22 +8,27 @@
s0 <- max(x0)-min(x0)
crit0 <- Inf
- fu <- function(x,...) .getBetaXiGPD(x,mu,xiGrid = xiGrid,withPos=withPos)
+ fu <- function(x,...) .getBetaXiGPD(x,0,xiGrid = xiGrid,withPos=withPos)
e0 <- NULL
es <- c(NA,NA)
### first try (to ensure global consistency): PickandsEstimator
try({mygpd <- GParetoFamily(loc=0,scale=1,shape=0.1, withPos=withPos)
- e1 <- medkMADhybr(x,ParamFamily=mygpd, k=10)
+ e1 <- medkMADhybr(x0,ParamFamily=mygpd, k=10)
if(.issueIntermediateParams){
cat("MedkMAD:\n");print(e1)}
e0 <- estimate(e1)}, silent=TRUE)
+ validi <- 0
+ es0 <- c(NA,NA)
if(!is.null(e0)) if(!is(e0,"try-error")){
+ if(!withMDE) {
+ names(e0) <- c("scale","shape")
+ return(e0)
+ }
mygpd <- GParetoFamily(loc=0,scale=e0[1],shape=e0[2], withPos=withPos,
start0Est = fu)
mde0 <- try(MDEstimator(x0, mygpd, distance=CvMDist, startPar=c("scale"=es0[1],"shape"=es0[2])),silent=TRUE)
- es0 <- c(NA,NA)
if(!is(mde0,"try-error")){
es <- estimate(mde0)
crit1 <- criterion(mde0)
@@ -31,12 +36,15 @@
cat("1st candidate:\n", round(es,6), " crit:", round(crit1,6), , " ")
}
if(quantile(1+es[2]*x0/es[1], epsn/n)>0){
+ validi <- 1
mdeb <- mde0
crit0 <- crit1
es0 <- es
if(.issueIntermediateParams){
cat("side condition '1+sc/sh (x-mu) >0' fulfilled;\n")
}
+ names(es) <- c("scale","shape")
+ return(es)
}else{
if(.issueIntermediateParams){
cat("side condition '1+sc/sh (x-mu) >0' violated;\n")
@@ -46,8 +54,13 @@
}
i <- 0
- for(xi in xiGrid){
+ sd0 <- c(Inf,Inf)
+ esS <- matrix(NA,length(xiGrid)+validi,2)
+ if(validi>0) esS[1,] <- es0
+
+ while((i<length(xiGrid))&&max(sd0)>1e-3){
i <- i + 1
+ xi <- xiGrid[i]
funl <- function(sig){
mygpd1 <- GPareto(loc=0,scale=sig,shape=xi)
CvMDist(x0,mygpd1)
@@ -66,6 +79,9 @@
cat("candidate no",i+1, ":\n", round(es,6), " crit:", round(crit1,6), " ")
}
if(quantile(1+es[2]*x0/es[1], epsn/n)>0){
+ validi <- validi+1
+ esS[validi,] <- es
+ if(validi>2) sd0 <- apply(esS[1:validi,],2,sd)
if(.issueIntermediateParams){
cat("side condition '1+sc/sh (x-mu) >0' fulfilled;\n")
}
@@ -81,6 +97,9 @@
es[1] <- intlo+1e-5
mygpd2 <- GPareto(loc=0,scale=es[1],shape=es[2])
crit1 <- CvMDist(x0,mygpd2)
+ validi <- validi+1
+ esS[validi,] <- es
+ if(validi>2) sd0 <- apply(esS[1:validi,],2,sd)
if(.issueIntermediateParams){
cat("candidate no",i+1, "(b):\n", round(es,6), " crit:", round(crit1,6), " ")
}
Added: branches/robast-1.1/pkg/RobExtremes/inst/scripts/RobFitsAtRealData.R
===================================================================
--- branches/robast-1.1/pkg/RobExtremes/inst/scripts/RobFitsAtRealData.R (rev 0)
+++ branches/robast-1.1/pkg/RobExtremes/inst/scripts/RobFitsAtRealData.R 2018-07-22 11:33:47 UTC (rev 1016)
@@ -0,0 +1,118 @@
+#######################################################################
+# example of robust fits at real data not included in Rd due to timings
+#######################################################################
+
+require(RobExtremes)
+require(ismev)
+data(portpirie)
+data(rain)
+detach(package:ismev)
+raini <- rain[rain>10]
+rainc <- c(raini,1000,10000)
+portpiriei <- portpirie[,2]
+portpiriec <- c(portpiriei,100)
+
+##--------------------------
+## GEV
+##--------------------------
+
+ppfiti <- ismev::gev.fit(portpiriei)
+gev.diag(ppfiti)
+ ##
+mlEi <- MLEstimator(portpiriei, GEVFamilyMuUnknown(withPos=FALSE))
+system.time(MBRi <- roptest(portpiriei, GEVFamilyMuUnknown(withPos=FALSE),risk=MBRRisk()))
+mlEi
+MBRi
+estimate(mlEi)
+estimate(MBRi)
+gev.diag(mlEi)
+gev.diag(MBRi)
+gev.prof(mlEi, m = 10, 4.1, 5)
+gev.profxi(mBRi, -0.3, 0.3)
+plot(MBRi at pIC)
+
+## contaminated:
+ppfitc <- ismev::gev.fit(portpiriec)
+gev.diag(ppfitc)
+ ##
+mlEc <- MLEstimator(portpiriec, GEVFamilyMuUnknown(withPos=FALSE))
+system.time(MBRc <- roptest(portpiriec, GEVFamilyMuUnknown(withPos=FALSE),risk=MBRRisk()))
+mlEc
+MBRc
+estimate(mlEc)
+estimate(MBRc)
+estimate(mlEi)
+estimate(MBRi)
+gev.diag(mlEc)
+gev.diag(MBRc)
+gev.prof(mlEc, m = 10, 4.1, 5)
+gev.profxi(mlEc, -0.3, 0.3)
+qqplot(portpiriec,MBRc)
+qqplot(portpiriec,MBRc,ylim=c(3.5,5))
+returnlevelplot(portpiriec,MBRc)
+returnlevelplot(portpiriec,MBRc,ylim=c(3.5,5))
+
+plot(MBRc at pIC)
+
+##--------------------------
+## GPD
+##--------------------------
+
+rnfiti <- ismev::gpd.fit(rain,10)
+gpd.diag(rnfiti)
+ ##
+mlE2i <- MLEstimator(raini, GParetoFamily(loc=10))
+gpd.diag(mlE2i)
+system.time(MBR2i <- roptest(raini, GParetoFamily(loc=10),risk=MBRRisk()))
+mlE2i
+MBR2i
+estimate(mlE2i)
+estimate(MBR2i)
+plot(MBR2i at pIC)
+gpd.diag(mlE2i)
+gpd.diag(MBR2i)
+gpd.prof(mlE2i, m = 10, 55, 77)
+gpd.profxi(MBR2i, -0.02, 0.02)
+GP <- GParetoFamily(scale=rnfiti$mle[1],shape=rnfiti$mle[2],loc=10)
+returnlevelplot(rain, GP, MaxOrPOT="POT", xlim=c(1e-1,1e3))
+
+## contaminated:
+
+rnfitc <- ismev::gpd.fit(c(rain,1000,10000),10)
+gpd.diag(rnfitc)
+ ##
+mlE2c <- MLEstimator(rainc, GParetoFamily(loc=10))
+gpd.diag(mlE2c)
+system.time(MBR2c <- roptest(rainc, GParetoFamily(loc=10),risk=MBRRisk()))
+mlE2c
+MBR2c
+estimate(mlE2c)
+estimate(MBR2c)
+estimate(mlE2i)
+estimate(MBR2i)
+gpd.diag(mlE2c)
+gpd.diag(MBR2c)
+gpd.prof(mlE2c, m = 10, 55, 77)
+gpd.profxi(mlE2c, -0.02, 0.02)
+plot(MBR2c at pIC)
+## to be fixed
+qqplot(rainc,MBR2c)
+qqplot(rainc,MBR2c,ylim=c(5,100))
+qqplot(rainc,MBR2c,xlim=c(5,100),ylim=c(5,100),log="xy")
+## to be fixed
+returnlevelplot(raini,MBR2i,MaxOrPot="POT",threshold=0)
+returnlevelplot(rainc,MBR2c,MaxOrPot="POT",threshold=0)
+returnlevelplot(rainc,MBR2c,ylim=c(10,100),MaxOrPot="POT",threshold=0)
+#
+L2F <- eval(MBR2c at pIC@CallL2Fam)
+dI2c <- L2F at distribution
+#loc(dI2c) <- 0
+qqplot(rainc,dI2c)
+rainc.10 <- rainc-10
+qqplot(rainc.10,dI2c-10)
+## to be fixed
+returnlevelplot(rainc.10,dI2c-10,MaxOrPot="POT",threshold=0)
+dI2i <- distribution(eval(MBR2i at pIC@CallL2Fam))
+loc(dI2i) <- 0
+qqplot(portpiriei-10,dI2i)
+qqplot(portpiriec,MBR2c)
Modified: branches/robast-1.1/pkg/RobExtremes/man/GEVFamily.Rd
===================================================================
--- branches/robast-1.1/pkg/RobExtremes/man/GEVFamily.Rd 2018-07-22 11:19:26 UTC (rev 1015)
+++ branches/robast-1.1/pkg/RobExtremes/man/GEVFamily.Rd 2018-07-22 11:33:47 UTC (rev 1016)
@@ -10,7 +10,7 @@
GEVFamily(loc = 0, scale = 1, shape = 0.5, of.interest = c("scale", "shape"),
p = NULL, N = NULL, trafo = NULL, start0Est = NULL, withPos = TRUE,
secLevel = 0.7, withCentL2 = FALSE, withL2derivDistr = FALSE,
- ..ignoreTrafo = FALSE, ..withWarningGEV = TRUE)
+ withMDE = FALSE, ..ignoreTrafo = FALSE, ..withWarningGEV = TRUE)
}
\arguments{
\item{loc}{ real: known/fixed threshold/location parameter }
@@ -35,7 +35,10 @@
the E()? Defaults to \code{FALSE}, but higher accuracy can be achieved
when set to \code{TRUE}.}
\item{withL2derivDistr}{logical: shall the distribution of the L2 derivative
- be computed? Defaults to \code{FALSE} (to speeds up computations).}
+ be computed? Defaults to \code{FALSE} (to speed up computations).}
+ \item{withMDE}{logical: should Minimum Distance Estimators be used to
+ find a good starting value for the parameter search?
+ Defaults to \code{FALSE} (to speed up computations).}
\item{..ignoreTrafo}{logical: only used internally in \code{kStepEstimator}; do not change this.}
\item{..withWarningGEV}{logical: shall warnings be issued if shape is large?}
}
Modified: branches/robast-1.1/pkg/RobExtremes/man/GEVFamilyMuUnknown.Rd
===================================================================
--- branches/robast-1.1/pkg/RobExtremes/man/GEVFamilyMuUnknown.Rd 2018-07-22 11:19:26 UTC (rev 1015)
+++ branches/robast-1.1/pkg/RobExtremes/man/GEVFamilyMuUnknown.Rd 2018-07-22 11:33:47 UTC (rev 1016)
@@ -10,7 +10,7 @@
GEVFamilyMuUnknown(loc = 0, scale = 1, shape = 0.5, of.interest = c("loc",
"scale", "shape"), p = NULL, N = NULL, trafo = NULL,
start0Est = NULL, withPos = TRUE, secLevel = 0.7,
- withCentL2 = FALSE, withL2derivDistr = FALSE,
+ withCentL2 = FALSE, withL2derivDistr = FALSE, withMDE = FALSE,
..ignoreTrafo = FALSE, ..withWarningGEV = TRUE, ..name = "")
}
\arguments{
@@ -36,7 +36,10 @@
the E()? Defaults to \code{FALSE}, but higher accuracy can be achieved
when set to \code{TRUE}.}
\item{withL2derivDistr}{logical: shall the distribution of the L2 derivative
- be computed? Defaults to \code{FALSE} (to speeds up computations).}
+ be computed? Defaults to \code{FALSE} (to speed up computations).}
+ \item{withMDE}{logical: should Minimum Distance Estimators be used to
+ find a good starting value for the parameter search?
+ Defaults to \code{FALSE} (to speed up computations).}
\item{..ignoreTrafo}{logical: only used internally in \code{kStepEstimator}; do not change this.}
\item{..withWarningGEV}{logical: shall warnings be issued if shape is large?}
\item{..name}{character: optional alternative name for the parametric family;
Modified: branches/robast-1.1/pkg/RobExtremes/man/GParetoFamily.Rd
===================================================================
--- branches/robast-1.1/pkg/RobExtremes/man/GParetoFamily.Rd 2018-07-22 11:19:26 UTC (rev 1015)
+++ branches/robast-1.1/pkg/RobExtremes/man/GParetoFamily.Rd 2018-07-22 11:33:47 UTC (rev 1016)
@@ -10,7 +10,7 @@
GParetoFamily(loc = 0, scale = 1, shape = 0.5, of.interest = c("scale", "shape"),
p = NULL, N = NULL, trafo = NULL, start0Est = NULL, withPos = TRUE,
secLevel = 0.7, withCentL2 = FALSE, withL2derivDistr = FALSE,
- ..ignoreTrafo = FALSE)
+ withMDE = FALSE, ..ignoreTrafo = FALSE)
}
\arguments{
\item{loc}{ real: known/fixed threshold/location parameter }
@@ -35,7 +35,10 @@
the E()? Defaults to \code{FALSE}, but higher accuracy can be achieved
when set to \code{TRUE}.}
\item{withL2derivDistr}{logical: shall the distribution of the L2 derivative
- be computed? Defaults to \code{FALSE} (to speeds up computations).}
+ be computed? Defaults to \code{FALSE} (to speed up computations).}
+ \item{withMDE}{logical: should Minimum Distance Estimators be used to
+ find a good starting value for the parameter search?
+ Defaults to \code{FALSE} (to speed up computations).}
\item{..ignoreTrafo}{logical: only used internally in \code{kStepEstimator}; do not change this.}
}
\details{
Added: branches/robast-1.1/pkg/RobExtremes/man/checkmakeIC-methods.Rd
===================================================================
--- branches/robast-1.1/pkg/RobExtremes/man/checkmakeIC-methods.Rd (rev 0)
+++ branches/robast-1.1/pkg/RobExtremes/man/checkmakeIC-methods.Rd 2018-07-22 11:33:47 UTC (rev 1016)
@@ -0,0 +1,58 @@
+\name{checkmakeIC-methods}
+\docType{methods}
+\alias{checkIC,IC,ParetoFamily-method}
+\alias{checkIC,IC,GParetoFamily-method}
+\alias{checkIC,IC,GEVFamily-method}
+\alias{checkIC,IC,GEVFamilyMuUnknown-method}
+\alias{makeIC,IC,ParetoFamily-method}
+\alias{makeIC,IC,GParetoFamily-method}
+\alias{makeIC,IC,GEVFamily-method}
+\alias{makeIC,IC,GEVFamilyMuUnknown-method}
+
+\title{Methods for Functions checkIC and makeIC in Package `RobExtremes' }
+
+\description{\code{checkIC} checks accuracy of the centering
+and Fisher consistency condition of an IC, \code{makeIC},
+by centering and restandardizing warrants these conditions.}
+
+\section{Methods}{\describe{
+\item{checkIC}{\code{signature(IC="IC", L2Fam = "ParetoFamily")}:
+ To enhance accuracy, the method for \code{"ParetoFamily"} uses
+ integration via the quantile transform, i.e., \eqn{E[h(X)]}
+ for a random variable \eqn{X\sim F}{X~F} with quantil function \eqn{q}
+ is computed as \eqn{\int_0^1 h(q(s))\,ds}{integral(from=0, to=1, integrand=h(q(s)))}
+}
+\item{checkIC}{\code{signature(IC="IC", L2Fam = "GParetoFamily")}:
+ As for \code{"ParetoFamily"}, to enhance accuracy,
+ the method for \code{"GParetoFamily"} uses
+ integration via the quantile transform.}
+\item{checkIC}{\code{signature(IC="IC", L2Fam = "GEVFamily")}:
+ As for \code{"ParetoFamily"}, to enhance accuracy,
+ the method for \code{"GEVFamily"} uses
+ integration via the quantile transform.}
+\item{checkIC}{\code{signature(IC="IC", L2Fam = "GEVFamilyMuUnknown")}:
+ As for \code{"ParetoFamily"}, to enhance accuracy,
+ the method for \code{"GEVFamilyMuUnknown"} uses
+ integration via the quantile transform.}
+\item{makeIC}{\code{signature(IC="IC", L2Fam = "ParetoFamily")}:
+ As with \code{"checkIC"}, to enhance accuracy,
+ the method for \code{"makeIC"} for \code{"ParetoFamily"} uses
+ integration via the quantile transform.}
+\item{makeIC}{\code{signature(IC="IC", L2Fam = "GParetoFamily")}:
+ As for \code{"ParetoFamily"}, to enhance accuracy,
+ the method for \code{"GParetoFamily"} uses
+ integration via the quantile transform.}
+\item{makeIC}{\code{signature(IC="IC", L2Fam = "GEVFamily")}:
+ As for \code{"ParetoFamily"}, to enhance accuracy,
+ the method for \code{"GEVFamily"} uses
+ integration via the quantile transform.}
+\item{makeIC}{\code{signature(IC="IC", L2Fam = "GEVFamilyMuUnknown")}:
+ As for \code{"ParetoFamily"}, to enhance accuracy,
+ the method for \code{"GEVFamilyMuUnknown"} uses
+ integration via the quantile transform.}
+}}
+
+\author{Peter Ruckdeschel \email{peter.ruckdeschel at uni-oldenburg.de}}
+\seealso{\code{\link[RobAStBase]{checkIC}},\code{\link[RobAStBase:makeIC-methods]{makeIC}}}
+\concept{influence curve}
+\keyword{robust}
Modified: branches/robast-1.1/pkg/RobExtremes/man/getStartIC-methods.Rd
===================================================================
--- branches/robast-1.1/pkg/RobExtremes/man/getStartIC-methods.Rd 2018-07-22 11:19:26 UTC (rev 1015)
+++ branches/robast-1.1/pkg/RobExtremes/man/getStartIC-methods.Rd 2018-07-22 11:33:47 UTC (rev 1016)
@@ -29,7 +29,9 @@
\item{getStartIC}{\code{signature(model = "ParetoFamily", risk = "interpolRisk")}:
computes the optimally robust influence function by interpolation
on a grid (using internal helper function \code{.getPsi.P}).}
-}}
+}
+All of these methods recenter and restandardize the obtained ICs to
+warrant centeredness and Fisher consistency.}
\value{
An IC of type \code{HampIC}.
}
Modified: branches/robast-1.1/pkg/RobExtremes/man/ismevgpdgevdiag-methods.Rd
===================================================================
--- branches/robast-1.1/pkg/RobExtremes/man/ismevgpdgevdiag-methods.Rd 2018-07-22 11:19:26 UTC (rev 1015)
+++ branches/robast-1.1/pkg/RobExtremes/man/ismevgpdgevdiag-methods.Rd 2018-07-22 11:33:47 UTC (rev 1016)
@@ -120,14 +120,11 @@
ppfit <- ismev::gev.fit(portpirie[,2])
gev.diag(ppfit)
##
- ## not tested on CRAN because it takes some time...
- \donttest{
- mlE <- MLEstimator(portpirie[,2], GEVFamilyMuUnknown(withPos=FALSE))
- gev.diag(mlE)
+ mlE <- MLEstimator(portpirie[,2], GEVFamilyMuUnknown(withPos=FALSE))
+ gev.diag(mlE)
- gev.prof(mlE, m = 10, 4.1, 5)
- gev.profxi(mlE, -0.3, 0.3)
- }
+ gev.prof(mlE, m = 10, 4.1, 5)
+ gev.profxi(mlE, -0.3, 0.3)
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/robast -r 1016
More information about the Robast-commits
mailing list