[Robast-commits] r547 - branches/robast-0.9/pkg/RobExtremes/R

noreply at r-forge.r-project.org noreply at r-forge.r-project.org
Mon Jan 21 19:55:57 CET 2013


Author: ruckdeschel
Date: 2013-01-21 19:55:57 +0100 (Mon, 21 Jan 2013)
New Revision: 547

Modified:
   branches/robast-0.9/pkg/RobExtremes/R/GEVFamily.R
   branches/robast-0.9/pkg/RobExtremes/R/GParetoFamily.R
   branches/robast-0.9/pkg/RobExtremes/R/LDEstimator.R
   branches/robast-0.9/pkg/RobExtremes/R/PickandsEstimator.R
   branches/robast-0.9/pkg/RobExtremes/R/WeibullFamily.R
Log:
RobExtremes: modularized  generating functions GEVFamily, GParetoFamily, and WeibullFamily [ as to defining trafo, tau, Dtau ]; deleted stuff from "simple" estimators (LDE, Pickands, QuantileBCC) which is delegated to Estimator in distrMod() (i.e., asVar-calculations)

Modified: branches/robast-0.9/pkg/RobExtremes/R/GEVFamily.R
===================================================================
--- branches/robast-0.9/pkg/RobExtremes/R/GEVFamily.R	2013-01-19 02:48:08 UTC (rev 546)
+++ branches/robast-0.9/pkg/RobExtremes/R/GEVFamily.R	2013-01-21 18:55:57 UTC (rev 547)
@@ -4,33 +4,10 @@
 ##
 ################################
 
+### some reusable blocks of code (to avoid redundancy):
 
-## 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 (any(param[2] <= tol))
-                 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){
+### pretreatment of of.interest
+.pretreat.of.interest <- function(of.interest,trafo){
     if(is.null(trafo)){
         of.interest <- unique(of.interest)
         if(length(of.interest) > 2)
@@ -60,62 +37,11 @@
             of.interest[1] <- "expected shortfall"
         }
     }
-    theta <- c(loc, scale, shape)
+  return(of.interest)
+}
 
-    ##symmetry
-    distrSymm <- NoSymmetry()
-
-    ## parameters
-    names(theta) <- c("loc", "scale", "shape")
-    scaleshapename <- c("scale", "shape")
-
-
-    if(!is.null(p)){
-       btq <- substitute({ q <- loc0 + theta[1]*((-log(1-p0))^(-theta[2])-1)/theta[2]
-                           names(q) <- "quantile"
-                           }, list(loc0 = loc, p0 = p))
-
-       bDq <- substitute({ scale <- theta[1];  shape <- theta[2]
-                        D1 <- ((-log(1-p0))^(-shape)-1)/shape
-                        D2 <- -scale/shape*(D1 + log(-log(1-p0))*(-log(1-p0))^(-shape))
-                        D <- t(c(D1, D2))
-                        rownames(D) <- "quantile"; colnames(D) <- NULL
-                        D }, list(p0 = p))
-       btes <- substitute({ if(theta[2]>=1L) 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 <- theta[1]*D1/theta[2]
-                            D22 <- theta[1]*dd/(1-p0)/theta[2]
-                            D2 <- -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) el <- NA else{
-                            el <- N0*(loc0+theta[1]*gamma(1-theta[2])/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]
-                            D1 <- N0*gamma(1-shape)/shape
-                            D2 <- -N0*theta[1]*digamma(1-theta[2])/theta[2]-
-                                   D1*scale/(1-shape)}
-                            D <- t(c(D1, D2))
-                            rownames(D) <- "expected loss"
-                            colnames(D) <- NULL
-                            D }, list(loc0 = loc, N0 = N))
-    }
+.define.tau.Dtau <- function( trafo, of.interest, btq, bDq, btes,
+                     bDes, btel, bDel, p, N){
     if(is.null(trafo)){
         tau <- NULL
         if("scale" %in% of.interest){
@@ -185,8 +111,100 @@
     }else{
         if(is.matrix(trafo) & nrow(trafo) > 2) stop("number of rows of 'trafo' > 2")
     }
+    return(list(trafo = trafo, tau = tau, Dtau = Dtau))
+}
 
+## 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 (any(param[2] <= tol))
+                 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){
+    theta <- c(loc, scale, shape)
+
+    of.interest <- .pretreat.of.interest(of.interest,trafo)
+
+    ##symmetry
+    distrSymm <- NoSymmetry()
+
+    ## parameters
+    names(theta) <- c("loc", "scale", "shape")
+    scaleshapename <- c("scale", "shape")
+
+
+    btq <- bDq <- btes <- bDes <- btel <- bDel <- NULL
+    if(!is.null(p)){
+       btq <- substitute({ q <- loc0 + theta[1]*((-log(1-p0))^(-theta[2])-1)/theta[2]
+                           names(q) <- "quantile"
+                           }, list(loc0 = loc, p0 = p))
+
+       bDq <- substitute({ scale <- theta[1];  shape <- theta[2]
+                        D1 <- ((-log(1-p0))^(-shape)-1)/shape
+                        D2 <- -scale/shape*(D1 + log(-log(1-p0))*(-log(1-p0))^(-shape))
+                        D <- t(c(D1, D2))
+                        rownames(D) <- "quantile"; colnames(D) <- NULL
+                        D }, list(p0 = p))
+       btes <- substitute({ if(theta[2]>=1L) 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 <- theta[1]*D1/theta[2]
+                            D22 <- theta[1]*dd/(1-p0)/theta[2]
+                            D2 <- -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) el <- NA else{
+                            el <- N0*(loc0+theta[1]*gamma(1-theta[2])/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]
+                            D1 <- N0*gamma(1-shape)/shape
+                            D2 <- -N0*theta[1]*digamma(1-theta[2])/theta[2]-
+                                   D1*scale/(1-shape)}
+                            D <- t(c(D1, D2))
+                            rownames(D) <- "expected loss"
+                            colnames(D) <- NULL
+                            D }, list(loc0 = loc, N0 = N))
+    }
+
+    def <- .define.tau.Dtau( trafo, of.interest, btq, bDq, btes,
+                             bDes, btel, bDel, p, N)
+    trafo <- def$trafo; tau <- def$tau; Dtau <- def$Dtau
+
+####
     param <- ParamFamParameter(name = "theta", main = c(theta[2],theta[3]),
                                fixed = theta[1],
                                trafo = trafo, withPosRestr = withPos,

Modified: branches/robast-0.9/pkg/RobExtremes/R/GParetoFamily.R
===================================================================
--- branches/robast-0.9/pkg/RobExtremes/R/GParetoFamily.R	2013-01-19 02:48:08 UTC (rev 546)
+++ branches/robast-0.9/pkg/RobExtremes/R/GParetoFamily.R	2013-01-21 18:55:57 UTC (rev 547)
@@ -38,37 +38,11 @@
                           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"
-        }
-    }
     theta <- c(loc, scale, shape)
 
