[Robast-commits] r525 - in branches/robast-0.9/pkg/RobExtremes: R tests tests/TestSuite

noreply at r-forge.r-project.org noreply at r-forge.r-project.org
Mon Nov 5 08:50:52 CET 2012


Author: horbenko
Date: 2012-11-05 08:50:51 +0100 (Mon, 05 Nov 2012)
New Revision: 525

Modified:
   branches/robast-0.9/pkg/RobExtremes/R/GEVFamily.R
   branches/robast-0.9/pkg/RobExtremes/R/drawPointsIC.R
   branches/robast-0.9/pkg/RobExtremes/R/plotOutlyingness.R
   branches/robast-0.9/pkg/RobExtremes/tests/
   branches/robast-0.9/pkg/RobExtremes/tests/TestSuite.R
   branches/robast-0.9/pkg/RobExtremes/tests/TestSuite/TestExpectation.R
   branches/robast-0.9/pkg/RobExtremes/tests/testResult.html
   branches/robast-0.9/pkg/RobExtremes/tests/testResult.txt
Log:
GEV-Family und weitere ?\195?\132nderungen

Modified: branches/robast-0.9/pkg/RobExtremes/R/GEVFamily.R
===================================================================
--- branches/robast-0.9/pkg/RobExtremes/R/GEVFamily.R	2012-10-31 08:23:45 UTC (rev 524)
+++ branches/robast-0.9/pkg/RobExtremes/R/GEVFamily.R	2012-11-05 07:50:51 UTC (rev 525)
@@ -1,22 +1,23 @@
 #################################
 ##
-## Class: GParetoFamily
+## Class: GEVFamily for positive shape
 ##
 ################################
 
+## class
+setClass("GEVFamily", contains="L2ParamFamily")
 
 ## methods
-setMethod("validParameter",signature(object="GParetoFamily"),
+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[1] <= tol)) 
+                 return(FALSE)
+             if (any(param[2] <= tol))
+                 return(FALSE)
              return(TRUE)
            })
 
@@ -25,48 +26,12 @@
 ## loc: known/fixed threshold/location parameter
 ## scale: scale parameter
 ## shape: shape parameter
-## of.interest: which parameters, transformations are of interest
-##              posibilites are: scale, shape, quantile, expected loss, expected shortfall
-##              a maximum number of two of these may be selected
-## p: probability needed for quantile and expected shortfall
-## N: expected frequency for expected loss
 ## trafo: optional parameter transformation
 ## start0Est: startEstimator for MLE and MDE --- if NULL HybridEstimator is used;
-### now uses exp-Trafo for scale!
 
