[Distr-commits] r1232 - in branches/distr-2.8/pkg/distrMod: . R man

noreply at r-forge.r-project.org noreply at r-forge.r-project.org
Mon Jul 30 19:47:38 CEST 2018


Author: ruckdeschel
Date: 2018-07-30 19:47:38 +0200 (Mon, 30 Jul 2018)
New Revision: 1232

Modified:
   branches/distr-2.8/pkg/distrMod/NAMESPACE
   branches/distr-2.8/pkg/distrMod/R/AllClass.R
   branches/distr-2.8/pkg/distrMod/R/AllGeneric.R
   branches/distr-2.8/pkg/distrMod/R/MCEstimate.R
   branches/distr-2.8/pkg/distrMod/R/internalMleCalc.R
   branches/distr-2.8/pkg/distrMod/R/mleCalc-methods.R
   branches/distr-2.8/pkg/distrMod/man/MCEstimate-class.Rd
   branches/distr-2.8/pkg/distrMod/man/meRes.Rd
Log:
[distrMod] _only_ in branch 2.8: 

(a) "cleverer parsing of dots" in M[L,C,D]Estimator:

Triggered by an error observed by Kornelius Rohmeyer, we now parse the dots (...) argument 
of M[L,C,D]Estimator a  little more closely and filter out obvious clashes: We filter out taboo arguments 
(causing clashes with formals in optim/optimize and calling arguments); from the remaining ones only 
[exactly] _named_ arguments of the optimizer (optim/optimize) and those matching either exactly arguments 
of the criterion or all remaining ones (if ... is a formal of the criterion) are passed on: 
This way, one can, if desired, do evaluation-wise validity checks in the criterion. 

(b) for diagnostic purposes, MCEstimate-class gains a slot optimReturn (of class "ANY" and filled by
NULL by default) which is filled by the return value of the optimizer in "mceCalc" -- it has a
corresponding accessor

As this is more for internal purposes, example code is in lines ll 256--294 in mleCalc-methods.R
(wrapped in a if(FALSE){ }).  


Modified: branches/distr-2.8/pkg/distrMod/NAMESPACE
===================================================================
--- branches/distr-2.8/pkg/distrMod/NAMESPACE	2018-07-30 14:47:31 UTC (rev 1231)
+++ branches/distr-2.8/pkg/distrMod/NAMESPACE	2018-07-30 17:47:38 UTC (rev 1232)
@@ -60,7 +60,7 @@
               "name.estimate", "trafo.estimate", "nuisance.estimate", 
               "fixed.estimate", "Infos", "Infos<-", "addInfo<-",
               "criterion", "criterion<-", "criterion.fct", "method",
-              "samplesize", "asvar", "asvar<-", "optimwarn",
+              "samplesize", "asvar", "asvar<-", "optimwarn", "optimReturn",
 			  "withPosRestr", "withPosRestr<-")
 exportMethods("untransformed.estimate", "untransformed.asvar")
 exportMethods("confint")

Modified: branches/distr-2.8/pkg/distrMod/R/AllClass.R
===================================================================
--- branches/distr-2.8/pkg/distrMod/R/AllClass.R	2018-07-30 14:47:31 UTC (rev 1231)
+++ branches/distr-2.8/pkg/distrMod/R/AllClass.R	2018-07-30 17:47:38 UTC (rev 1232)
@@ -465,7 +465,8 @@
                         criterion.fct = "function",
                         method = "character",
                         optimwarn = "character",
-                        startPar = "ANY"),
+                        startPar = "ANY",
+                        optimReturn = "ANY"),
          prototype(name = "Minimum criterion estimate",
                    estimate = numeric(0),
                    samplesize = numeric(0),
@@ -482,7 +483,8 @@
                                       list(fval = x, mat = matrix(1))},
                                 mat = matrix(1)),
                    optimwarn = "",
