[Robast-commits] r1033 - in pkg/RobExtremes: . R inst inst/AddMaterial inst/AddMaterial/interpolation inst/scripts man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Mon Jul 23 22:19:56 CEST 2018
Author: ruckdeschel
Date: 2018-07-23 22:19:55 +0200 (Mon, 23 Jul 2018)
New Revision: 1033
Added:
pkg/RobExtremes/R/L2LocScaleShapeUnion-methods.R
pkg/RobExtremes/R/checkIC.R
pkg/RobExtremes/R/getStartICPareto.R
pkg/RobExtremes/R/makeIC.R
pkg/RobExtremes/R/sysdata.rda
pkg/RobExtremes/inst/AddMaterial/LDE4Pareto/
pkg/RobExtremes/inst/AddMaterial/getLMPareto.R
pkg/RobExtremes/inst/scripts/RobFitsAtRealData.R
pkg/RobExtremes/man/checkmakeIC-methods.Rd
pkg/RobExtremes/man/internalEstimatorReturnClasses-class.Rd
pkg/RobExtremes/man/internalProbFamilyClasses-class.Rd
pkg/RobExtremes/man/internalProbFamilyReturnClasses-class.Rd
Removed:
pkg/RobExtremes/man/InternalReturnClasses-class.Rd
Modified:
pkg/RobExtremes/DESCRIPTION
pkg/RobExtremes/NAMESPACE
pkg/RobExtremes/R/AllClass.R
pkg/RobExtremes/R/AllGeneric.R
pkg/RobExtremes/R/Expectation.R
pkg/RobExtremes/R/GEVFamily.R
pkg/RobExtremes/R/GEVFamilyMuUnknown.R
pkg/RobExtremes/R/GParetoFamily.R
pkg/RobExtremes/R/Pareto.R
pkg/RobExtremes/R/ParetoFamily.R
pkg/RobExtremes/R/SnQn.R
pkg/RobExtremes/R/asvarMedkMAD.R
pkg/RobExtremes/R/asvarPickands.R
pkg/RobExtremes/R/checkEstClassForParamFamiliyMethods.R
pkg/RobExtremes/R/getStartIC.R
pkg/RobExtremes/R/gevgpddiag.R
pkg/RobExtremes/R/internal-getpsi.R
pkg/RobExtremes/R/kMAD.R
pkg/RobExtremes/R/plotOutlyingness.R
pkg/RobExtremes/R/rescaleFct.R
pkg/RobExtremes/R/startEstGEV.R
pkg/RobExtremes/R/startEstGPD.R
pkg/RobExtremes/inst/AddMaterial/interpolation/WriteUp-Interpolators.txt
pkg/RobExtremes/inst/NEWS
pkg/RobExtremes/inst/TOBEDONE
pkg/RobExtremes/man/0RobExtremes-package.Rd
pkg/RobExtremes/man/GEVFamily.Rd
pkg/RobExtremes/man/GEVFamilyMuUnknown.Rd
pkg/RobExtremes/man/GParetoFamily.Rd
pkg/RobExtremes/man/LDEstimator.Rd
pkg/RobExtremes/man/Pareto-class.Rd
pkg/RobExtremes/man/ParetoParameter-class.Rd
pkg/RobExtremes/man/Var.Rd
pkg/RobExtremes/man/getCVaR.Rd
pkg/RobExtremes/man/getStartIC-methods.Rd
pkg/RobExtremes/man/internal-interpolate.Rd
pkg/RobExtremes/man/internal-methods.Rd
pkg/RobExtremes/man/ismevgpdgevdiag-methods.Rd
Log:
[RobExtremes] merged branch 1.1 to trunk
Modified: pkg/RobExtremes/DESCRIPTION
===================================================================
--- pkg/RobExtremes/DESCRIPTION 2018-07-23 20:17:24 UTC (rev 1032)
+++ pkg/RobExtremes/DESCRIPTION 2018-07-23 20:19:55 UTC (rev 1033)
@@ -1,21 +1,23 @@
Package: RobExtremes
-Version: 1.0
-Date: 2017-04-25
+Version: 1.1.0
+Date: 2018-07-19
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), methods, distrMod(>= 2.5.2), ROptEst(>= 1.0), robustbase(>= 0.8-0), evd, actuar
-Suggests: 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", "aut", "cph"), email="peter.ruckdeschel at uni-oldenburg.de"))
+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), methods, distrMod(>= 2.5.2), ROptEst(>= 1.0), robustbase, evd
+Suggests: RUnit (>= 0.4.26), ismev (>= 1.39)
+Imports: RobAStRDA, distr, distrEx, RandVar, RobAStBase, startupmsg, actuar
+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", "aut", "cph"),
+ email="peter.ruckdeschel at uni-oldenburg.de"))
ByteCompile: yes
LazyLoad: yes
License: LGPL-3
@@ -23,4 +25,4 @@
URL: http://robast.r-forge.r-project.org/
LastChangedDate: {$LastChangedDate: 2011-09-30 11:10:33 +0200 (Fr, 30 Sep 2011) $}
LastChangedRevision: {$LastChangedRevision: 453 $}
-SVNRevision: 930
+VCS/SVNRevision: 990
Modified: pkg/RobExtremes/NAMESPACE
===================================================================
--- pkg/RobExtremes/NAMESPACE 2018-07-23 20:17:24 UTC (rev 1032)
+++ pkg/RobExtremes/NAMESPACE 2018-07-23 20:19:55 UTC (rev 1033)
@@ -9,8 +9,10 @@
import("robustbase")
import("RobAStBase")
import("ROptEst")
-import("actuar")
-import("evd")
+importFrom("actuar", "qpareto1", "ppareto1", "dpareto1", "rpareto1")
+importFrom("evd", "qgumbel", "pgumbel", "dgumbel", "rgumbel",
+ "qgpd", "pgpd", "dgpd", "rgpd",
+ "qgev", "pgev", "dgev", "rgev")
importFrom("startupmsg", "buildStartupMessage", "infoShow")
importFrom("stats", "dunif", "integrate", "optimize", "pnorm",
"predict", "qnorm", "quantile", "smooth.spline", "uniroot")
@@ -25,6 +27,11 @@
exportClasses("GParetoFamily", "GumbelLocationFamily", "WeibullFamily",
"ParetoFamily", "GEVFamily", "GEVFamilyMuUnknown")
exportClasses("DistributionsIntegratingByQuantiles")
+exportClasses("ParamWithLocAndScaleAndShapeFamParameter")
+exportClasses("L2LocScaleShapeUnion")
+exportClasses("GPDEstimate","GPDMCEstimate","GPDLDEstimate",
+ "GPDkStepEstimate","GEVEstimate","GEVLDEstimate",
+ "GEVkStepEstimate","GEVMCEstimate")
exportMethods("initialize", "show", "rescaleFunction")
exportMethods("loc", "loc<-", "kMAD", "Sn", "Qn")
exportMethods("validParameter",
@@ -34,10 +41,13 @@
"+", "*",
"Min", "Min<-",
"E", "var", "IQR", "skewness", "kurtosis", "median", "dispersion")
+exportMethods(".checkEstClassForParamFamily")
+exportMethods("locscaleshapename","locscalename","scaleshapename",
+ "locationname","scalename","shapename","locscaleshapename<-")
exportMethods("modifyModel", "getStartIC")
exportMethods("moveL2Fam2RefParam",
"moveICBackFromRefParam")
-
+exportMethods("checkIC", "makeIC")
export("EULERMASCHERONICONSTANT","APERYCONSTANT")
export("getCVaR", "getVaR", "getEL")
export("Gumbel", "Pareto", "GPareto", "GEV")
@@ -54,4 +64,4 @@
"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/AllClass.R
===================================================================
--- pkg/RobExtremes/R/AllClass.R 2018-07-23 20:17:24 UTC (rev 1032)
+++ pkg/RobExtremes/R/AllClass.R 2018-07-23 20:19:55 UTC (rev 1033)
@@ -23,10 +23,14 @@
}
+#setClassUnion("ParamWithLocAndScaleAndShapeFamParameterUnion",
+# c("ParamWithScaleFamParameter",
+# "ParamWithShapeFamParameter")
+# )
+
setClass("ParamWithLocAndScaleAndShapeFamParameter",
- contains = c("ParamWithScaleFamParameter",
- "ParamWithShapeFamParameter")
- )
+ contains = c("ParamWithScaleAndShapeFamParameter")
+)
# parameter of Gumbel distribution
@@ -47,10 +51,10 @@
# Gumbel distribution
setClass("Gumbel",
- prototype = prototype(r = function(n){ rgumbel(n, loc = 0, scale = 1) },
- d = function(x, log){ dgumbel(x, loc = 0, scale = 1, log = FALSE) },
+ prototype = prototype(r = function(n){ evd::rgumbel(n, loc = 0, scale = 1) },
+ d = function(x, log){ evd::dgumbel(x, loc = 0, scale = 1, log = FALSE) },
p = function(q, lower.tail = TRUE, log.p = FALSE){
- p0 <- pgumbel(q, loc = 0, scale = 1, lower.tail = lower.tail)
+ p0 <- evd::pgumbel(q, loc = 0, scale = 1, lower.tail = lower.tail)
if(log.p) return(log(p0)) else return(p0)
},
q = function(p, loc = 0, scale = 1, lower.tail = TRUE, log.p = FALSE){
@@ -66,7 +70,7 @@
p0 <- p
p0[ii01] <- if(log.p) log(0.5) else 0.5
- q1 <- qgumbel(p0, loc = 0, scale = 1,
+ q1 <- evd::qgumbel(p0, loc = 0, scale = 1,
lower.tail = lower.tail)
q1[i0] <- if(lower.tail) -Inf else Inf
q1[i1] <- if(!lower.tail) -Inf else Inf
@@ -198,7 +202,7 @@
prototype = prototype(
r = function(n){ rgev(n,loc = 0, scale = 1, shape = 0.5) },
d = function(x, log = FALSE){
- dgev(x, loc = 0, scale = 1, shape = 0.5, log = log)
+ dgev(x, loc = 0, scale = 1, shape = 0.5, log = log)
},
p = function(q, lower.tail = TRUE, log.p = FALSE ){
p0 <- pgev(q, loc = 0, scale = 1, shape = 0.5)
@@ -250,7 +254,7 @@
setClass("WeibullFamily", contains="L2ScaleShapeUnion")
## virtual in-between class for common parts in modifyModel - method
-setClass("L2LocScaleShapeUnion", representation(scaleshapename ="character"),
+setClass("L2LocScaleShapeUnion", representation(locscaleshapename = "character"),
contains = c("L2GroupParamFamily","VIRTUAL")
)
Modified: pkg/RobExtremes/R/AllGeneric.R
===================================================================
--- pkg/RobExtremes/R/AllGeneric.R 2018-07-23 20:17:24 UTC (rev 1032)
+++ pkg/RobExtremes/R/AllGeneric.R 2018-07-23 20:19:55 UTC (rev 1033)
@@ -46,3 +46,16 @@
if(!isGeneric("gpd.profxi")){
setGeneric("gpd.profxi", function(z, ...) standardGeneric("gpd.profxi"))
}
+if(!isGeneric("locscaleshapename")){
+ setGeneric("locscaleshapename", function(object) standardGeneric("locscaleshapename"))
+}
+if(!isGeneric("locscaleshapename<-")){
+ setGeneric("locscaleshapename<-", function(object,value) standardGeneric("locscaleshapename<-"))
+}
+if(!isGeneric("shapename")){
+ setGeneric("shapename", function(object) standardGeneric("shapename"))
+}
+if(!isGeneric("locationname")){
+ setGeneric("locationname", function(object) standardGeneric("locationname"))
+}
+
Modified: pkg/RobExtremes/R/Expectation.R
===================================================================
--- pkg/RobExtremes/R/Expectation.R 2018-07-23 20:17:24 UTC (rev 1032)
+++ pkg/RobExtremes/R/Expectation.R 2018-07-23 20:19:55 UTC (rev 1033)
@@ -57,9 +57,12 @@
dots.withoutUseApply <- dots
useApply <- TRUE
if(!is.null(dots$useApply)) useApply <- dots$useApply
+
dots.withoutUseApply$useApply <- NULL
+ dots.withoutUseApply$stop.on.error <- NULL
+
integrand <- function(x, dfun, ...){ di <- dim(x)
- y <- q(object)(x)##quantile transformation
+ y <- q.l(object)(x)##quantile transformation
if(useApply){
funy <- sapply(y,fun, ...)
dim(y) <- di
@@ -76,12 +79,29 @@
upp <- p(object)(Ib["upp"])
if(is.nan(low)) low <- 0
if(is.nan(upp)) upp <- 1
- return(do.call(distrExIntegrate, c(list(f = integrand,
+
+ if(upp < 0.98){
+ int <- do.call(distrExIntegrate, c(list(f = integrand,
lower = low,
upper = upp,
- rel.tol = rel.tol,
- distr = object, dfun = dunif), dots.withoutUseApply)))
+ rel.tol = rel.tol, stop.on.error = FALSE,
+ distr = object, dfun = dunif), dots.withoutUseApply))
+ }else{
+ int1 <- do.call(distrExIntegrate, c(list(f = integrand,
+ lower = low,
+ upper = 0.98,
+ rel.tol = rel.tol, stop.on.error = FALSE,
+ distr = object, dfun = dunif), dots.withoutUseApply))
+ int2 <- do.call(distrExIntegrate, c(list(f = integrand,
+ lower = 0.98,
+ upper = upp,
+ rel.tol = rel.tol, stop.on.error = FALSE,
+ distr = object, dfun = dunif), dots.withoutUseApply))
+ int <- int1+int2
+ }
+ return(int)
+
})
setMethod("E", signature(object = "GPareto",
Modified: pkg/RobExtremes/R/GEVFamily.R
===================================================================
--- pkg/RobExtremes/R/GEVFamily.R 2018-07-23 20:17:24 UTC (rev 1032)
+++ pkg/RobExtremes/R/GEVFamily.R 2018-07-23 20:19:55 UTC (rev 1033)
@@ -1,465 +1,472 @@
-#################################
-##
-## 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,withMu=FALSE){
- if(is.null(trafo)){
- of.interest <- unique(of.interest)
- if(!withMu && length(of.interest) > 2)
- stop("A maximum number of two parameters resp. parameter transformations may be selected.")
- 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
- 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) && ("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+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+muAdd])){
- of.interest[2+muAdd] <- of.interest[1+muAdd]
- of.interest[1+muAdd] <- "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.")
- if(is.null(tau)){
- tau <- function(theta){ }; body(tau) <- btq
- Dtau <- function(theta){ };body(Dtau) <- bDq
- }else{
- tau1 <- tau
- tau <- function(theta){ }
- 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
- 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){
- if(is.null(p)) stop("Probability 'p' has to be specified.")
- if(is.null(tau)){
- tau <- function(theta){ }; body(tau) <- btes
- Dtau <- function(theta){ }; body(Dtau) <- bDes
- }else{
- tau1 <- tau
- tau <- function(theta){ }
- 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
- 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){
- if(is.null(N)) stop("Expected frequency 'N' has to be specified.")
- if(is.null(tau)){
- tau <- function(theta){ }; body(tau) <- btel
- Dtau <- function(theta){ }; body(Dtau) <- bDel
- }else{
- tau1 <- tau
- tau <- function(theta){ }
- 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
- 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)) }
- return(trafo)
-}
-
-## methods
-setMethod("validParameter",signature(object="GEVFamily"),
- function(object, param, tol =.Machine$double.eps){
- if (is(param, "ParamFamParameter"))
- param <- main(param)
- if (!all(is.finite(param)))
- return(FALSE)
- if (any(param[1] <= tol))
- return(FALSE)
- if(object at param@withPosRestr) if (any(param[2] <= tol))
- return(FALSE)
- if (any(param[2] <= -1/2))
- return(FALSE)
- return(TRUE)
- })
-
-
-## generating function
-## loc: known/fixed threshold/location parameter
-## scale: scale parameter
-## shape: shape parameter
-## trafo: optional parameter transformation
-## start0Est: startEstimator for MLE and MDE --- if NULL HybridEstimator is used;
-
-GEVFamily <- function(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){
- theta <- c(loc, scale, shape)
- if(..withWarningGEV).warningGEVShapeLarge(shape)
-
- of.interest <- .pretreat.of.interest(of.interest,trafo)
-
- ##symmetry
- distrSymm <- NoSymmetry()
-
- ## parameters
- names(theta) <- c("loc", "scale", "shape")
- scaleshapename <- c("scale"="scale", "shape"="shape")
-
-
- btq <- bDq <- btes <- bDes <- btel <- bDel <- NULL
- if(!is.null(p)){
- btq <- substitute({ q <- loc0 + theta[1]*((-log(p0))^(-theta[2])-1)/theta[2]
- names(q) <- "quantile"
- q
- }, list(loc0 = loc, p0 = p))
-
- bDq <- substitute({ scale <- theta[1]; shape <- theta[2]
- D1 <- ((-log(p0))^(-shape)-1)/shape
- D2 <- -scale/shape*(D1 + log(-log(p0))*(-log(p0))^(-shape))
- D <- t(c(D1, D2))
- rownames(D) <- "quantile"; colnames(D) <- NULL
- D }, list(p0 = p))
- btes <- substitute({ if(theta[2]>=1L){
- warning("Expected value is infinite for shape > 1")
- es <- NA
- }else{
- pg <- pgamma(-log(p0),1-theta[2], lower.tail = TRUE)
- es <- theta[1] * (gamma(1-theta[2]) * pg/ (1-p0) - 1 )/
- theta[2] + loc0 }
- names(es) <- "expected shortfall"
- es }, list(loc0 = loc, p0 = p))
- bDes <- substitute({ if(theta[2]>=1L){ D1 <- D2 <- NA} else {
- scale <- theta[1]; shape <- theta[2]
- pg <- pgamma(-log(p0), 1-theta[2], lower.tail = TRUE)
- dd <- ddigamma(-log(p0),1-theta[2])
- g0 <- gamma(1-theta[2])
- D1 <- (g0*pg/(1-p0)-1)/theta[2]
- D21 <- D1/theta[2]
- D22 <- dd/(1-p0)/theta[2]
- D2 <- -theta[1]*(D21+D22)}
- D <- t(c(D1, D2))
- rownames(D) <- "expected shortfall"
- colnames(D) <- NULL
- D }, list(loc0 = loc, p0 = p))
- }
- if(!is.null(N)){
- btel <- substitute({ if(theta[2]>=1L){
- warning("Expected value is infinite for shape > 1")
- el <- NA
- }else{
- el <- N0*(loc0+theta[1]*(gamma(1-theta[2])-1)/theta[2])}
- names(el) <- "expected loss"
- el }, list(loc0 = loc,N0 = N))
- bDel <- substitute({ if(theta[2]>=1L){ D1 <- D2 <- NA}else{
- scale <- theta[1]; shape <- theta[2]
- ga <- gamma(1-shape)
- D1 <- N0*(ga-1)/shape
- D2 <- -N0*scale*ga*digamma(1-shape)/shape-
- D1*scale/shape}
- D <- t(c(D1, D2))
- rownames(D) <- "expected loss"
- colnames(D) <- NULL
- D }, list(loc0 = loc, N0 = N))
- }
-
- fromOfInt <- FALSE
- if(is.null(trafo)||..ignoreTrafo){fromOfInt <- TRUE
- trafo <- .define.tau.Dtau(of.interest, btq, bDq, btes, bDes,
- btel, bDel, p, N)
- }else if(is.matrix(trafo) & nrow(trafo) > 2)
- stop("number of rows of 'trafo' > 2")
-####
- param <- ParamFamParameter(name = "theta", main = c(theta[2],theta[3]),
- fixed = theta[1],
- trafo = trafo, withPosRestr = withPos,
- .returnClsName ="ParamWithScaleAndShapeFamParameter")
-
- ## distribution
- distribution <- GEV(loc = loc, scale = scale, shape = shape)
-
- ## 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")
- ### 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, ...)
- e0 <- if(is(e1,"Estimate")) estimate(e1) else e1
- }else stop("Argument 'start0Est' must be a function or NULL.")
- if(!is.null(names(e0)))
- e0 <- e0[c("scale", "shape")]
- }
-# print(e0); print(str(x)); print(head(summary(x))); print(mu)
- 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)
- }
-
- ## what to do in case of leaving the parameter domain
- makeOKPar <- function(theta) {
- if(withPos){
- 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)
- }
-
- modifyPar <- function(theta){
- theta <- makeOKPar(theta)
- 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{
- # changed 20140402: theta <- abs(theta)
- sc <- theta[1]
- sh <- theta[2]
- }
- GEV(loc = loc, scale = theta[1], shape = theta[2])
- }
-
-
- ## L2-derivative of the distribution
- L2deriv.fct <- function(param) {
- sc <- force(main(param)[1])
- k <- force(main(param)[2])
- tr <- fixed(param)[1]
- if(..withWarningGEV).warningGEVShapeLarge(k)
-
- Lambda1 <- function(x) {
- y <- x*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
-# xi*(-1/xi-1)*(x[ind]-mu)/beta^2/(1+xi*(x[ind]-mu)/beta) - (x[ind]-mu)*(1+xi*(x[ind]-mu)/beta)^(-1/xi-1)/beta^2
- return(y)
- }
- Lambda2 <- function(x) {
- y <- x*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
- y[ind]<- (1-x1^(-1/k))/k*(log(x1)/k-x2)-x2
-# log(1+xi*(x[ind]-mu)/beta)/xi^2+(-1/xi-1)*(x[ind]-mu)/beta/(1+xi*(x[ind]-mu)/beta) - (1+xi*(x[ind]-mu)/beta)^(-1/xi)*log(1+xi*(x[ind]-mu)/beta)/xi^2 + (1+xi*(x[ind]-mu)/beta)^(-1/xi-1)*(x[ind]-mu)/beta/xi
- return(y)
- }
- ## additional centering of scores to increase numerical precision!
- if(withCentL2){
- dist0 <- GEV(scale = sc, shape = k, loc = tr)
- suppressWarnings({
- z1 <- E(dist0, fun=Lambda1)
- z2 <- E(dist0, fun=Lambda2)
- })
- }else{z1 <- z2 <- 0}
- return(list(function(x){ Lambda1(x)-z1 },function(x){ Lambda2(x)-z2 }))
- }
-
- ## Fisher Information matrix as a function of parameters
- FisherInfo.fct <- function(param) {
- sc <- force(main(param)[1])
- k <- force(main(param)[2])
- if(..withWarningGEV).warningGEVShapeLarge(k)
- G20 <- gamma(2*k)
- G10 <- gamma(k)
- G11 <- digamma(k)*gamma(k)
- G01 <- -0.57721566490153 # digamma(1)
- G02 <- 1.9781119906559 #trigamma(1)+digamma(1)^2
- x0 <- (k+1)^2*2*k
- I11 <- G20*x0-2*G10*k*(k+1)+1
- I11 <- I11/sc^2/k^2
- I12 <- G20*(-x0)+ G10*(k^3+4*k^2+3*k) - k -1
- I12 <- I12 + G11*(k^3+k^2) -G01*k
- I12 <- I12/sc/k^3
- I22 <- G20*x0 +(k+1)^2 -G10*(x0+2*k*(k+1))
- I22 <- I22 - G11*2*k^2*(k+1) + G01*2*k*(1+k)+k^2 *G02
- I22 <- I22 /k^4
- mat <- PosSemDefSymmMatrix(matrix(c(I11,I12,I12,I22),2,2))
- dimnames(mat) <- list(scaleshapename,scaleshapename)
- return(mat)
- }
-
-
-
- FisherInfo <- FisherInfo.fct(param)
- name <- "GEV Family"
-
- ## initializing the GPareto family with components of L2-family
- L2Fam <- new("GEVFamily")
- L2Fam at scaleshapename <- scaleshapename
- L2Fam at name <- name
- L2Fam at param <- param
- L2Fam at distribution <- distribution
- L2Fam at L2deriv.fct <- L2deriv.fct
- L2Fam at FisherInfo.fct <- FisherInfo.fct
- L2Fam at FisherInfo <- FisherInfo
- L2Fam at startPar <- startPar
- L2Fam at makeOKPar <- makeOKPar
- L2Fam at modifyParam <- modifyPar
- L2Fam at L2derivSymm <- FunSymmList(NonSymmetric(), NonSymmetric())
- L2Fam at L2derivDistrSymm <- DistrSymmList(NoSymmetry(), NoSymmetry())
-
- L2deriv <- EuclRandVarList(RealRandVariable(L2deriv.fct(param),
- Domain = Reals()))
- L2derivDistr <- NULL
- if(withL2derivDistr){
- suppressWarnings(L2derivDistr <-
- imageDistr(RandVar = L2deriv, distr = distribution))
- }
-
- if(fromOfInt){
- L2Fam at fam.call <- substitute(GEVFamily(loc = loc0, scale = scale0,
- shape = shape0, of.interest = of.interest0,
- p = p0, N = N0,
- withPos = withPos0, withCentL2 = FALSE,
- withL2derivDistr = FALSE, ..ignoreTrafo = TRUE),
- list(loc0 = loc, scale0 = scale, shape0 = shape,
- of.interest0 = of.interest, p0 = p, N0 = N,
- withPos0 = withPos))
- }else{
- L2Fam at fam.call <- substitute(GEVFamily(loc = loc0, scale = scale0,
- shape = shape0, of.interest = NULL,
- p = p0, N = N0, trafo = trafo0,
- withPos = withPos0, withCentL2 = FALSE,
- withL2derivDistr = FALSE),
- list(loc0 = loc, scale0 = scale, shape0 = shape,
- p0 = p, N0 = N,
- withPos0 = withPos, trafo0 = trafo))
- }
-
- L2Fam at LogDeriv <- function(x){
- x0 <- (x-loc)/scale
- x1 <- 1 + x0 * shape
- (shape+1)/scale/x1 + x1^(-1-1/shape)/scale
- }
-
- L2Fam at L2deriv <- L2deriv
- L2Fam at L2derivDistr <- L2derivDistr
- L2Fam at .withMDE <- FALSE
- L2Fam at .withEvalAsVar <- FALSE
- L2Fam at .withEvalL2derivDistr <- FALSE
- return(L2Fam)
-}
-
-#ddigamma(t,s) is d/ds \int_0^t exp(-x) x^(s-1) dx
-
-ddigamma <- function(t,s){
- int <- function(x) exp(-x)*(log(x))*x^(s-1)
- integrate(int, lower=0, upper=t)$value
- }
-
\ No newline at end of file
+#################################
+##
+## 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,withMu=FALSE){
+ if(is.null(trafo)){
+ of.interest <- unique(of.interest)
+ if(!withMu && length(of.interest) > 2)
+ stop("A maximum number of two parameters resp. parameter transformations may be selected.")
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/robast -r 1033
More information about the Robast-commits
mailing list