-GParetoFamily <- function(loc = 0, scale = 1, shape = 0.5, 
-                          of.interest = c("scale", "shape"), 
-                          p = NULL, N = NULL, trafo = NULL,
-                          start0Est = NULL, withPos = TRUE){
-    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"
-        }
-    }
+GEVFamily <- function(loc = 0, scale = 1, shape = 0.5,trafo = NULL,start0Est = NULL){
+    if(is.null(trafo)) trafo = diag(2)
+    
     theta <- c(loc, scale, shape)
 
     ##symmetry
@@ -74,196 +39,22 @@
 
     ## parameters
     names(theta) <- c("loc", "scale", "shape")
-    scaleshapename <- c("scale", "shape")
 
-    if(is.null(trafo)){
-        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) <- substitute({ q <- loc0 + theta[1]*((1-p0)^(-theta[2])-1)/theta[2]
-                                          names(q) <- "quantile"
-                                          q },
-                                        list(loc0 = loc, p0 = p))
-                Dtau <- function(theta){ }
-                body(Dtau) <- substitute({ scale <- theta[1]
-                                           shape <- theta[2]
-                                           D1 <- ((1-p0)^(-shape)-1)/shape
-                                           D2 <- -scale/shape*(D1 + log(1-p0)*(1-p0)^(-shape))
-                                           D <- t(c(D1, D2))
-                                           rownames(D) <- "quantile"
-                                           colnames(D) <- NULL
-                                           D },
-                                         list(p0 = p))
-            }else{
-                tau1 <- tau
-                tau <- function(theta){ }
-                body(tau) <- substitute({ q <- loc0 + theta[1]*((1-p0)^(-theta[2])-1)/theta[2]
-                                          names(q) <- "quantile"
-                                          c(tau0(theta), q) },
-                                        list(tau0 = tau1, loc0 = loc, p0 = p))
-                Dtau1 <- Dtau
-                Dtau <- function(theta){}
-                body(Dtau) <- substitute({ scale <- theta[1]
-                                           shape <- theta[2]
-                                           D1 <- ((1-p0)^(-shape)-1)/shape
-                                           D2 <- -scale/shape*(D1 + log(1-p0)*(1-p0)^(-shape))
-                                           D <- t(c(D1, D2))
-                                           rownames(D) <- "quantile"
-                                           colnames(D) <- NULL
-                                           rbind(Dtau0(theta), D) },
-                                         list(Dtau0 = Dtau1, p0 = p))
-            }
-        }
-        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) <- substitute({ q <- loc0 + theta[1]*((1-p0)^(-theta[2])-1)/theta[2]
-                                          es <- (q + theta[1] - theta[2]*loc0)/(1-theta[2]) 
-                                          names(es) <- "expected shortfall"
-                                          es }, 
-                                        list(loc0 = loc, p0 = p))
-                Dtau <- function(theta){ }
-                body(Dtau) <- substitute({ scale <- theta[1]
-                                           shape <- theta[2]
-                                           q <- loc0 + theta[1]*((1-p0)^(-theta[2])-1)/theta[2]
-                                           dq1 <- ((1-p0)^(-shape)-1)/shape
-                                           dq2 <- -scale/shape*(dq1 + log(1-p0)*(1-p0)^(-shape))
-                                           D1 <- (dq1 + 1)/(1-shape)
-                                           D2 <- (dq2 - loc0)/(1-shape) + (q + scale - loc0*shape)/(1-shape)^2
-                                           D <- t(c(D1, D2))
-                                           rownames(D) <- "expected shortfall"
-                                           colnames(D) <- NULL
-                                           D },
-                                         list(loc0 = loc, p0 = p))
-            }else{
-                tau1 <- tau
-                tau <- function(theta){ }
-                body(tau) <- substitute({ q <- loc0 + theta[1]*((1-p0)^(-theta[2])-1)/theta[2]
-                                          es <- (q + theta[1] - theta[2]*loc0)/(1-theta[2]) 
-                                          names(es) <- "expected shortfall"
-                                          c(tau0(theta), es) }, 
-                                        list(tau0 = tau1, loc0 = loc, p0 = p))
-                Dtau1 <- Dtau
-                Dtau <- function(theta){}
-                body(Dtau) <- substitute({ scale <- theta[1]
-                                           shape <- theta[2]
-                                           q <- loc0 + theta[1]*((1-p0)^(-theta[2])-1)/theta[2]
-                                           dq1 <- ((1-p0)^(-shape)-1)/shape
-                                           dq2 <- -scale/shape*(dq1 + log(1-p0)*(1-p0)^(-shape))
-                                           D1 <- (dq1 + 1)/(1-shape)
-                                           D2 <- (dq2 - loc0)/(1-shape) + (q + scale - loc0*shape)/(1-shape)^2
-                                           D <- t(c(D1, D2))
-                                           rownames(D) <- "expected shortfall"
-                                           colnames(D) <- NULL
-                                           rbind(Dtau0(theta), D) },
-                                         list(Dtau0 = Dtau1, loc0 = loc, p0 = p))
-            }
-        }
-        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) <- substitute({ el <- N0*(loc0 + theta[1]*gamma(1/theta[2]-1)/(theta[2]^2*gamma(1/theta[2]+1)))
-                                          names(el) <- "expected loss"
-                                          el },
-                                        list(loc0 = loc,N0 = N))
-                Dtau <- function(theta){ }
-                body(Dtau) <- substitute({ scale <- theta[1]
-                                           shape <- theta[2]
-                                           Gneg <- gamma(1/shape-1)
-                                           Gpos <- gamma(1/shape+1)
-                                           D1 <- N0*Gneg/(shape^2*Gpos)
-                                           D2 <- N0*scale*Gneg*(digamma(1/shape+1) - 2*shape - digamma(1/shape-1))/(shape^4*Gpos)
-                                           D <- t(c(D1, D2))
-                                           rownames(D) <- "expected loss"
-                                           colnames(D) <- NULL
-                                           D },
-                                         list(loc0 = loc, N0 = N))
-            }else{
-                tau1 <- tau
-                tau <- function(theta){ }
-                body(tau) <- substitute({ el <- N0*(loc0 + theta[1]*gamma(1/theta[2]-1)/(theta[2]^2*gamma(1/theta[2]+1)))
-                                          names(el) <- "expected loss"
-                                          c(tau0(theta), el) },
-                                        list(tau0 = tau1, loc0 = loc,N0 = N))
-                Dtau1 <- Dtau
-                Dtau <- function(theta){}
-                body(Dtau) <- substitute({ scale <- theta[1]
-                                           shape <- theta[2]
-                                           Gneg <- gamma(1/shape-1)
-                                           Gpos <- gamma(1/shape+1)
-                                           D1 <- N0*Gneg/(shape^2*Gpos)
-                                           D2 <- N0*scale*Gneg*(digamma(1/shape+1) - 2*shape - digamma(1/shape-1))/(shape^4*Gpos)
-                                           D <- t(c(D1, D2))
-                                           rownames(D) <- "expected loss"
-                                           colnames(D) <- NULL
-                                           rbind(Dtau0(theta), D) },
-                                         list(Dtau0 = Dtau1, loc0 = loc, N0 = N))
-            }
-        }
-        trafo <- function(x){ list(fval = tau(x), mat = Dtau(x)) }
-    }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]),
+    param <- ParamFamParameter(name = "theta", main = theta[2:3], 
                                fixed = theta[1],
-                               trafo = trafo, withPosRestr = withPos,
-                               .returnClsName ="ParamWithScaleAndShapeFamParameter")
+                               trafo = trafo)
 
     ## distribution