-                   startPar = NULL
+                   startPar = NULL,
+                   optimReturn = NULL
                    ),
          contains = "Estimate")
 

Modified: branches/distr-2.8/pkg/distrMod/R/AllGeneric.R
===================================================================
--- branches/distr-2.8/pkg/distrMod/R/AllGeneric.R	2018-07-30 14:47:31 UTC (rev 1231)
+++ branches/distr-2.8/pkg/distrMod/R/AllGeneric.R	2018-07-30 17:47:38 UTC (rev 1232)
@@ -242,6 +242,9 @@
 if(!isGeneric("startPar")){
     setGeneric("startPar", function(object, ...) standardGeneric("startPar"))
 }
+if(!isGeneric("optimReturn")){
+    setGeneric("optimReturn", function(object, ...) standardGeneric("optimReturn"))
+}
 if(!isGeneric("makeOKPar")){
     setGeneric("makeOKPar", function(object, ...) standardGeneric("makeOKPar"))
 }              

Modified: branches/distr-2.8/pkg/distrMod/R/MCEstimate.R
===================================================================
--- branches/distr-2.8/pkg/distrMod/R/MCEstimate.R	2018-07-30 14:47:31 UTC (rev 1231)
+++ branches/distr-2.8/pkg/distrMod/R/MCEstimate.R	2018-07-30 17:47:38 UTC (rev 1232)
@@ -3,6 +3,7 @@
 ###############################################################################
 
 
+setMethod("optimReturn", "MCEstimate", function(object) object at optimReturn)
 setMethod("optimwarn", "MCEstimate", function(object) object at optimwarn)
 setMethod("criterion", "MCEstimate", function(object) object at criterion)
 setMethod("criterion.fct", "MCEstimate", function(object) object at criterion.fct)

Modified: branches/distr-2.8/pkg/distrMod/R/internalMleCalc.R
===================================================================
--- branches/distr-2.8/pkg/distrMod/R/internalMleCalc.R	2018-07-30 14:47:31 UTC (rev 1231)
+++ branches/distr-2.8/pkg/distrMod/R/internalMleCalc.R	2018-07-30 17:47:38 UTC (rev 1232)
@@ -140,6 +140,7 @@
                   untransformed.asvar = untransformed.asvar,
                   criterion.fct = res$crit.fct, method = res$method,
                   fixed = fixed(param), optimwarn = res$warns,
+                  optimReturn = res$optReturn,
                   startPar = res$startPar)
     return(res.me)
 }

