add a new argument

lbbe-software · Jun 26, 2024 · 398a0f9 · 398a0f9
1 parent 3d0b3ce
commit 398a0f9
Show file tree

Hide file tree

Showing 25 changed files with 821 additions and 618 deletions.
diff --git a/NEWS.md b/NEWS.md
@@ -9,6 +9,7 @@ NEW FEATURES
 - make `gofstat()` work with `fitdistcens` objects (giving only AIC and BIC values).
 - add calculation of the hessian using `optimHess` within `fitdist` when it is not given by `optim`.
 - compute the asymptotic covariance matrix with MME : Now the theoretical moments `m<dist>` should be defined up to an order which equals to twice the maximal order given `order`.
+- add a new argument `calcvcov` in order to (dis)able the computation of covariance matrix for any method.
 - graphics function `*comp()` now return a list of drawn points and/or lines when `plotstyle == "graphics"`.
 - add a `density` function for `bootdist(cens)` objects.
 - add DOIs in man pages.

diff --git a/R/fitdist.R b/R/fitdist.R
@@ -22,8 +22,8 @@
 ###         R functions
 ### 
 
-fitdist <- function (data, distr, method = c("mle", "mme", "qme", "mge", "mse"), start=NULL, 
-                     fix.arg=NULL, discrete, keepdata = TRUE, keepdata.nb=100, ...) 
+fitdist <- function (data, distr, method = c("mle", "mme", "qme", "mge", "mse"), start = NULL, 
+                     fix.arg = NULL, discrete, keepdata = TRUE, keepdata.nb = 100, calcvcov = TRUE, ...) 
 {
     #check argument distr
     if (!is.character(distr)) 
@@ -49,6 +49,8 @@ fitdist <- function (data, distr, method = c("mle", "mme", "qme", "mge", "mse"),
       stop("wrong argument 'discrete'.")
     if(!is.logical(keepdata) || !is.numeric(keepdata.nb) || keepdata.nb < 2)
       stop("wrong arguments 'keepdata' and 'keepdata.nb'")
+    if(!is.logical(calcvcov) || length(calcvcov) != 1)
+      stop("wrong argument 'calcvcov'.")
     #check argument method
     if(any(method == "mom"))
         warning("the name \"mom\" for matching moments is NO MORE used and is replaced by \"mme\"")
@@ -104,11 +106,16 @@ fitdist <- function (data, distr, method = c("mle", "mme", "qme", "mge", "mse"),
                 stop("moment matching estimation needs an 'order' argument")   
 
         mme <- mmedist(data, distname, start=arg_startfix$start.arg, 
-                       fix.arg=arg_startfix$fix.arg, checkstartfix=TRUE, ...)
+                       fix.arg=arg_startfix$fix.arg, checkstartfix=TRUE, 
+                       calcvcov=calcvcov, ...)
 
         varcovar <- mme$vcov
-        correl <- cov2cor(varcovar)
-        sd <- sqrt(diag(varcovar))
+        if(!is.null(varcovar))
+        {
+          correl <- cov2cor(varcovar)
+          sd <- sqrt(diag(varcovar))
+        }else
+          correl <- sd <- NULL
 
         estimate <- mme$estimate
         loglik <- mme$loglik
@@ -121,29 +128,20 @@ fitdist <- function (data, distr, method = c("mle", "mme", "qme", "mge", "mse"),
     }else if (method == "mle")
     {
         mle <- mledist(data, distname, start=arg_startfix$start.arg, 
-                       fix.arg=arg_startfix$fix.arg, checkstartfix=TRUE, ...)
+                       fix.arg=arg_startfix$fix.arg, checkstartfix=TRUE, 
+                       calcvcov=calcvcov, ...)
         if (mle$convergence>0) 
            stop("the function mle failed to estimate the parameters, 
                 with the error code ", mle$convergence, "\n") 
         estimate <- mle$estimate
-        if(!is.null(mle$hessian)){
-            #check for NA values and invertible Hessian
-            if(all(!is.na(mle$hessian)) && qr(mle$hessian)$rank == NCOL(mle$hessian))
-            {
-                varcovar <- solve(mle$hessian)
-                sd <- sqrt(diag(varcovar))
-                correl <- cov2cor(varcovar)
-            }else
-            {
-                varcovar <- NA
-                sd <- NA
-                correl <- NA                            
-            }
-        }else{
-            varcovar <- NA
-            sd <- NA
-            correl <- NA            
-        }
+        varcovar <- mle$vcov
+        if(!is.null(varcovar))
+        {
+          correl <- cov2cor(varcovar)
+          sd <- sqrt(diag(varcovar))
+        }else
+          correl <- sd <- NULL
+
         loglik <- mle$loglik
         npar <- length(estimate)
         aic <- -2*loglik+2*npar
@@ -157,15 +155,15 @@ fitdist <- function (data, distr, method = c("mle", "mme", "qme", "mge", "mse"),
             stop("quantile matching estimation needs an 'probs' argument") 
 
         qme <- qmedist(data, distname, start=arg_startfix$start.arg, 
-                       fix.arg=arg_startfix$fix.arg, checkstartfix=TRUE, ...)
+                       fix.arg=arg_startfix$fix.arg, checkstartfix=TRUE, 
+                       calcvcov=calcvcov, ...)
 
         estimate <- qme$estimate
-        sd <- NULL
         loglik <- qme$loglik
         npar <- length(estimate)
         aic <- -2*loglik+2*npar
         bic <- -2*loglik+log(n)*npar
-        correl <- varcovar <- NULL
+        sd <- correl <- varcovar <- NULL
 
         convergence <- qme$convergence   
         fix.arg <- qme$fix.arg
@@ -176,38 +174,38 @@ fitdist <- function (data, distr, method = c("mle", "mme", "qme", "mge", "mse"),
             warning("maximum GOF estimation has a default 'gof' argument set to 'CvM'")    
 
         mge <- mgedist(data, distname, start=arg_startfix$start.arg, 
-                       fix.arg=arg_startfix$fix.arg, checkstartfix=TRUE, ...)
+                       fix.arg=arg_startfix$fix.arg, checkstartfix=TRUE, 
+                       calcvcov=calcvcov, ...)
 
         estimate <- mge$estimate
-        sd <- NULL
         loglik <- mge$loglik
         npar <- length(estimate)
         aic <- -2*loglik+2*npar
         bic <- -2*loglik+log(n)*npar
-        correl <- varcovar <- NULL
+        sd <- correl <- varcovar <- NULL
 
         convergence <- mge$convergence
         fix.arg <- mge$fix.arg
         weights <- NULL
     }else if (method == "mse")
     {
       mse <- msedist(data, distname, start=arg_startfix$start.arg, 
-                     fix.arg=arg_startfix$fix.arg, checkstartfix=TRUE, ...)
+                     fix.arg=arg_startfix$fix.arg, checkstartfix=TRUE, 
+                     calcvcov=calcvcov, ...)
 
       estimate <- mse$estimate
-      sd <- NULL
       loglik <- mse$loglik
       npar <- length(estimate)
       aic <- -2*loglik+2*npar
       bic <- -2*loglik+log(n)*npar
-      correl <- varcovar <- NULL
+      sd <- correl <- varcovar <- NULL
 
       convergence <- mse$convergence
       fix.arg <- mse$fix.arg
       weights <- mse$weights
     }else
     {
-        stop("match.arg() did not work correctly")
+        stop("match.arg() for 'method' did not work correctly")
     }
 
     #needed for bootstrap

diff --git a/R/mgedist.R b/R/mgedist.R
@@ -26,7 +26,7 @@
 
 mgedist <- function (data, distr, gof = "CvM", start=NULL, fix.arg=NULL, optim.method="default",
     lower=-Inf, upper=Inf, custom.optim=NULL, silent=TRUE, gradient=NULL, 
-    checkstartfix=FALSE, ...)
+    checkstartfix=FALSE, calcvcov=FALSE, ...)
     # data may correspond to a vector for non censored data or to
     # a dataframe of two columns named left and right for censored data 
 {
@@ -368,6 +368,5 @@ mgedist <- function (data, distr, gof = "CvM", start=NULL, fix.arg=NULL, optim.m
                     counts=opt$counts, optim.message=opt$message,
                     loglik=loglik(opt$par, fix.arg, data, ddistname), gof=gof)
     }   