-    distribution <- GPareto(loc = loc, scale = scale, shape = shape)
+    distribution <- GEV(loc = loc, scale = scale, shape = shape)
 
     ## starting parameters
     startPar <- function(x,...){
-        tr <- theta[1]
+        mu <- theta[1]
         
-        if(any(x < tr)) stop("some data smaller than 'loc' parameter")
-
         ## Pickand estimator
         if(is.null(start0Est)){
-           e0 <- estimate(medkMADhybr(x, k=10, ParamFamily=GParetoFamily(loc = theta[1],
-                            scale = theta[2], shape = theta[3]),
-                            q.lo = 1e-3, q.up = 15))
+        source("kMedMad_Qn_Estimators.R")
+           e0 <- EVTEst(x,est="kMAD")
         }else{
            if(is(start0Est,"function")){
               e1 <- start0Est(x, ...)
@@ -272,55 +63,41 @@
            if(!is.null(names(e0)))
                e0 <- e0[c("scale", "shape")]
         }
+        if(any(x < mu-e0["scale"]/e0["shape"])) stop("some data smaller than 'loc-scale/shape' ")
+
         names(e0) <- NULL
         return(e0)
     }
 
+    modifyPar <- function(theta){
+        theta <- abs(theta)
+        GEV(loc = loc, scale = theta[1], shape = theta[2])
+    }
 
     ## what to do in case of leaving the parameter domain
     makeOKPar <- function(theta) {
-        if(withPos){
-           if(!is.null(names(theta)))
-                 theta["shape"] <- abs(theta["shape"])
-           else  theta[2] <- abs(theta[2])
-        }
+        theta <- abs(theta)
+        theta[2] <- pmin(theta[2],10)
         return(theta)
     }
 
-    modifyPar <- function(theta){
-        theta <- makeOKPar(theta)
-
-        if(!is.null(names(theta))){
-            sc <- theta["scale"]
-            sh <- theta["shape"]
-        }else{
-            theta <- abs(theta)
-            sc <- theta[1]
-            sh <- theta[2]
-        }
-        GPareto(loc = loc, scale = sc, shape = sh)
-    }
-
-
     ## L2-derivative of the distribution
     L2deriv.fct <- function(param) {
-        sc <- force(main(param)[1])
-        k <- force(main(param)[2])
-        tr <- fixed(param)[1] 
+        beta <- force(main(param)[1])
+        xi <- force(main(param)[2])
+        mu <- fixed(param)[1] 
 
         Lambda1 <- function(x) {
-            y <- x*0
-            x0 <- (x-tr)/sc
-            x1 <- x0[x0>0]
-            y[x0>0] <- -1/sc + (1+k)/(1+k*x1)*x1/sc
-            return(y)
+         y <- x*0
+         ind <- x>(mu-beta/xi)
+         y[ind] <- -1/beta - 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
-            x0 <- (x-tr)/sc
-            x1 <- x0[x0>0]
-            y[x0>0] <- log(1+k*x1)/k^2 - (1/k+1)*x1/(1+k*x1)
-            return(y)
+         y <- x*0
+         ind <- x>(mu-beta/xi)
+         y[ind]<- 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!
         z1 <- E(distribution, fun=Lambda1)
@@ -330,54 +107,33 @@
 
     ## Fisher Information matrix as a function of parameters
     FisherInfo.fct <- function(param) {
-        sc <- force(main(param)[1])
-        k <- force(main(param)[2])
-#        tr <- force(fixed(param)[1])
-#        fct <- L2deriv.fct(param)
-#        P2 <-  GPareto(loc = tr, scale = sc, shape = k)
-        E11 <- sc^-2
-        E12 <- (sc*(1+k))^-1
-        E22 <- 2/(1+k)
-        mat <- PosSemDefSymmMatrix(matrix(c(E11,E12,E12,E22)/(1+2*k),2,2))
-        dimnames(mat) <- list(scaleshapename,scaleshapename)
-        return(mat)
+        beta <- force(main(param)[1])
+        xi <- force(main(param)[2])
+        mu <- force(fixed(param)[1])
+        fct <- L2deriv.fct(param)
+        P <-  GEV(loc = mu, scale = beta, shape = xi)
+        E11 <- E(P,function(x)fct[[1]](x)*fct[[1]](x))
+        E12 <- E(P,function(x)fct[[1]](x)*fct[[2]](x))
+        E22 <- E(P,function(x)fct[[2]](x)*fct[[2]](x))
+        return(PosSemDefSymmMatrix(matrix(c(E11,E12,E12,E22),2,2)))
     }
 
     FisherInfo <- FisherInfo.fct(param)
-    name <- "Generalized Pareto Family"
+    name <- "Generalized Extreme Value Family with positive shape parameter: Frechet Family"
 
     ## initializing the GPareto family with components of L2-family
-    L2Fam <- new("GParetoFamily")
-    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()))
-
-    L2Fam at fam.call <- substitute(GParetoFamily(loc = loc0, scale = scale0,
-                                 shape = shape0, of.interest = of.interest0,
-                                 p = p0, N = N0, trafo = trafo0,
-                                 withPos = withPos0),
-                         list(loc0 = loc, scale0 = scale, shape0 = shape,
-                              of.interest0 = of.interest, p0 = p, N0 = N,
-                              trafo0 = trafo, withPos0 = withPos))
-
-    L2Fam at LogDeriv <- function(x) (shape+1)/(shape*(scale+(x-loc)))
-    L2Fam at L2deriv <- L2deriv
-
-    L2Fam at L2derivDistr <- imageDistr(RandVar = L2deriv, distr = distribution)
-
-
-    return(L2Fam)
+    res <- L2ParamFamily(name = name, param = param, 
+                         distribution = distribution, 
+                         L2deriv.fct = L2deriv.fct, 
+                         FisherInfo.fct = FisherInfo.fct,
+                         FisherInfo = FisherInfo,
+                         startPar = startPar,
+                         makeOKPar = makeOKPar,
+                         modifyParam = modifyPar,
+                         .returnClsName = "GEVFamily")
+    f.call <- substitute(GEVFamily(loc = loc0, scale = scale0, shape = shape0,trafo = trafo0), 
+                         list(loc0 = loc, scale0 = scale, shape0 = shape,trafo0 = trafo))
+    res at fam.call <- f.call
+    return(res)
 }
 