Modified: branches/distr-2.8/pkg/distrMod/R/mleCalc-methods.R
===================================================================
--- branches/distr-2.8/pkg/distrMod/R/mleCalc-methods.R	2018-07-30 14:47:31 UTC (rev 1231)
+++ branches/distr-2.8/pkg/distrMod/R/mleCalc-methods.R	2018-07-30 17:47:38 UTC (rev 1232)
@@ -14,12 +14,12 @@
 meRes <- function(x, estimate, criterion.value, param, crit.fct,
                   method = "explicit solution",
                   crit.name = "Maximum Likelihood", Infos, warns = "",
-                  startPar = NULL)
+                  startPar = NULL, optReturn = NULL)
         return(list(estimate = estimate, criterion = criterion.value,
                     param = param, crit.fct = crit.fct, method = method, 
                     crit.name = crit.name, Infos = Infos, 
                     samplesize = samplesize(x), warns = warns,
-                    startPar = startPar))
+                    startPar = startPar, optReturn = optReturn))
 
 get.criterion.fct <- function(theta, Data, ParamFam, criterion.ff, fun, ...){
 
@@ -104,14 +104,66 @@
         lnx <- length(nuisance(PFam))
         fixed <- fixed(PFam)
 
+## added 2018 07 30: parse the dots argument a little to filter
+## out arguments which could be of interest to optim/optimize resp. to criterion
+## consequence: MCEstimator can be called with (exactly) named arguments of
+## optimize/optim and with (exactly) named additional arguments (from position
+## 3 on) of the criterion function; in particular if check.validity is a formal.
+## All other arguments are not passed on / in particular if ... is not a formal
+## and if check.validity is not a formal of the criterion, it is simply ignored.
+
+#       mceCalcDots1 <- match.call(call = sys.call(sys.parent(1)),
+#                                 expand.dots = FALSE)$"..."
+       mceCalcDots <- match.call(expand.dots = FALSE)$"..."
+       filterDots <- function(dots){
+          if(length(dots)){
+               dotsOptIz <- NULL
+               nfmlsOptiz <- NULL
+
+               dotsNames <- names(dots)
+               if(length(param(PFam)) == 1){
+                  nfmlsOptIz <- names(formals(optimize))
+                  nOptProh <- c("f","interval")
+               }else{
+                  nfmlsOptIz <- names(formals(optim))
+                  nOptProh <- c("fn","par")
+               }
+               nfmlsOptIz <- nfmlsOptIz[! (nfmlsOptIz %in% nOptProh)]
+               if(length(nfmlsOptIz))
+                  dotsOptIz <- dots[dotsNames %in% nfmlsOptIz]
+               if(length(dotsOptIz)==0) dotsOptIz <- NULL
+
+               dotsRest <- dots[!(dotsNames %in% nfmlsOptIz)]
+               nfmlsCrit <- names(formals(criterion))
+               dotsForCrit <- NULL
+
+               if(length(dotsRest)&&length(nfmlsCrit)>2){
+                  dotsRestNames <- names(dotsRest)
+                  nmprohib <- c("Data", "theta", "ParamFamily",
+                                "criterionF", nfmlsOptIz)
+                  nmprohib <- nmprohib[nmprohib!="..."]
+                  nfmlsCrit <- nfmlsCrit[! (nfmlsCrit %in% nmprohib)]
+                  dotsRestNamesProhib <- dotsRestNames %in% nmprohib
+
+                  dotsRest <- dotsRest[!dotsRestNamesProhib]
+                  if(length(dotsRest)){
+                     dotsRestNames <- names(dotsRest)
+                     critL <- (dotsRestNames %in% nfmlsCrit)
+                     critL <- critL | ("..." %in% nfmlsCrit)
+                     dotsForCrit <- dotsRest[critL]
+                  }
+                  if(length(dotsForCrit)==0) dotsForCrit <- NULL
+               }
+               dotsForOpt <- c(dotsOptIz,dotsForCrit[!names(dotsForCrit)%in% nOptProh])
+               return(list(dotsForOpt=dotsForOpt, dotsCrit=dotsForCrit))
+          }else return(NULL)
+       }
+
+       dotsToPass <- do.call(filterDots, list(mceCalcDots))
        allwarns <- character(0)
        fun <- function(theta, Data, ParamFamily, criterionF, ...){
                vP <- TRUE
                if(validity.check) vP <- validParameter(ParamFamily, theta)
-               dots <- list(...)
-               dots$trafo <- NULL
-               dots$penalty <- NULL
-               dots$withBiasC <- NULL
                if(is.function(penalty)) penalty <- penalty(theta)
                if(!vP) {crit0 <- penalty; theta <- mO(theta)
                }else{
@@ -120,8 +172,7 @@
                                        names(nuisance(ParamFamily)))
                   else  names(theta) <- names(main(ParamFamily))
                   distr.new <- try(ParamFamily at modifyParam(theta), silent = TRUE)
-                  argList <- c(list(Data, distr.new), dots)
-                  argList$check.validity <- NULL
+                  argList <- c(list(Data, distr.new), ... )
                   if(withthetaPar) argList <- c(argList, list(thetaPar = theta))
                   if(is(distr.new,"try.error")){
                       crit0 <- penalty
@@ -148,20 +199,22 @@
                return(critP)}
 
     if(length(param(PFam)) == 1){
-        res <- optimize(f = fun, interval = startPar, Data = x,
-                      ParamFamily = PFam, criterionF = criterion, ...)
-        theta <- res$minimum
+        optres <- do.call(optimize, c(list(f = fun, interval = startPar, Data = x,
+                      ParamFamily = PFam, criterionF = criterion),
+                      dotsToPass$dotsForOpt))
+        theta <- optres$minimum
         names(theta) <- names(main(PFam))
-        crit <- res$objectiv
+        crit <- optres$objective
         method <- "optimize"
     }else{
         if(is(startPar,"Estimate")) startPar <- untransformed.estimate(startPar)
-        res <- optim(par = startPar, fn = fun, Data = x,
-                   ParamFamily = PFam, criterionF = criterion, ...)
-        theta <- as.numeric(res$par)
+        optres <- do.call(optim, c(list(par = startPar, fn = fun, Data = x,
+                   ParamFamily = PFam, criterionF = criterion),
+                   dotsToPass$dotsForOpt))
+        theta <- as.numeric(optres$par)
         names(theta) <- c(names(main(PFam)),names(nuisance(PFam)))
         method <- "optim"
-        crit <- res$value
+        crit <- optres$value
     }
 
     vP <- TRUE
@@ -175,6 +228,7 @@
     param <- .callParamFamParameter(PFam, theta, idx, nuis, fixed)
 
     fun2 <- function(theta, Data, ParamFamily, criterion, ...){
+               dotsTP <- filterDots(list(...))
                vP <- TRUE
                if(validity.check) vP <- validParameter(ParamFamily, theta)
                if(!vP) theta <- makeOKPar(ParamFamily)(theta)
@@ -183,17 +237,66 @@
                                        names(nuisance(ParamFamily)))
                else  names(theta) <- names(main(ParamFamily))
                distr.new <- ParamFamily at modifyParam(theta)
-               crit1 <- criterion(Data, distr.new, ...)
+               crit1 <- do.call(criterion, c(list(Data, distr.new),
+                                dotsToPass$dotsCrit))
                return(crit1)}
 
-    crit.fct <- get.criterion.fct(theta, Data = x, ParamFam = PFam, 
+    crit.fct <- get.criterion.fct(theta, Data = x, ParamFam = PFam,
                                    criterion.ff = criterion, fun2, ...)
-    
+
     return(meRes(x, theta, crit, param, crit.fct, method = method,
                  crit.name = crit.name, Infos = Infos, warns= allwarns,
-                 startPar = startPar))
+                 startPar = startPar, optReturn = optres))
            })
 