-
     return(res)                
 }
diff --git a/R/mledist.R b/R/mledist.R
@@ -26,7 +26,7 @@
 
 mledist <- function (data, distr, start=NULL, fix.arg=NULL, optim.method="default", 
     lower=-Inf, upper=Inf, custom.optim=NULL, weights=NULL, silent=TRUE, gradient=NULL, 
-    checkstartfix=FALSE, ...)
+    checkstartfix=FALSE, calcvcov=FALSE, ...)
     # data may correspond to a vector for non censored data or to
     # a dataframe of two columns named left and right for censored data 
 {
@@ -298,7 +298,7 @@ mledist <- function (data, distr, start=NULL, fix.arg=NULL, optim.method="defaul
             #recompute Hessian if Hessian matrix has NA values
             if(!inherits(opttryerror, "try-error"))
             {
-              if(any(is.na(opt$hessian)) || is.null(opt$hessian)) 
+              if(calcvcov && any(is.na(opt$hessian)) || is.null(opt$hessian)) 
               {
                 opthessian <- try(optimHess(par=vstart, fn=fnobj, fix.arg=fix.arg, obs=data, gr=gradient, ddistnam=ddistname), silent=TRUE)
                 if(!inherits(opthessian, "try-error"))
@@ -320,7 +320,7 @@ mledist <- function (data, distr, start=NULL, fix.arg=NULL, optim.method="defaul
             #recompute Hessian if Hessian matrix has NA values
             if(!inherits(opttryerror, "try-error"))
             {
-              if(any(is.na(opt$hessian)) || is.null(opt$hessian)) 
+              if(calcvcov && any(is.na(opt$hessian)) || is.null(opt$hessian)) 
               {
                 opthessian <- try(optimHess(par=vstart, fn=fnobj, fix.arg=fix.arg, obs=data, gr=gradient, ddistnam=ddistname), 
                                   silent=TRUE)
@@ -346,7 +346,7 @@ mledist <- function (data, distr, start=NULL, fix.arg=NULL, optim.method="defaul
           #recompute Hessian if Hessian matrix has NA values
           if(!inherits(opttryerror, "try-error"))
           {
-            if(any(is.na(opt$hessian)) || is.null(opt$hessian)) 
+            if(calcvcov && any(is.na(opt$hessian)) || is.null(opt$hessian)) 
             {
               opthessian <- try(optimHess(par=vstart, fn=fnobj, fix.arg=fix.arg, obs=data, gr=gradient, ddistnam=ddistname),
                                 silent=TRUE)
@@ -362,8 +362,9 @@ mledist <- function (data, distr, start=NULL, fix.arg=NULL, optim.method="defaul
             warnings("The function optim encountered an error and stopped.")
             if(getOption("show.error.messages")) print(attr(opttryerror, "condition"))          
             return(list(estimate = rep(NA, length(vstart)), convergence = 100, loglik = NA, 
-                        hessian = NA, optim.function=opt.fun, fix.arg = fix.arg, 
-                        optim.method=meth, fix.arg.fun = fix.arg.fun, counts=c(NA, NA)))
+                        hessian = NA, optim.function = opt.fun, fix.arg = fix.arg, 
+                        optim.method=meth, fix.arg.fun = fix.arg.fun, counts=c(NA, NA), 
+                        vcov = NULL))
         }
 
         if (opt$convergence>0) {
@@ -372,10 +373,20 @@ mledist <- function (data, distr, start=NULL, fix.arg=NULL, optim.method="defaul
         }
         if(is.null(names(opt$par)))
           names(opt$par) <- names(vstart)
+        if(calcvcov && !is.null(opt$hessian))
+        {
+          #see R/util-mledist-vcov.R
+          varcovar <- mle.vcov(opt$hessian)
+          #add names
+          if(!is.null(varcovar))
+            colnames(varcovar) <- rownames(varcovar) <- names(opt$par)
+        }else
+          varcovar <- NULL
         res <- list(estimate = opt$par, convergence = opt$convergence, value=opt$value,  
                     hessian = opt$hessian, optim.function=opt.fun, optim.method=meth, 
                     fix.arg = fix.arg, fix.arg.fun = fix.arg.fun, weights = weights, 
-                    counts=opt$counts, optim.message=opt$message, loglik = -opt$value)
+                    counts=opt$counts, optim.message=opt$message, loglik = -opt$value,
+                    vcov=varcovar)
     }
     else # Try to minimize the minus (log-)likelihood using a user-supplied optim function 
     {
@@ -384,7 +395,7 @@ mledist <- function (data, distr, start=NULL, fix.arg=NULL, optim.method="defaul
             opttryerror <- try(opt <- custom.optim(fn=fnobj, fix.arg=fix.arg, obs=data, 
                 ddistnam=ddistname, par=vstart, ...), silent=TRUE)
         else
-            opttryerror <-try(opt<-custom.optim(fn=fnobjcens, fix.arg=fix.arg, rcens=rcens, 
+            opttryerror <-try(opt <- custom.optim(fn=fnobjcens, fix.arg=fix.arg, rcens=rcens, 
                 lcens=lcens, icens=icens, ncens=ncens, ddistnam=ddistname, pdistnam=pdistname, 
                 par=vstart, ...), silent=TRUE)              
         options(warn=owarn)
@@ -404,12 +415,21 @@ mledist <- function (data, distr, start=NULL, fix.arg=NULL, optim.method="defaul
         }
         if(is.null(names(opt$par)))
           names(opt$par) <- names(vstart)
+        if(calcvcov && !is.null(opt$hessian))
+        {
+          varcovar <- mle.vcov(opt$hessian)
+          #add names
+          if(!is.null(varcovar))
+            colnames(varcovar) <- rownames(varcovar) <- names(opt$par)
+        }else
+          varcovar <- NULL
         argdot <- list(...)
         method.cust <- argdot$method
         res <- list(estimate = opt$par, convergence = opt$convergence, value=opt$value, 
                     hessian = opt$hessian, optim.function = custom.optim, optim.method = method.cust, 
                     fix.arg = fix.arg, fix.arg.fun = fix.arg.fun, weights = weights, 
-                    counts=opt$counts, optim.message=opt$message, loglik = -opt$value)        
+                    counts=opt$counts, optim.message=opt$message, loglik = -opt$value,
+                    vcov=varcovar)
     }   
 
     return(res) 
@@ -418,5 +438,5 @@ mledist <- function (data, distr, start=NULL, fix.arg=NULL, optim.method="defaul
 ## old function with previous name for censored data
 mledistcens <- function(censdata, distr, start=NULL, optim.method="default", lower=-Inf, upper=Inf)
 {
-    stop("The function \"mledistcens\" is no more used. Now the same function \"mledist\" must be used both for censored and non censored data.")
+    stop("The function \"mledistcens\" is deprecated. Now the same function \"mledist\" must be used both for censored and non censored data.")
 }