Modified: branches/robast-0.9/pkg/RobExtremes/R/drawPointsIC.R
===================================================================
--- branches/robast-0.9/pkg/RobExtremes/R/drawPointsIC.R	2012-10-31 08:23:45 UTC (rev 524)
+++ branches/robast-0.9/pkg/RobExtremes/R/drawPointsIC.R	2012-11-05 07:50:51 UTC (rev 525)
@@ -17,9 +17,9 @@
 ##evaluation of IC on observations
 IC0 = evalIC(IC,as.matrix(x)) 
 
-if (is.null(col)) col = rgb(192,192,192,maxColorValue=255,alpha=0.8) 
-if (is.null(bg)) bg = c("#70707050")
-if (is.null(bg)) xlab = "x"
+if (is.null(col)) col = rep(rgb(192,192,192,maxColorValue=255,alpha=0.8), length= dim(IC0)[1])
+if (is.null(bg)) bg = rep(c("#70707050"), length = dim(IC0)[1])
+if (is.null(xlab)) xlab = "x"
 if (is.null(ylab)) {ylab = numeric()
 for (i in 1:dim(IC0)[1]) ylab[i] = paste("IC of",names(L2fam at param@main[i]),"parameter")}
 
@@ -34,9 +34,9 @@
 ##plotting observations onto IC
 points(x,IC0[i,]
 ,pch = 21
-,cex = abs(IC0[i,]*sqrt(2))
-,col = col
-,bg = bg
+,cex = (colSums(IC0^2))^(1/2)*sqrt(2)
+,col = col[i]
+,bg = bg[i]
 )
 
 grid()