+
+##########################################################################
+# added 2018 07 30
+##########################################################################
+if(FALSE){
+require(distrMod)
+nF <- NormLocationScaleFamily()
+negLoglikelihood <- function(x, Distribution){
+  res <- -sum(log(Distribution at d(x)))
+  names(res) <- "Negative Log-Likelihood"
+  return(res)
+}
+negLoglikelihood2 <- function(x, Distribution, ...){
+  dots <- list(...)
+  print(dots)
+  res <- -sum(log(Distribution at d(x)))
+  names(res) <- "Negative Log-Likelihood"
+  return(res)
+}
+negLoglikelihood3 <- function(x, Distribution, check.validity, fn=3){
+  print(c(chk=check.validity))
+  print(c(fn=fn))
+  res <- -sum(log(Distribution at d(x)))
+  names(res) <- "Negative Log-Likelihood"
+  return(res)
+}
+set.seed(123)
+x <- rnorm(10)
+MCEstimator(x = x, ParamFamily = nF, criterion = negLoglikelihood)
+re <- MCEstimator(x = x, ParamFamily = nF, criterion = negLoglikelihood,
+            hessian = TRUE)
+re
+optimReturn(re)
+
+MCEstimator(x = x, ParamFamily = nF, criterion = negLoglikelihood2,
+            fups="fu")
+re2 <- MCEstimator(x = x, ParamFamily = nF, criterion = negLoglikelihood2,
+            fups="fu", hessian = TRUE, fn="LU")
+re2
+optimReturn(re2)
+MCEstimator(x = x, ParamFamily = nF, criterion = negLoglikelihood3)
+MCEstimator(x = x, ParamFamily = nF, criterion = negLoglikelihood3, fn="LU")
+}
+##########################################################################
+# end added 2018 07 30
+##########################################################################
+
+
 ################################################################################
 ####### particular methods
 ################################################################################