+    of.interest <- .pretreat.of.interest(of.interest,trafo)
+             ## code .pretreat.of.interest in GEV.family.R
+
     ##symmetry
     distrSymm <- NoSymmetry()
 
@@ -76,6 +50,7 @@
     names(theta) <- c("loc", "scale", "shape")
     scaleshapename <- c("scale", "shape")
 
+    btq <- bDq <- btes <- bDes <- btel <- bDel <- NULL
     if(!is.null(p)){
        btq <- substitute({ q <- loc0 + theta[1]*((1-p0)^(-theta[2])-1)/theta[2]
                            names(q) <- "quantile"
@@ -119,75 +94,12 @@
                             colnames(D) <- NULL
                             D }, list(loc0 = loc, N0 = N))
     }
-    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) <- btq
-                Dtau <- function(theta){ };body(Dtau) <- bDq
-            }else{
-                tau1 <- tau
-                tau <- function(theta){ }
-                body(tau) <- substitute({ btq0; c(tau0(theta), q) },
-                                        list(btq0=btq, tau0 = tau1))
-                Dtau1 <- Dtau
-                Dtau <- function(theta){}
-                body(Dtau) <- substitute({ bDq0; rbind(Dtau0(theta), D) },
-                                         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; c(tau0(theta), es) },
-                                        list(tau0 = tau1, btes0=btes))
-                Dtau1 <- Dtau
-                Dtau <- function(theta){}
-                body(Dtau) <- substitute({ bDes0; rbind(Dtau0(theta), D) },
-                                         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; c(tau0(theta), el) },
-                                        list(tau0 = tau1, btel0=btel))
-                Dtau1 <- Dtau
-                Dtau <- function(theta){}
-                body(Dtau) <- substitute({ bDel0; rbind(Dtau0(theta), D) },
-                                         list(Dtau0 = Dtau1, bDel0=bDel))
-            }
-        }
-        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")
-    }
+
+    def <- .define.tau.Dtau( trafo, of.interest, btq, bDq, btes,
+                             bDes, btel, bDel, p, N)
+           # code .define.tau.Dtau is in file GEVFamily.R
+    trafo <- def$trafo; tau <- def$tau; Dtau <- def$Dtau
+
     param <- ParamFamParameter(name = "theta", main = c(theta[2],theta[3]),
                                fixed = theta[1],
                                trafo = trafo, withPosRestr = withPos,

Modified: branches/robast-0.9/pkg/RobExtremes/R/LDEstimator.R
===================================================================
--- branches/robast-0.9/pkg/RobExtremes/R/LDEstimator.R	2013-01-19 02:48:08 UTC (rev 546)
+++ branches/robast-0.9/pkg/RobExtremes/R/LDEstimator.R	2013-01-21 18:55:57 UTC (rev 547)
@@ -99,35 +99,6 @@
                       na.rm = na.rm.0, ...,
                       ParamFamily = ParamFamily)
 