@@ -44,9 +44,37 @@
 }
 
 ##Example
-L2fam = GParetoFamily(shape = 0.7)
-IC = optIC(L2fam,risk=asCov())
-x = q(L2fam at distribution)(seq(0,1,length=50))
 
-##with probit transformation
-drawPointsIC(x,IC)
+##GPD family
+L2fam.GPD = GParetoFamily(shape = 0.7)
+IC.GPD.MLE = optIC(L2fam.GPD,risk=asCov())
+#IC.GPD.OMSE = optIC(L2fam.GPD,risk=asMSE())
+
+xgpd = q(L2fam.GPD at distribution)(seq(0,1,length=50))
+
+# par(mfcol=c(1,2))
+
+pdf("IC_GPD.pdf")
+drawPointsIC(xgpd,IC.GPD.MLE,col=c(rgb(7,154,121,maxColorValue=255,alpha=0.8),
+rgb(7,154,121,maxColorValue=255,alpha=0.8)),bg=c(c("#00640040"),c("#00640040")))
+
+# drawPointsIC(xgpd,IC.GPD.OMSE,col=c(rgb(7,154,121,maxColorValue=255,alpha=0.8),
+# rgb(7,154,121,maxColorValue=255,alpha=0.8)),bg=c(c("#ff100040"),c("#ff100040")))
+dev.off()
+
+##GEV family
+L2fam.GEV = GEVFamily(shape = 0.7)
+IC.GEV.MLE = optIC(L2fam.GEV,risk=asCov())
+#IC.GEV.OMSE = optIC(L2fam.GEV,risk=asMSE())
+
+xgev = q(L2fam.GEV at distribution)(seq(0,1,length=50))
+
+pdf("IC_GEV.pdf")
+# par(mfcol=c(1,2))
+drawPointsIC(xgev,IC.GEV.MLE,col=c(rgb(7,154,121,maxColorValue=255,alpha=0.8),
+rgb(7,154,121,maxColorValue=255,alpha=0.8)),bg=c(c("#00640040"),c("#00640040")))
+
+# drawPointsIC(xgpd,IC.GEV.OMSE,col=c(rgb(7,154,121,maxColorValue=255,alpha=0.8),
+# rgb(7,154,121,maxColorValue=255,alpha=0.8)),bg=c(c("#ff100040"),c("#ff100040")))
+dev.off()
+

