[Robast-commits] r904 - in pkg/RobExtremes: . R inst inst/AddMaterial/interpolation inst/scripts inst/unitTests man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sun Sep 4 19:09:57 CEST 2016
Author: ruckdeschel
Date: 2016-09-04 19:09:57 +0200 (Sun, 04 Sep 2016)
New Revision: 904
Added:
pkg/RobExtremes/R/checkEstClassForParamFamiliyMethods.R
pkg/RobExtremes/R/getCVaR.R
pkg/RobExtremes/R/gevgpddiag.R
pkg/RobExtremes/R/startEstGEV.R
pkg/RobExtremes/R/startEstGPD.R
pkg/RobExtremes/inst/scripts/GEVcheck.R
pkg/RobExtremes/man/getCVaR.Rd
pkg/RobExtremes/man/internal-methods.Rd
pkg/RobExtremes/man/ismevgpdgevdiag-methods.Rd
Modified:
pkg/RobExtremes/DESCRIPTION
pkg/RobExtremes/NAMESPACE
pkg/RobExtremes/R/00fromRobAStRDA.R
pkg/RobExtremes/R/AllClass.R
pkg/RobExtremes/R/AllGeneric.R
pkg/RobExtremes/R/AllInitialize.R
pkg/RobExtremes/R/Functionals.R
pkg/RobExtremes/R/GEV.R
pkg/RobExtremes/R/GEVFamily.R
pkg/RobExtremes/R/GEVFamily.R.bak
pkg/RobExtremes/R/GEVFamilyMuUnknown.R
pkg/RobExtremes/R/GParetoFamily.R
pkg/RobExtremes/R/LDEstimator.R
pkg/RobExtremes/R/Pareto.R
pkg/RobExtremes/R/ParetoFamily.R
pkg/RobExtremes/R/WeibullFamily.R
pkg/RobExtremes/R/bdpPickands.R
pkg/RobExtremes/R/getStartIC.R
pkg/RobExtremes/R/internal-getpsi.R
pkg/RobExtremes/R/interpolLM.R
pkg/RobExtremes/R/interpolSn.R
pkg/RobExtremes/R/plotOutlyingness.R
pkg/RobExtremes/inst/AddMaterial/interpolation/SnTest.Rdata
pkg/RobExtremes/inst/AddMaterial/interpolation/getLMInterpol.R
pkg/RobExtremes/inst/AddMaterial/interpolation/interpolationscripts.R
pkg/RobExtremes/inst/AddMaterial/interpolation/plotInterpol.R
pkg/RobExtremes/inst/NEWS
pkg/RobExtremes/inst/TOBEDONE
pkg/RobExtremes/inst/unitTests/runit.expectation.R
pkg/RobExtremes/man/0RobExtremes-package.Rd
pkg/RobExtremes/man/GEVFamily.Rd
pkg/RobExtremes/man/GEVFamilyMuUnknown.Rd
pkg/RobExtremes/man/GEVParameter-class.Rd
pkg/RobExtremes/man/GParetoFamily.Rd
pkg/RobExtremes/man/Var.Rd
pkg/RobExtremes/man/internal-interpolate.Rd
Log:
merged RobExtremes from branches 1.0 back into trunk
Modified: pkg/RobExtremes/DESCRIPTION
===================================================================
--- pkg/RobExtremes/DESCRIPTION 2016-09-04 16:52:44 UTC (rev 903)
+++ pkg/RobExtremes/DESCRIPTION 2016-09-04 17:09:57 UTC (rev 904)
@@ -1,17 +1,26 @@
Package: RobExtremes
-Version: 0.9
-Date: 2013-09-11
-Title: Optimally robust estimation for extreme value distributions
+Version: 1.0
+Date: 2016-09-04
+Title: Optimally Robust Estimation for Extreme Value Distributions
Description: Optimally robust estimation for extreme value distributions using S4 classes and methods (based on
- packages distr, distrEx, distrMod, RobAStBase, and ROptEst)
-Depends: R (>= 2.14.0), evd, RobAStBase(>= 0.9)
-Imports: methods, distr(>= 2.5.2), distrEx(>= 2.4), distrMod(>= 2.5.2), RobAStRDA, actuar, robustbase(>= 0.8-0), RandVar(>= 0.9.2), ROptEst(>= 0.9)
-Suggests: RUnit (>= 0.4.26)
-Author: Peter Ruckdeschel, Matthias Kohl, Nataliya Horbenko
-Maintainer: Peter Ruckdeschel <peter.ruckdeschel at itwm.fraunhofer.de>
+ packages distr, distrEx, distrMod, RobAStBase, and ROptEst).
+Depends: R (>= 2.14.0), methods, distrMod(>= 2.5.2), ROptEst(>= 1.0), robustbase(>= 0.8-0), evd, actuar
+Suggests: RUnit (>= 0.4.26), ismev (>= 1.39)
+Imports: RobAStRDA, distr, distrEx, RandVar, RobAStBase, startupmsg
+Authors at R: c(person("Nataliya", "Horbenko", role=c("aut","cph")),
+ person("Bernhard", "Spangl", role="ctb", comment="contributed smoothed grid values of the Lagrange multipliers"),
+ person("Sascha", "Desmettre", role="ctb", comment="contributed smoothed grid values of the Lagrange multipliers"),
+ person("Eugen", "Massini", role="ctb", comment="contributed an interactive smoothing routine for smoothing the Lagrange multipliers
+ and smoothed grid values of the Lagrange multipliers"),
+ person("Daria", "Pupashenko", role="ctb", comment="contributed MDE-estimation for GEV distribution in the framework of her PhD thesis 2011--14"),
+ person("Gerald", "Kroisandt", role="ctb", comment="contributed testing routines"),
+ person("Matthias", "Kohl", role=c("aut", "cph")),
+ person("Peter", "Ruckdeschel", role=c("cre", "cph"), email="peter.ruckdeschel at uni-oldenburg.de"))
+ByteCompile: yes
LazyLoad: yes
-ByteCompile: yes
License: LGPL-3
+Encoding: latin1
+Additional_repositories: http://stamats.github.io/robast/
URL: http://robast.r-forge.r-project.org/
LastChangedDate: {$LastChangedDate: 2011-09-30 11:10:33 +0200 (Fr, 30 Sep 2011) $}
LastChangedRevision: {$LastChangedRevision: 453 $}
Modified: pkg/RobExtremes/NAMESPACE
===================================================================
--- pkg/RobExtremes/NAMESPACE 2016-09-04 16:52:44 UTC (rev 903)
+++ pkg/RobExtremes/NAMESPACE 2016-09-04 17:09:57 UTC (rev 904)
@@ -11,12 +11,16 @@
import("ROptEst")
import("actuar")
import("evd")
+importFrom("startupmsg", "buildStartupMessage", "infoShow")
+importFrom("stats", "dunif", "integrate", "optimize", "pnorm",
+ "predict", "qnorm", "quantile", "smooth.spline", "uniroot")
+importFrom("utils", "getFromNamespace")
exportClasses("GumbelParameter",
"ParetoParameter",
- "GParetoParameter",
- "GEVParameter",
- "LDEstimate")
+ "GParetoParameter",
+ "GEVParameter",
+ "LDEstimate")
exportClasses("Gumbel", "Pareto", "GPareto", "GEV")
exportClasses("GParetoFamily", "GumbelLocationFamily", "WeibullFamily",
"ParetoFamily", "GEVFamily", "GEVFamilyMuUnknown")
@@ -32,9 +36,10 @@
"E", "var", "IQR", "skewness", "kurtosis", "median", "dispersion")
exportMethods("modifyModel", "getStartIC")
exportMethods("moveL2Fam2RefParam",
- "moveICBackFromRefParam")
+ "moveICBackFromRefParam")
export("EULERMASCHERONICONSTANT","APERYCONSTANT")
+export("getCVaR", "getVaR", "getEL")
export("Gumbel", "Pareto", "GPareto", "GEV")
export("GParetoFamily", "GumbelLocationFamily", "WeibullFamily", "GEVFamily",
"ParetoFamily", "GEVFamilyMuUnknown")
@@ -44,3 +49,9 @@
export("loc", "loc<-", "kMAD", "Sn", "Qn",
"asvarMedkMAD","asvarPickands", "asvarQBCC")
exportMethods("rescaleFunction")
+S3method(print, riskMeasure)
+exportMethods("gev.diag", "gpd.diag","gev.prof", "gpd.prof",
+ "gev.profxi", "gpd.profxi")
+export("gev.diag", "gpd.diag","gev.prof", "gpd.prof",
+ "gev.profxi", "gpd.profxi")
+
\ No newline at end of file
Modified: pkg/RobExtremes/R/00fromRobAStRDA.R
===================================================================
--- pkg/RobExtremes/R/00fromRobAStRDA.R 2016-09-04 16:52:44 UTC (rev 903)
+++ pkg/RobExtremes/R/00fromRobAStRDA.R 2016-09-04 17:09:57 UTC (rev 904)
@@ -1,3 +1,66 @@
-.versionSuff <- RobAStRDA:::.versionSuff
-.MakeSmoothGridList <- RobAStRDA:::.MakeSmoothGridList
+#------------------------------------
+#### utilities copied from package RobAStRDA v.1.0 svn-rev 767
+#------------------------------------
+.versionSuff <- function(name){
+ paste(sep="", name, if(getRversion()<"2.16") ".O" else ".N")
+}
+.MakeSmoothGridList <- function(thGrid, Y, df = NULL,
+ gridRestrForSmooth = NULL){
+ if(length(dim(Y))==3)
+ LMGrid <- Y[,1,,drop=TRUE]
+ else LMGrid <- Y[,drop=FALSE]
+
+ if(!is.null(df)){
+ df0 <- vector("list",ncol(LMGrid))
+ if(is.numeric(df)){
+ df <- rep(df,length.out=ncol(LMGrid))
+ for(i in 1:ncol(LMGrid)) df0[[i]] <- df[i]
+ df <- df0
+ }
+ }else{
+ df0 <- vector("list",ncol(LMGrid)+1)
+ df0[[ncol(LMGrid)+1]] <- NULL
+ df <- df0
+ }
+
+ iNA <- apply(LMGrid,1, function(u) any(is.na(u)))
+ LMGrid <- LMGrid[!iNA,,drop=FALSE]
+ thGrid <- thGrid[!iNA]
+ oG <- order(thGrid)
+ thGrid <- thGrid[oG]
+ LMGrid <- LMGrid[oG,,drop=FALSE]
+
+ if(is.null(gridRestrForSmooth))
+ gridRestrForSmooth <- as.data.frame(matrix(TRUE,nrow(LMGrid),ncol(LMGrid)))
+ if((is.vector(gridRestrForSmooth)&&!is.list(gridRestrForSmooth))||
+ is.matrix(gridRestrForSmooth))
+ gridRestrForSmooth <- as.data.frame(gridRestrForSmooth)
+
+ if(is.list(gridRestrForSmooth)){
+ gm <- vector("list",ncol(LMGrid))
+ idx <- rep(1:length(gridRestrForSmooth), length.out=ncol(LMGrid))
+ for (i in 1:ncol(LMGrid)){
+ if(!is.null(gridRestrForSmooth[[idx[i]]])){
+ gm[[i]] <- gridRestrForSmooth[[idx[i]]]
+ }else{
+ gm[[i]] <- rep(TRUE,nrow(LMGrid))
+ }
+ }
+ gridRestrForSmooth <- gm
+ }
+
+ for(i in 1:ncol(LMGrid)){
+ gmi <- gridRestrForSmooth[[i]]
+ if(is.null(df[[i]])){
+ SmoothSpline <- smooth.spline(thGrid[gmi], LMGrid[gmi, i])
+ LMGrid[, i] <- predict(SmoothSpline, thGrid)$y
+ } else {
+ SmoothSpline <- smooth.spline(thGrid[gmi], LMGrid[gmi, i],
+ df = df[[i]])
+ LMGrid[, i] <- predict(SmoothSpline, thGrid)$y
+ }
+ }
+ return(cbind(xi=thGrid,LM=LMGrid))
+}
+
Modified: pkg/RobExtremes/R/AllClass.R
===================================================================
--- pkg/RobExtremes/R/AllClass.R 2016-09-04 16:52:44 UTC (rev 903)
+++ pkg/RobExtremes/R/AllClass.R 2016-09-04 17:09:57 UTC (rev 904)
@@ -277,3 +277,14 @@
mat = matrix(1))
),
contains = "Estimate")
+
+setOldClass("gev.fit")
+setOldClass("gpd.fit")
+setClass("GPDEstimate", contains="Estimate")
+setClass("GPDMCEstimate", contains=c("MCEstimate", "GPDEstimate"))
+setClass("GPDLDEstimate", contains=c("LDEstimate", "GPDEstimate"))
+setClass("GPDkStepEstimate", contains=c("kStepEstimate", "GPDEstimate"))
+setClass("GEVEstimate", contains="Estimate")
+setClass("GEVLDEstimate", contains=c("LDEstimate", "GEVEstimate"))
+setClass("GEVkStepEstimate", contains=c("kStepEstimate", "GEVEstimate"))
+setClass("GEVMCEstimate", contains=c("MCEstimate", "GEVEstimate"))
Modified: pkg/RobExtremes/R/AllGeneric.R
===================================================================
--- pkg/RobExtremes/R/AllGeneric.R 2016-09-04 16:52:44 UTC (rev 903)
+++ pkg/RobExtremes/R/AllGeneric.R 2016-09-04 17:09:57 UTC (rev 904)
@@ -28,3 +28,21 @@
if(!isGeneric(".loc")){
setGeneric(".loc", function(L2Fam, ...) standardGeneric(".loc"))
}
+if(!isGeneric("gev.diag")){
+ setGeneric("gev.diag", function(z) standardGeneric("gev.diag"))
+}
+if(!isGeneric("gev.prof")){
+ setGeneric("gev.prof", function(z, ...) standardGeneric("gev.prof"))
+}
+if(!isGeneric("gev.profxi")){
+ setGeneric("gev.profxi", function(z, ...) standardGeneric("gev.profxi"))
+}
+if(!isGeneric("gpd.diag")){
+ setGeneric("gpd.diag", function(z,...) standardGeneric("gpd.diag"))
+}
+if(!isGeneric("gpd.prof")){
+ setGeneric("gpd.prof", function(z, ...) standardGeneric("gpd.prof"))
+}
+if(!isGeneric("gpd.profxi")){
+ setGeneric("gpd.profxi", function(z, ...) standardGeneric("gpd.profxi"))
+}
Modified: pkg/RobExtremes/R/AllInitialize.R
===================================================================
--- pkg/RobExtremes/R/AllInitialize.R 2016-09-04 16:52:44 UTC (rev 903)
+++ pkg/RobExtremes/R/AllInitialize.R 2016-09-04 17:09:57 UTC (rev 904)
@@ -14,7 +14,7 @@
body(.Object at p) <- substitute({p1 <- pgumbel(q, loc = loc1, scale = scale1, lower.tail = lower.tail)
return(if(log.p) log(p1) else p1)},
list(loc1 = loc, scale1 = scale))
- .Object at q <- function(p, loc = loc1, scale = scale1, lower.tail = TRUE, log.p = FALSE){}
+ .Object at q <- function(p, lower.tail = TRUE, log.p = FALSE){}
body(.Object at q) <- substitute({
## P.R.: changed to vectorized form
p1 <- if(log.p) exp(p) else p
Modified: pkg/RobExtremes/R/Functionals.R
===================================================================
--- pkg/RobExtremes/R/Functionals.R 2016-09-04 16:52:44 UTC (rev 903)
+++ pkg/RobExtremes/R/Functionals.R 2016-09-04 17:09:57 UTC (rev 904)
@@ -55,7 +55,7 @@
return(var(as(x,"AbscontDistribution"),...))
else{ xi <- shape(x); sigma <- scale(x)
if(xi>=1/2) return(NA)
- if(xi==0) return(pi^2/6)
+ if(xi==0) return(sigma^2*pi^2/6)
if((xi!=0)&&(xi<1/2))return(sigma^2*(gamma(1-2*xi)-gamma(1-xi)^2)/xi^2)
}})
### http://en.wikipedia.org/wiki/Generalized_extreme_value_distribution
Modified: pkg/RobExtremes/R/GEV.R
===================================================================
--- pkg/RobExtremes/R/GEV.R 2016-09-04 16:52:44 UTC (rev 903)
+++ pkg/RobExtremes/R/GEV.R 2016-09-04 17:09:57 UTC (rev 904)
@@ -1,88 +1,88 @@
-###########################################
-## Class: Generalized Extreme Value Distribution
-##
-## @param: location, scale, shape
-###########################################
-
-## access methods
-setMethod("location", "GEVParameter", function(object) object at loc)
-setMethod("loc", "GEVParameter", function(object) object at loc)
-setMethod("scale", "GEVParameter",
- function(x, center = TRUE, scale = TRUE) x at scale)
- ### odd arg-list due to existing function in base package
-setMethod("shape", "GEVParameter", function(object) object at shape)
-
-## replace Methods
-setReplaceMethod("loc", "GEVParameter",
- function(object, value){ object at loc <- value; object })
-setReplaceMethod("location", "GEVParameter",
- function(object, value){ object at loc <- value; object })
-setReplaceMethod("scale", "GEVParameter",
- function(object, value){
- if(length(value) != 1 || value <= 0)
- stop("'value' has to be a single positive real number!")
- object at scale <- value; object})
-setReplaceMethod("shape", "GEVParameter",
- function(object, value){ object at shape <- value; object})
-
-## wrapped access methods
-setMethod("location", "GEV", function(object) loc(object at param))
-setMethod("loc", "GEV", function(object) loc(object at param))
-setMethod("scale", "GEV",
- function(x, center = TRUE, scale = TRUE)
- scale(x at param))
-setMethod("shape", "GEV", function(object) shape(object at param))
-
-
-## wrapped replace methods
-setMethod("loc<-", "GEV", function(object, value)
- new("GEV", loc = value, scale = scale(object), shape = shape(object)))
-setMethod("location<-", "GEV", function(object, value)
- new("GEV", loc = value, scale = scale(object), shape = shape(object)))
-setMethod("scale<-", "GEV", function(object, value)
- new("GEV", loc = loc(object), scale = value, shape = shape(object)))
-setMethod("shape<-", "GEV", function(object, value)
- new("GEV", loc = loc(object), scale = scale(object), shape = value))
-
-setValidity("GEVParameter", function(object){
- if(length(loc(object)) != 1)
- stop("location has to be a numeric of length 1")
- if(length(scale(object)) != 1)
- stop("scale has to be a numeric of length 1")
- if(scale(object) <= 0)
- stop("scale has to be positive")
- if(length(shape(object)) != 1)
- stop("shape has to be a numeric of length 1")
-# if(shape(object) < 0)
-# stop("shape has to be non-negative")
- else return(TRUE)
-})
-
-## generating function
-GEV <- function(loc = 0, scale = 1, shape = 0, location = loc){
- if(!missing(loc)&&!missing(location))
- if(!isTRUE(all.equal(loc,location)))
- stop("Only one of arguments 'loc' and 'location' may be used.")
- if(!missing(location)) loc <- location
- #if(abs(shape) < .Machine$double.eps) return(Gumbel(loc=loc,scale=scale))
- new("GEV", loc = loc, scale = scale, shape = shape) }
-
-## extra methods for GEV distribution
-setMethod("+", c("GEV","numeric"),
- function(e1, e2){
- if (length(e2)>1) stop("length of operator must be 1")
- new("GEV", loc = loc(e1) + e2, scale = scale(e1), shape=shape(e1))
- })
-
-setMethod("*", c("GEV","numeric"),
- function(e1, e2){
- if (length(e2)>1) stop("length of operator must be 1")
- if (isTRUE(all.equal(e2,0)))
- return(new("Dirac", location = 0, .withArith = TRUE))
- GEV <- new("GEV", loc = loc(e1) * abs(e2),
- scale = scale(e1)*abs(e2), shape=shape(e1))
- if(e2<0) GEV <- (-1)*as(GEV,"AbscontDistribution")
- return(GEV)
- })
-
-
+###########################################
+## Class: Generalized Extreme Value Distribution
+##
+## @param: location, scale, shape
+###########################################
+
+## access methods
+setMethod("location", "GEVParameter", function(object) object at loc)
+setMethod("loc", "GEVParameter", function(object) object at loc)
+setMethod("scale", "GEVParameter",
+ function(x, center = TRUE, scale = TRUE) x at scale)
+ ### odd arg-list due to existing function in base package
+setMethod("shape", "GEVParameter", function(object) object at shape)
+
+## replace Methods
+setReplaceMethod("loc", "GEVParameter",
+ function(object, value){ object at loc <- value; object })
+setReplaceMethod("location", "GEVParameter",
+ function(object, value){ object at loc <- value; object })
+setReplaceMethod("scale", "GEVParameter",
+ function(object, value){
+ if(length(value) != 1 || value <= 0)
+ stop("'value' has to be a single positive real number!")
+ object at scale <- value; object})
+setReplaceMethod("shape", "GEVParameter",
+ function(object, value){ object at shape <- value; object})
+
+## wrapped access methods
+setMethod("location", "GEV", function(object) loc(object at param))
+setMethod("loc", "GEV", function(object) loc(object at param))
+setMethod("scale", "GEV",
+ function(x, center = TRUE, scale = TRUE)
+ scale(x at param))
+setMethod("shape", "GEV", function(object) shape(object at param))
+
+
+## wrapped replace methods
+setMethod("loc<-", "GEV", function(object, value)
+ new("GEV", loc = value, scale = scale(object), shape = shape(object)))
+setMethod("location<-", "GEV", function(object, value)
+ new("GEV", loc = value, scale = scale(object), shape = shape(object)))
+setMethod("scale<-", "GEV", function(object, value)
+ new("GEV", loc = loc(object), scale = value, shape = shape(object)))
+setMethod("shape<-", "GEV", function(object, value)
+ new("GEV", loc = loc(object), scale = scale(object), shape = value))
+
+setValidity("GEVParameter", function(object){
+ if(length(loc(object)) != 1)
+ stop("location has to be a numeric of length 1")
+ if(length(scale(object)) != 1)
+ stop("scale has to be a numeric of length 1")
+ if(scale(object) <= 0)
+ stop("scale has to be positive")
+ if(length(shape(object)) != 1)
+ stop("shape has to be a numeric of length 1")
+# if(shape(object) < 0)
+# stop("shape has to be non-negative")
+ else return(TRUE)
+})
+
+## generating function
+GEV <- function(loc = 0, scale = 1, shape = 0, location = loc){
+ if(!missing(loc)&&!missing(location))
+ if(!isTRUE(all.equal(loc,location)))
+ stop("Only one of arguments 'loc' and 'location' may be used.")
+ if(!missing(location)) loc <- location
+ #if(abs(shape) < .Machine$double.eps) return(Gumbel(loc=loc,scale=scale))
+ new("GEV", loc = loc, scale = scale, shape = shape) }
+
+## extra methods for GEV distribution
+setMethod("+", c("GEV","numeric"),
+ function(e1, e2){
+ if (length(e2)>1) stop("length of operator must be 1")
+ new("GEV", loc = loc(e1) + e2, scale = scale(e1), shape=shape(e1))
+ })
+
+setMethod("*", c("GEV","numeric"),
+ function(e1, e2){
+ if (length(e2)>1) stop("length of operator must be 1")
+ if (isTRUE(all.equal(e2,0)))
+ return(new("Dirac", location = 0, .withArith = TRUE))
+ GEV <- new("GEV", loc = loc(e1) * abs(e2),
+ scale = scale(e1)*abs(e2), shape=shape(e1))
+ if(e2<0) GEV <- (-1)*as(GEV,"AbscontDistribution")
+ return(GEV)
+ })
+
+
Modified: pkg/RobExtremes/R/GEVFamily.R
===================================================================
--- pkg/RobExtremes/R/GEVFamily.R 2016-09-04 16:52:44 UTC (rev 903)
+++ pkg/RobExtremes/R/GEVFamily.R 2016-09-04 17:09:57 UTC (rev 904)
@@ -15,34 +15,46 @@
}
### pretreatment of of.interest
-.pretreat.of.interest <- function(of.interest,trafo){
+.pretreat.of.interest <- function(of.interest,trafo,withMu=FALSE){
if(is.null(trafo)){
of.interest <- unique(of.interest)
- if(length(of.interest) > 2)
+ if(!withMu && length(of.interest) > 2)
stop("A maximum number of two parameters resp. parameter transformations may be selected.")
- if(!all(of.interest %in% c("scale", "shape", "quantile", "expected loss", "expected shortfall")))
+ if(withMu && length(of.interest) > 3)
+ stop("A maximum number of three parameters resp. parameter transformations may be selected.")
+ if(!withMu && !all(of.interest %in% c("scale", "shape", "quantile", "expected loss", "expected shortfall")))
stop("Parameters resp. transformations of interest have to be selected from: ",
"'scale', 'shape', 'quantile', 'expected loss', 'expected shortfall'.")
+ if(withMu && !all(of.interest %in% c("loc", "scale", "shape", "quantile", "expected loss", "expected shortfall")))
+ stop("Parameters resp. transformations of interest have to be selected from: ",
+ "'loc', 'scale', 'shape', 'quantile', 'expected loss', 'expected shortfall'.")
## reordering of of.interest
- if(("scale" %in% of.interest) && ("scale" != of.interest[1])){
- of.interest[2] <- of.interest[1]
- of.interest[1] <- "scale"
+ muAdd <- 0
+ if(withMu & "loc" %in% of.interest){
+ muAdd <- 1
+ muWhich <- which(of.interest=="loc")
+ notmuWhich <- which(!of.interest %in% "loc")
+ of.interest <- of.interest[c(muWhich,notmuWhich)]
}
- if(!("scale" %in% of.interest) && ("shape" %in% of.interest) && ("shape" != of.interest[1])){
- of.interest[2] <- of.interest[1]
- of.interest[1] <- "shape"
+ if(("scale" %in% of.interest) && ("scale" != of.interest[1+muAdd])){
+ of.interest[2+muAdd] <- of.interest[1+muAdd]
+ of.interest[1+muAdd] <- "scale"
}
+ if(!("scale" %in% of.interest) && ("shape" %in% of.interest) && ("shape" != of.interest[1+muAdd])){
+ of.interest[2+muAdd] <- of.interest[1+muAdd]
+ of.interest[1+muAdd] <- "shape"
+ }
if(!any(c("scale", "shape") %in% of.interest) && ("quantile" %in% of.interest)
- && ("quantile" != of.interest[1])){
- of.interest[2] <- of.interest[1]
- of.interest[1] <- "quantile"
+ && ("quantile" != of.interest[1+muAdd])){
+ of.interest[2+muAdd] <- of.interest[1+muAdd]
+ of.interest[1+muAdd] <- "quantile"
}
if(!any(c("scale", "shape", "quantile") %in% of.interest)
&& ("expected shortfall" %in% of.interest)
- && ("expected shortfall" != of.interest[1])){
- of.interest[2] <- of.interest[1]
- of.interest[1] <- "expected shortfall"
+ && ("expected shortfall" != of.interest[1+muAdd])){
+ of.interest[2+muAdd] <- of.interest[1+muAdd]
+ of.interest[1+muAdd] <- "expected shortfall"
}
}
return(of.interest)
@@ -74,12 +86,20 @@
}else{
tau1 <- tau
tau <- function(theta){ }
- body(tau) <- substitute({ btq0; c(tau0(theta), q) },
- list(btq0=btq, tau0 = tau1))
+ body(tau) <- substitute({ btq0
+ th0 <- tau0(theta)
+ th <- c(th0, q)
+ names(th) <- c(names(th0),"quantile")
+ th
+ }, list(btq0=btq, tau0 = tau1))
Dtau1 <- Dtau
Dtau <- function(theta){}
- body(Dtau) <- substitute({ bDq0; rbind(Dtau0(theta), D) },
- list(Dtau0 = Dtau1, bDq0 = bDq))
+ body(Dtau) <- substitute({ bDq0
+ D0 <- Dtau0(theta)
+ D1 <- rbind(D0, D)
+ rownames(D1) <- c(rownames(D0),"quantile")
+ D1
+ }, list(Dtau0 = Dtau1, bDq0 = bDq))
}
}
if("expected shortfall" %in% of.interest){
@@ -90,12 +110,18 @@
}else{
tau1 <- tau
tau <- function(theta){ }
- body(tau) <- substitute({ btes0; c(tau0(theta), es) },
- list(tau0 = tau1, btes0=btes))
+ body(tau) <- substitute({ btes0
+ th0 <- tau0(theta)
+ th <- c(th0, es)
+ names(th) <- c(names(th0),"expected shortfall")
+ th}, list(tau0 = tau1, btes0=btes))
Dtau1 <- Dtau
Dtau <- function(theta){}
- body(Dtau) <- substitute({ bDes0; rbind(Dtau0(theta), D) },
- list(Dtau0 = Dtau1, bDes0=bDes))
+ body(Dtau) <- substitute({ bDes0
+ D0 <- Dtau0(theta)
+ D1 <- rbind(D0, D)
+ rownames(D1) <- c(rownames(D0),"expected shortfall")
+ D1}, list(Dtau0 = Dtau1, bDes0=bDes))
}
}
if("expected loss" %in% of.interest){
@@ -106,12 +132,18 @@
}else{
tau1 <- tau
tau <- function(theta){ }
- body(tau) <- substitute({ btel0; c(tau0(theta), el) },
- list(tau0 = tau1, btel0=btel))
+ body(tau) <- substitute({ btel0
+ th0 <- tau0(theta)
+ th <- c(th0, el)
+ names(th) <- c(names(th0),"expected los")
+ th}, list(tau0 = tau1, btel0=btel))
Dtau1 <- Dtau
Dtau <- function(theta){}
- body(Dtau) <- substitute({ bDel0; rbind(Dtau0(theta), D) },
- list(Dtau0 = Dtau1, bDel0=bDel))
+ body(Dtau) <- substitute({ bDel0
+ D0 <- Dtau0(theta)
+ D1 <- rbind(D0, D)
+ rownames(D1) <- c(rownames(D0),"expected loss")
+ D1}, list(Dtau0 = Dtau1, bDel0=bDel))
}
}
trafo <- function(x){ list(fval = tau(x), mat = Dtau(x)) }
@@ -127,8 +159,10 @@
return(FALSE)
if (any(param[1] <= tol))
return(FALSE)
- if (any(param[2] <= tol))
+ if(object at param@withPosRestr) if (any(param[2] <= tol))
return(FALSE)
+ if (any(param[2] <= -1/2))
+ return(FALSE)
return(TRUE)
})
@@ -144,11 +178,13 @@
of.interest = c("scale", "shape"),
p = NULL, N = NULL, trafo = NULL,
start0Est = NULL, withPos = TRUE,
+ secLevel = 0.7,
withCentL2 = FALSE,
withL2derivDistr = FALSE,
- ..ignoreTrafo = FALSE){
+ ..ignoreTrafo = FALSE,
+ ..withWarningGEV = TRUE){
theta <- c(loc, scale, shape)
- .warningGEVShapeLarge(shape)
+ if(..withWarningGEV).warningGEVShapeLarge(shape)
of.interest <- .pretreat.of.interest(of.interest,trafo)
@@ -234,13 +270,18 @@
## starting parameters
startPar <- function(x,...){
mu <- theta[1]
-
+ n <- length(x)
+ epsn <- min(floor(secLevel*sqrt(n))+1,n)
+
## Pickand estimator
if(is.null(start0Est)){
+ ### replaced 20140402: CvMMDE-with xi on Grid
#source("kMedMad_Qn_Estimators.R")
- PF <- GEVFamily(loc = theta[1], scale = theta[2], shape = theta[3])
- e1 <- PickandsEstimator(x,ParamFamily=PF)
- e0 <- estimate(e1)
+ ### replaced 20140402:
+ # 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)
}else{
if(is(start0Est,"function")){
e1 <- start0Est(x, ...)
@@ -250,9 +291,12 @@
e0 <- e0[c("scale", "shape")]
}
# print(e0); print(str(x)); print(head(summary(x))); print(mu)
- if(any(x < mu-e0["scale"]/e0["shape"]))
- stop("some data smaller than 'loc-scale/shape' ")
-
+ if(quantile(e0[2]/e0[1]*(x-mu), epsn/n)< (-1)){
+ if(e0[2]>0)
+ stop("shape is positive and some data smaller than 'loc-scale/shape' ")
+ else
+ stop("shape is negative and some data larger than 'loc-scale/shape' ")
+ }
names(e0) <- NULL
return(e0)
}
@@ -263,9 +307,11 @@
theta <- abs(theta)
}else{
if(!is.null(names(theta))){
+ if(theta["shape"]< (-1/2)) theta["shape"] <- -1/2+1e-4
theta["scale"] <- abs(theta["scale"])
}else{
theta[1] <- abs(theta[1])
+ if(theta[2]< (-1/2)) theta[2] <- -1/2+1e-4
}
}
return(theta)
@@ -273,12 +319,13 @@
modifyPar <- function(theta){
theta <- makeOKPar(theta)
- .warningGEVShapeLarge(theta["shape"])
+ sh <- if(!is.null(names(theta))) theta["shape"] else theta[2]
+ if(..withWarningGEV).warningGEVShapeLarge(sh)
if(!is.null(names(theta))){
sc <- theta["scale"]
sh <- theta["shape"]
}else{
- theta <- abs(theta)
+ # changed 20140402: theta <- abs(theta)
sc <- theta[1]
sh <- theta[2]
}
@@ -291,11 +338,11 @@
sc <- force(main(param)[1])
k <- force(main(param)[2])
tr <- fixed(param)[1]
- .warningGEVShapeLarge(k)
+ if(..withWarningGEV).warningGEVShapeLarge(k)
Lambda1 <- function(x) {
y <- x*0
- ind <- (x > tr-sc/k) # = [later] (x1>0)
+ ind <- if(k>0)(x > tr-sc/k) else (x<tr-sc/k)# = [later] (x1>0)
x <- (x[ind]-tr)/sc
x1 <- 1 + k * x
y[ind] <- (x*(1-x1^(-1/k))-1)/x1/sc
@@ -304,7 +351,7 @@
}
Lambda2 <- function(x) {
y <- x*0
- ind <- (x > tr-sc/k) # = [later] (x1>0)
+ ind <- if(k>0)(x > tr-sc/k) else (x<tr-sc/k)# = [later] (x1>0)
x <- (x[ind]-tr)/sc
x1 <- 1 + k * x
x2 <- x / x1
@@ -327,7 +374,7 @@
FisherInfo.fct <- function(param) {
sc <- force(main(param)[1])
k <- force(main(param)[2])
- .warningGEVShapeLarge(k)
+ if(..withWarningGEV).warningGEVShapeLarge(k)
G20 <- gamma(2*k)
G10 <- gamma(k)
G11 <- digamma(k)*gamma(k)
Modified: pkg/RobExtremes/R/GEVFamily.R.bak
===================================================================
--- pkg/RobExtremes/R/GEVFamily.R.bak 2016-09-04 16:52:44 UTC (rev 903)
+++ pkg/RobExtremes/R/GEVFamily.R.bak 2016-09-04 17:09:57 UTC (rev 904)
@@ -1,412 +1,412 @@
-#################################
-##
-## Class: GEVFamily for positive shape
-##
-################################
-
-### some reusable blocks of code (to avoid redundancy):
-
-.warningGEVShapeLarge <- function(xi){
- if(xi>=4.5)
- warning("A shape estimate larger than 4.5 was produced.\n",
- "Shape parameter values larger than 4.5 are critical\n",
- "in the GEV family as to numerical issues. Be careful with \n",
- "ALE results obtained here; they might be unreliable.")
- }
-
-### pretreatment of of.interest
-.pretreat.of.interest <- function(of.interest,trafo){
- if(is.null(trafo)){
- of.interest <- unique(of.interest)
- if(length(of.interest) > 2)
- stop("A maximum number of two parameters resp. parameter transformations may be selected.")
- if(!all(of.interest %in% c("scale", "shape", "quantile", "expected loss", "expected shortfall")))
- stop("Parameters resp. transformations of interest have to be selected from: ",
- "'scale', 'shape', 'quantile', 'expected loss', 'expected shortfall'.")
-
- ## reordering of of.interest
- if(("scale" %in% of.interest) && ("scale" != of.interest[1])){
- of.interest[2] <- of.interest[1]
- of.interest[1] <- "scale"
- }
- if(!("scale" %in% of.interest) && ("shape" %in% of.interest) && ("shape" != of.interest[1])){
- of.interest[2] <- of.interest[1]
- of.interest[1] <- "shape"
- }
- if(!any(c("scale", "shape") %in% of.interest) && ("quantile" %in% of.interest)
- && ("quantile" != of.interest[1])){
- of.interest[2] <- of.interest[1]
- of.interest[1] <- "quantile"
- }
- if(!any(c("scale", "shape", "quantile") %in% of.interest)
- && ("expected shortfall" %in% of.interest)
- && ("expected shortfall" != of.interest[1])){
- of.interest[2] <- of.interest[1]
- of.interest[1] <- "expected shortfall"
- }
- }
- return(of.interest)
-}
-
-.define.tau.Dtau <- function(of.interest, btq, bDq, btes,
- bDes, btel, bDel, p, N){
- tau <- NULL
- if("scale" %in% of.interest){
- tau <- function(theta){ th <- theta[1]; names(th) <- "scale"; th}
- Dtau <- function(theta){ D <- t(c(1, 0)); rownames(D) <- "scale"; D}
- }
- if("shape" %in% of.interest){
- if(is.null(tau)){
- tau <- function(theta){th <- theta[2]; names(th) <- "shape"; th}
- Dtau <- function(theta){D <- t(c(0,1));rownames(D) <- "shape";D}
- }else{
- tau <- function(theta){th <- theta
- names(th) <- c("scale", "shape"); th}
- Dtau <- function(theta){ D <- diag(2);
- rownames(D) <- c("scale", "shape");D}
- }
- }
- if("quantile" %in% of.interest){
- if(is.null(p)) stop("Probability 'p' has to be specified.")
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/robast -r 904
More information about the Robast-commits
mailing list