[Robast-commits] r906 - in branches/robast-1.1/pkg/RobExtremes: . R inst/AddMaterial/interpolation inst/unitTests

noreply at r-forge.r-project.org noreply at r-forge.r-project.org
Sun Sep 4 19:18:41 CEST 2016


Author: ruckdeschel
Date: 2016-09-04 19:18:41 +0200 (Sun, 04 Sep 2016)
New Revision: 906

Modified:
   branches/robast-1.1/pkg/RobExtremes/DESCRIPTION
   branches/robast-1.1/pkg/RobExtremes/NAMESPACE
   branches/robast-1.1/pkg/RobExtremes/R/AllInitialize.R
   branches/robast-1.1/pkg/RobExtremes/R/GEV.R
   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/Pareto.R
   branches/robast-1.1/pkg/RobExtremes/R/ParetoFamily.R
   branches/robast-1.1/pkg/RobExtremes/R/WeibullFamily.R
   branches/robast-1.1/pkg/RobExtremes/R/bdpPickands.R
   branches/robast-1.1/pkg/RobExtremes/R/internal-getpsi.R
   branches/robast-1.1/pkg/RobExtremes/R/interpolSn.R
   branches/robast-1.1/pkg/RobExtremes/R/plotOutlyingness.R
   branches/robast-1.1/pkg/RobExtremes/inst/AddMaterial/interpolation/SnTest.Rdata
   branches/robast-1.1/pkg/RobExtremes/inst/unitTests/runit.expectation.R
Log:
RobExtremes: merged back branch 1.0 into 1.1

Modified: branches/robast-1.1/pkg/RobExtremes/DESCRIPTION
===================================================================
--- branches/robast-1.1/pkg/RobExtremes/DESCRIPTION	2016-09-04 17:14:02 UTC (rev 905)
+++ branches/robast-1.1/pkg/RobExtremes/DESCRIPTION	2016-09-04 17:18:41 UTC (rev 906)
@@ -8,10 +8,14 @@
 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 itwm.fraunhofer.de"))
+		person("Peter", "Ruckdeschel", role=c("cre", "cph"), email="peter.ruckdeschel at uni-oldenburg.de"))
 ByteCompile: yes
 License: LGPL-3
 Encoding: latin1

Modified: branches/robast-1.1/pkg/RobExtremes/NAMESPACE
===================================================================
--- branches/robast-1.1/pkg/RobExtremes/NAMESPACE	2016-09-04 17:14:02 UTC (rev 905)
+++ branches/robast-1.1/pkg/RobExtremes/NAMESPACE	2016-09-04 17:18:41 UTC (rev 906)
@@ -12,6 +12,9 @@
 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",

Modified: branches/robast-1.1/pkg/RobExtremes/R/AllInitialize.R
===================================================================
--- branches/robast-1.1/pkg/RobExtremes/R/AllInitialize.R	2016-09-04 17:14:02 UTC (rev 905)
+++ branches/robast-1.1/pkg/RobExtremes/R/AllInitialize.R	2016-09-04 17:18:41 UTC (rev 906)
@@ -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: branches/robast-1.1/pkg/RobExtremes/R/GEV.R
===================================================================
--- branches/robast-1.1/pkg/RobExtremes/R/GEV.R	2016-09-04 17:14:02 UTC (rev 905)
+++ branches/robast-1.1/pkg/RobExtremes/R/GEV.R	2016-09-04 17:18:41 UTC (rev 906)
@@ -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: branches/robast-1.1/pkg/RobExtremes/R/GEVFamily.R
===================================================================
--- branches/robast-1.1/pkg/RobExtremes/R/GEVFamily.R	2016-09-04 17:14:02 UTC (rev 905)
+++ branches/robast-1.1/pkg/RobExtremes/R/GEVFamily.R	2016-09-04 17:18:41 UTC (rev 906)
@@ -1,465 +1,465 @@
-#################################
-##
-## 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.")
+        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")
[TRUNCATED]

To get the complete diff run:
    svnlook diff /svnroot/robast -r 906


More information about the Robast-commits mailing list