Modified: branches/robast-0.9/pkg/RobExtremes/R/plotOutlyingness.R
===================================================================
--- branches/robast-0.9/pkg/RobExtremes/R/plotOutlyingness.R	2012-10-31 08:23:45 UTC (rev 524)
+++ branches/robast-0.9/pkg/RobExtremes/R/plotOutlyingness.R	2012-11-05 07:50:51 UTC (rev 525)
@@ -14,7 +14,7 @@
 ##@fam - parameter family
 ##@alpha - confidence level for quantile
 #
-plotOutlyingness = function(x,X=GPareto(),fam=GParetoFamily(),alpha=0.95){
+plotOutlyingness = function(x,alpha=0.99,X=GPareto(),fam=GParetoFamily()){
 
   ##logarithmic representation (for distributions with positive support)
   fam at distribution = log(fam at distribution)
@@ -39,8 +39,8 @@
  ,cex.idn = 1.7
  ,col.cutoff = rgb(202,202,202,maxColorValue=255) 
  ,offset = 0
- ,cutoff.quantile.y = alpha
- ,cutoff.quantile.x = alpha
+ ,cutoff.quantile.y = 0.99
+ ,cutoff.quantile.x = 0.99
  ,cutoff.x = cutoff()
  ,cutoff.y = cutoff.sememp()
  ,robCov.x = TRUE
@@ -59,3 +59,9 @@
  ,ylab="Mahalanobis distance"
 )
 }
+
+##Example
+X= GPareto()
+fam = GParetoFamily()
+x = r(X)(1000)
+plotOutlyingness(x,alpha=0.9,X=X,fam=fam)


Property changes on: branches/robast-0.9/pkg/RobExtremes/tests
___________________________________________________________________
Added: svn:ignore
   + Tests


Modified: branches/robast-0.9/pkg/RobExtremes/tests/TestSuite/TestExpectation.R
===================================================================
--- branches/robast-0.9/pkg/RobExtremes/tests/TestSuite/TestExpectation.R	2012-10-31 08:23:45 UTC (rev 524)
+++ branches/robast-0.9/pkg/RobExtremes/tests/TestSuite/TestExpectation.R	2012-11-05 07:50:51 UTC (rev 525)
@@ -8,8 +8,9 @@
 #author, take precedence over the dummy definitions provided by the
 #RUnit package and are called once for every test case identified.
 