-    cat("\n asvar",estimate at asvar,"\n")
-
-
-##->
-if(FALSE){
-    if(missing(asvar)) asvar <- NULL
-
-    if((is.null(asvar))&&(!missing(asvar.fct))&&(!is.null(asvar.fct)))
-          asvar <- asvar.fct(ParamFamily, estimate, ...)
-
-    estimate at untransformed.asvar <- asvar
-
-    l.e <- length(estimate at untransformed.estimate)
-    idx <- NULL
-    idm <- 1:l.e
-    if(!is.null(nuis.idx))
-        {idx <- nuis.idx
-         idm <- idm[-idx]
-         mat <- diag(length(idm))}
-
-    if(!.isUnitMatrix(estimate at trafo$mat)){
-       estimate at estimate <- estimate at trafo$fct(estimate)
-       
-       if(!is.null(asvar))
-           estimate at asvar <- estimate at trafo$mat%*%asvar[idm,idm]%*%t(estimate at trafo$mat)
-    }
-
-}
-##<-
     estimate at estimate.call <- es.call
 
     if(missing(Infos))

Modified: branches/robast-0.9/pkg/RobExtremes/R/PickandsEstimator.R
===================================================================
--- branches/robast-0.9/pkg/RobExtremes/R/PickandsEstimator.R	2013-01-19 02:48:08 UTC (rev 546)
+++ branches/robast-0.9/pkg/RobExtremes/R/PickandsEstimator.R	2013-01-21 18:55:57 UTC (rev 547)
@@ -66,27 +66,6 @@
                           nuis.idx = nuis.idx.0, trafo = trafo.0,
                           fixed = fixed.0, na.rm = na.rm.0, ...,
                           ParamFamily = ParamFamily)
-##->
-if(FALSE){
-    estimate at untransformed.asvar <- asvar(estimate)
-    estimate at asvar <- asvar
-
-
-    l.e <- length(estimate at untransformed.estimate)
-    idx <- NULL
-    idm <- 1:l.e
-    if(!is.null(nuis.idx))
-        {idx <- nuis.idx
-         idm <- idm[-idx]
-         mat <- diag(length(idm))}
-
-    if(!.isUnitMatrix(estimate at trafo$mat)){
-       estimate at estimate <- estimate at trafo$fct(estimate)
-       if(!is.null(asvar))
-           estimate at asvar <- estimate at trafo$mat%*%asvar[idm,idm]%*%t(estimate at trafo$mat)
-    }
-}
-## <-
     estimate at estimate.call <- es.call
 
     if(missing(Infos))

Modified: branches/robast-0.9/pkg/RobExtremes/R/WeibullFamily.R
===================================================================
--- branches/robast-0.9/pkg/RobExtremes/R/WeibullFamily.R	2013-01-19 02:48:08 UTC (rev 546)
+++ branches/robast-0.9/pkg/RobExtremes/R/WeibullFamily.R	2013-01-21 18:55:57 UTC (rev 547)
@@ -37,37 +37,11 @@
                           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"
-        }
-    }
     theta <- c(scale, shape)
 
+    of.interest <- .pretreat.of.interest(of.interest,trafo)
+             ## code .pretreat.of.interest in GEV.family.R
+
     ##symmetry
     distrSymm <- NoSymmetry()
 
@@ -75,6 +49,7 @@
     names(theta) <- c("scale", "shape")
     scaleshapename <- c("scale", "shape")
 
+    btq <- bDq <- btes <- bDes <- btel <- bDel <- NULL
     if(!is.null(p)){
        btq <- substitute({ q <- theta[1]*(-log(1-p0))^(1/theta[2])
                            names(q) <- "quantile"
@@ -120,75 +95,10 @@
                             D }, list(loc0 = loc, N0 = N))
     }
 
-    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) <- btq
-                Dtau <- function(theta){ };body(Dtau) <- bDq
-            }else{
-                tau1 <- tau
-                tau <- function(theta){ }
-                body(tau) <- substitute({ btq0; c(tau0(theta), q) },
-                                        list(btq0=btq, tau0 = tau1))
-                Dtau1 <- Dtau
-                Dtau <- function(theta){}
-                body(Dtau) <- substitute({ bDq0; rbind(Dtau0(theta), D) },
-                                         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; c(tau0(theta), es) },
-                                        list(tau0 = tau1, btes0=btes))
-                Dtau1 <- Dtau
-                Dtau <- function(theta){}
-                body(Dtau) <- substitute({ bDes0; rbind(Dtau0(theta), D) },
-                                         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; c(tau0(theta), el) },
-                                        list(tau0 = tau1, btel0=btel))
-                Dtau1 <- Dtau
-                Dtau <- function(theta){}
-                body(Dtau) <- substitute({ bDel0; rbind(Dtau0(theta), D) },
-                                         list(Dtau0 = Dtau1, bDel0=bDel))
-            }
-        }
-        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")
-    }
+    def <- .define.tau.Dtau( trafo, of.interest, btq, bDq, btes,
+                             bDes, btel, bDel, p, N)
+           # code .define.tau.Dtau is in file GEVFamily.R
+    trafo <- def$trafo; tau <- def$tau; Dtau <- def$Dtau
 
 
     param <- ParamFamParameter(name = "theta", main = c(theta[1],theta[2]),



More information about the Robast-commits mailing list