Modified: branches/distr-2.8/pkg/distrMod/man/MCEstimate-class.Rd
===================================================================
--- branches/distr-2.8/pkg/distrMod/man/MCEstimate-class.Rd	2018-07-30 14:47:31 UTC (rev 1231)
+++ branches/distr-2.8/pkg/distrMod/man/MCEstimate-class.Rd	2018-07-30 17:47:38 UTC (rev 1232)
@@ -10,6 +10,8 @@
 \alias{method,MCEstimate-method}
 \alias{optimwarn}
 \alias{optimwarn,MCEstimate-method}
+\alias{optimReturn}
+\alias{optimReturn,MCEstimate-method}
 \alias{criterion<-}
 \alias{criterion<-,MCEstimate-method}
 \alias{coerce,MCEstimate,mle-method}
@@ -49,6 +51,9 @@
       additional informations. }
     \item{\code{optimwarn}}{ object of class \code{"character"}
       warnings issued during optimization. }
+    \item{\code{optimReturn}}{ object of class \code{"ANY"}
+      the return value of the optimizer (or \code{NULL} if, e.g.,
+      closed form solutions are used). }
     \item{\code{startPar}}{ --- object of class \code{"ANY"}; filled either
      with \code{NULL} (no starting value used) or with \code{"numeric"}
       --- the value of the starting parameter. }
@@ -88,6 +93,9 @@
     \item{optimwarn}{\code{signature(object = "MCEstimate")}:
       accessor function for slot \code{optimwarn}. }
 
+    \item{optimReturn}{\code{signature(object = "MCEstimate")}:
+      accessor function for slot \code{optimReturn}. }
+
     \item{startPar}{\code{signature(object = "MCEstimate")}:
       accessor function for slot \code{startPar}. }
 

Modified: branches/distr-2.8/pkg/distrMod/man/meRes.Rd
===================================================================
--- branches/distr-2.8/pkg/distrMod/man/meRes.Rd	2018-07-30 14:47:31 UTC (rev 1231)
+++ branches/distr-2.8/pkg/distrMod/man/meRes.Rd	2018-07-30 17:47:38 UTC (rev 1232)
@@ -9,7 +9,8 @@
 
 \usage{
 meRes(x, estimate, criterion.value, param, crit.fct, method = "explicit solution",
-      crit.name = "Maximum Likelihood", Infos, warns = "", startPar = NULL)
+      crit.name = "Maximum Likelihood", Infos, warns = "", startPar = NULL,
+      optReturn = NULL)
 get.criterion.fct(theta, Data, ParamFam, criterion.ff, fun, ...)
 \S4method{samplesize}{numeric}(object)
 }
@@ -39,7 +40,10 @@
    (with certain checking whether parameter value is permitted and possibly
     penalizing if not; see code to , for example.)}
   \item{startPar}{value of argument \code{StartPar} --- starting parameter used.}
-  \item{\dots}{further arguments to be passed 
+  \item{optReturn}{ object of class \code{"ANY"}
+      the return value of the optimizer (or \code{NULL} if, e.g.,
+      closed form solutions are used). }
+  \item{\dots}{further arguments to be passed
                to \code{optim} / \code{optimize}} 
   \item{object}{numeric; the data at which to evaluate the estimate}
   }



More information about the Distr-commits mailing list