- .setUp()
- {  
+ .setUp{ 
+
+ 
    ##expectation of Pareto distributed random variable
    expectation.Pareto = function(shape0=1,Min0=1){
     X = Pareto(shape=shape0,Min=Min0)
@@ -19,7 +20,7 @@
    test.expectationPareto = function(){
     checkEquals(expectation.Pareto(1,1), Inf)
     checkEquals(expectation.Pareto(2,1), 0)
-     }    
+   }    
   
 #   test.HTMLInfo.Pareto = function(){
 #    track <- tracker()
@@ -32,13 +33,14 @@
 #    resTrack <- track$getTrackInfo()
 #    }
 # 
-  }
+#  }
 
 #  .tearDown(){
 #   ##create HTML sites in folder ./results for all inspect calls
 #   printHTML.trackInfo(resTrack,"TestSuite/TestExpectation")
 #   }
 
+}
 
 #  
 # Beispiele

Modified: branches/robast-0.9/pkg/RobExtremes/tests/TestSuite.R
===================================================================
--- branches/robast-0.9/pkg/RobExtremes/tests/TestSuite.R	2012-10-31 08:23:45 UTC (rev 524)
+++ branches/robast-0.9/pkg/RobExtremes/tests/TestSuite.R	2012-11-05 07:50:51 UTC (rev 525)
@@ -21,7 +21,7 @@
 
 ## run test suite
 myTestSuite <- defineTestSuite("RUnit Tests for RobExtremes",dirs="./TestSuite",
-testFileRegexp = "^Test.+\\R$",testFuncRegexp = "^test.+")
+testFileRegexp = "^Test.+R$",testFuncRegexp = "^test.+")
 testResult <- runTestSuite(myTestSuite)
 
 ## is test suite valid?
@@ -34,4 +34,12 @@
 ## to standard out:
 printTextProtocol(testResult, "testResult.txt",showDetails = TRUE)
 printHTMLProtocol(testResult, "testResult.html")
- 
+
+
+#  
+# runTestFile(absFileName, useOwnErrorHandler = TRUE,
+#                  testFuncRegexp = "^test.+",
+#                  rngKind = "Marsaglia-Multicarry",
+#                  rngNormalKind = "Kinderman-Ramage",
+#                  verbose = getOption("RUnit")$verbose)
+

Modified: branches/robast-0.9/pkg/RobExtremes/tests/testResult.html
===================================================================
--- branches/robast-0.9/pkg/RobExtremes/tests/testResult.html	2012-10-31 08:23:45 UTC (rev 524)
+++ branches/robast-0.9/pkg/RobExtremes/tests/testResult.html	2012-11-05 07:50:51 UTC (rev 525)
@@ -1,11 +1,11 @@
 <!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN"
 "http://www.w3.org/TR/html4/transitional.dtd">
-<html><head><title>RUNIT TEST PROTOCOL--Thu Sep  6 16:40:17 2012</title>
+<html><head><title>RUNIT TEST PROTOCOL--Tue Sep 18 10:02:30 2012</title>
 </head>
-<body><h1 TRUE>RUNIT TEST PROTOCOL--Thu Sep  6 16:40:17 2012</h1>
-<p>Number of test functions: 2</p>
+<body><h1 TRUE>RUNIT TEST PROTOCOL--Tue Sep 18 10:02:30 2012</h1>
+<p>Number of test functions: 1</p>
 <p>Number of errors: 0</p>
-<p>Number of failures: 0</p>
+<p style=color:red>Number of failures: 1</p>
 <hr>
 <h3 TRUE>1 Test suite</h3>
 <table border="1" width="60%" >
@@ -15,15 +15,29 @@
 <th width="20%">Failures</th>
 </tr>
 <tr><td><a href="#RUnit Tests for RobExtremes">RUnit Tests for RobExtremes</a></td>
-<td>2</td>
+<td>1</td>
 <td>0</td>
-<td>0</td>
+<td bgcolor="red">1</td>
 </tr>
 </table>
 <hr>
+<h3 TRUE>Failures</h3>
+<table border="1" width="100%" >
+<tr><th width="30%">Test suite : test function</th>
+<th width="70%">message</th>
+</tr>
+<tr><td><a href="#RUnit Tests for RobExtremes_._TestSuite_TestExpectation.R_test.expectationPareto">RUnit Tests for RobExtremes : test.expectationPareto</a></td>
+<td>Error in checkEquals(expectation.Pareto(2, 1), 0) : 
+  Mean relative difference: 1
+</td>
+</tr>
+</table>
+<hr>
 <h3 TRUE>Details</h3>
 <p><a name="RUnit Tests for RobExtremes"><h5 TRUE>Test Suite: RUnit Tests for RobExtremes</h5>
-</a>Test function regexp: ^test.+<br/>Test file regexp: ^Test.+\R$<br/>Involved directory:<br/>./TestSuite<br/><ul><li><a href="./TestSuite/TestExpectation.R">Test file: TestExpectation.R</a><ul><li><a name="RUnit Tests for RobExtremes_._TestSuite_TestExpectation.R_test.expectationPareto">test.expectationPareto: (1 checks) ... OK (0 seconds)<br/></a></li><li><a name="RUnit Tests for RobExtremes_._TestSuite_TestExpectation.R_test.HTMLInfo.Pareto">test.HTMLInfo.Pareto: (0 checks) ... OK (0.01 seconds)<br/></a></li></ul></li></li></li></li></li></li></li></li></li></li></li></li></li></ul><hr>
+</a>Test function regexp: ^test.+<br/>Test file regexp: ^Test.+R$<br/>Involved directory:<br/>./TestSuite<br/><ul><li><a href="./TestSuite/TestExpectation.R">Test file: TestExpectation.R</a><ul><li><a name="RUnit Tests for RobExtremes_._TestSuite_TestExpectation.R_test.expectationPareto"><u style=color:red>test.expectationPareto: FAILURE !! (check number 2)  </u></a>Error in checkEquals(expectation.Pareto(2, 1), 0) : 
+  Mean relative difference: 1
+<br/></li></ul></li></li></li></li></li></li></li></li></li></li></li></li></li></ul><hr>
 <table border="0" width="80%" >
 <tr><th>Name</th>
 <th>Value</th>

Modified: branches/robast-0.9/pkg/RobExtremes/tests/testResult.txt
===================================================================
--- branches/robast-0.9/pkg/RobExtremes/tests/testResult.txt	2012-10-31 08:23:45 UTC (rev 524)
+++ branches/robast-0.9/pkg/RobExtremes/tests/testResult.txt	2012-11-05 07:50:51 UTC (rev 525)
@@ -1,12 +1,14 @@
-RUNIT TEST PROTOCOL -- Thu Sep  6 16:40:17 2012 
+RUNIT TEST PROTOCOL -- Tue Sep 18 10:02:30 2012 
 *********************************************** 
-Number of test functions: 2 
+Number of test functions: 1 
 Number of errors: 0 
-Number of failures: 0 
+Number of failures: 1 
 
  
 1 Test Suite : 
-RUnit Tests for RobExtremes - 2 test functions, 0 errors, 0 failures
+RUnit Tests for RobExtremes - 1 test function, 0 errors, 1 failure
+FAILURE in test.expectationPareto: Error in checkEquals(expectation.Pareto(2, 1), 0) : 
+  Mean relative difference: 1
 
 
 
@@ -14,10 +16,11 @@
 *************************** 
 Test Suite: RUnit Tests for RobExtremes 
 Test function regexp: ^test.+ 
-Test file regexp: ^Test.+\R$ 
+Test file regexp: ^Test.+R$ 
 Involved directory: 
 ./TestSuite 
 --------------------------- 
 Test file: ./TestSuite/TestExpectation.R 
-test.expectationPareto: (1 checks) ... OK (0 seconds)
-test.HTMLInfo.Pareto: (0 checks) ... OK (0.01 seconds)
+test.expectationPareto: FAILURE !! (check number 2)
+Error in checkEquals(expectation.Pareto(2, 1), 0) : 
+  Mean relative difference: 1



More information about the Robast-commits mailing list