From 3e157fddf1d00db8449b1163ba3c1eaf1bf03b13 Mon Sep 17 00:00:00 2001
From: Christophe Dutang <dutangc@gmail.com>
Date: Mon, 10 Jun 2024 11:38:43 +0200
Subject: [PATCH] add defensive programming

---
 R/fitdist.R              |   1 +
 R/fitdistcens.R          | 343 ++++++++++++++++++++-------------------
 R/mgedist.R              |   1 +
 R/mledist.R              |   1 +
 R/mmedist.R              |   1 +
 R/msedist.R              |   1 +
 R/qmedist.R              |   1 +
 R/util-check-NAInfNan.R  |  36 ++++
 share/todolist.md        |   4 +-
 tests/t-fitdist.R        |  11 ++
 tests/t-fitdistcens.R    |  19 +++
 tests/t-mgedist.R        |   7 +
 tests/t-mledist-nocens.R |  10 +-
 tests/t-mmedist.R        |   7 +
 tests/t-msedist.R        |  12 ++
 tests/t-qmedist.R        |  10 ++
 16 files changed, 291 insertions(+), 174 deletions(-)
 create mode 100644 R/util-check-NAInfNan.R

diff --git a/R/fitdist.R b/R/fitdist.R
index 6f4484dd..a213ce4b 100644
--- a/R/fitdist.R
+++ b/R/fitdist.R
@@ -61,6 +61,7 @@ fitdist <- function (data, distr, method = c("mle", "mme", "qme", "mge", "mse"),
     #check argument data
     if (!(is.vector(data) & is.numeric(data) & length(data)>1))
         stop("data must be a numeric vector of length greater than 1")
+    checkUncensoredNAInfNan(data)
  
     #encapsulate three dots arguments
     my3dots <- list(...)    
diff --git a/R/fitdistcens.R b/R/fitdistcens.R
index 906e01c5..57057be6 100644
--- a/R/fitdistcens.R
+++ b/R/fitdistcens.R
@@ -25,198 +25,199 @@
 fitdistcens <- function (censdata, distr, start=NULL, fix.arg=NULL, 
                          keepdata = TRUE, keepdata.nb=100, ...) 
 {
-    if (missing(censdata) ||
-        !(is.vector(censdata$left) & is.vector(censdata$right) & length(censdata[, 1])>1))
-        stop("datacens must be a dataframe with two columns named left 
+  if (missing(censdata) || !is.data.frame(censdata) ||
+      !(is.vector(censdata$left) & is.vector(censdata$right) & length(censdata[, 1])>1))
+    stop("censdata must be a dataframe with two columns named left 
             and right and more than one line")
-    leftsupright <- censdata$left>censdata$right
-    leftsupright <- leftsupright[!is.na(leftsupright)]
-    if (any(leftsupright))
-        stop("each left value must be less or equal to the corresponding right value")
-    if (!is.character(distr)) distname <- substring(as.character(match.call()$distr), 2)
-    else distname <- distr
-    ddistname <- paste("d", distname, sep="")
-    if (!exists(ddistname, mode="function"))
-        stop(paste("The ", ddistname, " function must be defined"))
-    pdistname <- paste("p", distname, sep="")
-    if (!exists(pdistname, mode="function"))
-        stop(paste("The ", pdistname, " function must be defined"))
-    if(!is.logical(keepdata) || !is.numeric(keepdata.nb) || keepdata.nb < 3)
-      stop("wrong arguments 'keepdata' and 'keepdata.nb'.")
-    
-    #encapsulate three dots arguments
-    my3dots <- list(...)    
-    if (length(my3dots) == 0) 
-      my3dots <- NULL
-    
-    #format data for calculation of starting values
-    pseudodata <- cens2pseudo(censdata)$pseudo
-    
-    # manage starting/fixed values: may raise errors or return two named list
-    arg_startfix <- manageparam(start.arg=start, fix.arg=fix.arg, obs=pseudodata, 
-                                distname=distname)
-    
-    #check inconsistent parameters
-    argddistname <- names(formals(ddistname))
-    hasnodefaultval <- sapply(formals(ddistname), is.name)
-    arg_startfix <- checkparamlist(arg_startfix$start.arg, arg_startfix$fix.arg, 
-                                   argddistname, hasnodefaultval)
-    #arg_startfix contains two names list (no longer NULL nor function)
-    
-    #store fix.arg.fun if supplied by the user
-    if(is.function(fix.arg))
-      fix.arg.fun <- fix.arg
-    else
-      fix.arg.fun <- NULL
-    
-    # check d, p, q, functions of distname
-    dpq2test <- c("d", "p")
-    resdpq <- testdpqfun(distname, dpq2test, start.arg=arg_startfix$start.arg, 
-                         fix.arg=arg_startfix$fix.arg, discrete=FALSE)
-    if(any(!resdpq$ok))
-    {
-      for(x in resdpq[!resdpq$ok, "txt"])
-        warning(x)
-    }
-    
-        
-    # MLE fit with mledist 
-    mle <- mledist(censdata, distname, start=arg_startfix$start.arg, 
-                   fix.arg=arg_startfix$fix.arg, checkstartfix=TRUE, ...)
-    if (mle$convergence>0) 
-        stop("the function mle failed to estimate the parameters, 
+  
+  checkCensoredDataFrameNAInfNan(censdata)
+  
+  if (!is.character(distr)) 
+    distname <- substring(as.character(match.call()$distr), 2)
+  else 
+    distname <- distr
+  ddistname <- paste("d", distname, sep="")
+  if (!exists(ddistname, mode="function"))
+    stop(paste("The ", ddistname, " function must be defined"))
+  pdistname <- paste("p", distname, sep="")
+  if (!exists(pdistname, mode="function"))
+    stop(paste("The ", pdistname, " function must be defined"))
+  if(!is.logical(keepdata) || !is.numeric(keepdata.nb) || keepdata.nb < 3)
+    stop("wrong arguments 'keepdata' and 'keepdata.nb'.")
+  
+  #encapsulate three dots arguments
+  my3dots <- list(...)    
+  if (length(my3dots) == 0) 
+    my3dots <- NULL
+  
+  #format data for calculation of starting values
+  pseudodata <- cens2pseudo(censdata)$pseudo
+  
+  # manage starting/fixed values: may raise errors or return two named list
+  arg_startfix <- manageparam(start.arg=start, fix.arg=fix.arg, obs=pseudodata, 
+                              distname=distname)
+  
+  #check inconsistent parameters
+  argddistname <- names(formals(ddistname))
+  hasnodefaultval <- sapply(formals(ddistname), is.name)
+  arg_startfix <- checkparamlist(arg_startfix$start.arg, arg_startfix$fix.arg, 
+                                 argddistname, hasnodefaultval)
+  #arg_startfix contains two names list (no longer NULL nor function)
+  
+  #store fix.arg.fun if supplied by the user
+  if(is.function(fix.arg))
+    fix.arg.fun <- fix.arg
+  else
+    fix.arg.fun <- NULL
+  
+  # check d, p, q, functions of distname
+  dpq2test <- c("d", "p")
+  resdpq <- testdpqfun(distname, dpq2test, start.arg=arg_startfix$start.arg, 
+                       fix.arg=arg_startfix$fix.arg, discrete=FALSE)
+  if(any(!resdpq$ok))
+  {
+    for(x in resdpq[!resdpq$ok, "txt"])
+      warning(x)
+  }
+  
+  
+  # MLE fit with mledist 
+  mle <- mledist(censdata, distname, start=arg_startfix$start.arg, 
+                 fix.arg=arg_startfix$fix.arg, checkstartfix=TRUE, ...)
+  if (mle$convergence>0) 
+    stop("the function mle failed to estimate the parameters, 
         with the error code ", mle$convergence) 
-    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
-        }
+  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
-    }
-    
-    loglik <- mle$loglik
-    n <- nrow(censdata)
-    npar <- length(estimate)
-    aic <- -2*loglik+2*npar
-    bic <- -2*loglik+log(n)*npar
-         
-    fix.arg <- mle$fix.arg
-    weights <- mle$weights
-    
-    #needed for bootstrap
-    if (!is.null(fix.arg)) 
-      fix.arg <- as.list(fix.arg)
-    
-    if(keepdata)
-    {
-      reslist <- list(estimate = estimate, method="mle", sd = sd, cor = correl, 
-                      vcov = varcovar, loglik = loglik, aic=aic, bic=bic, n=n, censdata=censdata, 
-                      distname=distname, fix.arg=fix.arg, fix.arg.fun = fix.arg.fun, 
-                      dots=my3dots, convergence=mle$convergence, discrete=FALSE,
-                      weights = weights)
-    }else
-    {
-      n2keep <- min(keepdata.nb, n)-4
-      imin <- unique(apply(censdata, 2, which.min))
-      imax <- unique(apply(censdata, 2, which.max))
-      subdata <- censdata[sample((1:n)[-c(imin, imax)], size=n2keep, replace=FALSE), ]
-      subdata <- rbind.data.frame(subdata, censdata[c(imin, imax), ])
-      
-      reslist <- list(estimate = estimate, method="mle", sd = sd, cor = correl, 
-                      vcov = varcovar, loglik = loglik, aic=aic, bic=bic, n=n, censdata=subdata, 
-                      distname=distname, fix.arg=fix.arg, fix.arg.fun = fix.arg.fun, 
-                      dots=my3dots, convergence=mle$convergence, discrete=FALSE,
-                      weights = weights)
+      varcovar <- NA
+      sd <- NA
+      correl <- NA
     }
-    
-    return(structure(reslist, class = "fitdistcens"))
-        
+  }else{
+    varcovar <- NA
+    sd <- NA
+    correl <- NA
+  }
+  
+  loglik <- mle$loglik
+  n <- nrow(censdata)
+  npar <- length(estimate)
+  aic <- -2*loglik+2*npar
+  bic <- -2*loglik+log(n)*npar
+  
+  fix.arg <- mle$fix.arg
+  weights <- mle$weights
+  
+  #needed for bootstrap
+  if (!is.null(fix.arg)) 
+    fix.arg <- as.list(fix.arg)
+  
+  if(keepdata)
+  {
+    reslist <- list(estimate = estimate, method="mle", sd = sd, cor = correl, 
+                    vcov = varcovar, loglik = loglik, aic=aic, bic=bic, n=n, censdata=censdata, 
+                    distname=distname, fix.arg=fix.arg, fix.arg.fun = fix.arg.fun, 
+                    dots=my3dots, convergence=mle$convergence, discrete=FALSE,
+                    weights = weights)
+  }else
+  {
+    n2keep <- min(keepdata.nb, n)-4
+    imin <- unique(apply(censdata, 2, which.min))
+    imax <- unique(apply(censdata, 2, which.max))
+    subdata <- censdata[sample((1:n)[-c(imin, imax)], size=n2keep, replace=FALSE), ]
+    subdata <- rbind.data.frame(subdata, censdata[c(imin, imax), ])
+    
+    reslist <- list(estimate = estimate, method="mle", sd = sd, cor = correl, 
+                    vcov = varcovar, loglik = loglik, aic=aic, bic=bic, n=n, censdata=subdata, 
+                    distname=distname, fix.arg=fix.arg, fix.arg.fun = fix.arg.fun, 
+                    dots=my3dots, convergence=mle$convergence, discrete=FALSE,
+                    weights = weights)
+  }
+  
+  return(structure(reslist, class = "fitdistcens"))
+  
 }
 
 print.fitdistcens <- function(x, ...){
-    if (!inherits(x, "fitdistcens"))
-        stop("Use only with 'fitdistcens' objects")
-    cat("Fitting of the distribution '", x$distname, "' on censored data by maximum likelihood \n")
-    cat("Parameters:\n")
-    print(cbind.data.frame("estimate" = x$estimate), ...)
-
-    if(!is.null(x$fix.arg))
+  if (!inherits(x, "fitdistcens"))
+    stop("Use only with 'fitdistcens' objects")
+  cat("Fitting of the distribution '", x$distname, "' on censored data by maximum likelihood \n")
+  cat("Parameters:\n")
+  print(cbind.data.frame("estimate" = x$estimate), ...)
+  
+  if(!is.null(x$fix.arg))
+  {
+    if(is.null(x$fix.arg.fun))
     {
-      if(is.null(x$fix.arg.fun))
-      {
-        cat("Fixed parameters:\n")
-      }else
-      {
-        cat("Fixed parameters (computed by a user-supplied function):\n")
-      }
-      print(cbind.data.frame("value" = unlist(x$fix.arg)), ...)
+      cat("Fixed parameters:\n")
+    }else
+    {
+      cat("Fixed parameters (computed by a user-supplied function):\n")
     }
-
+    print(cbind.data.frame("value" = unlist(x$fix.arg)), ...)
+  }
+  
 }
 
 plot.fitdistcens <- function(x, ...){
-    if (!inherits(x, "fitdistcens"))
-        stop("Use only with 'fitdistcens' objects")
-    if(!is.null(x$weights))
-      stop("The plot of the fit is not yet available when using weights")
+  if (!inherits(x, "fitdistcens"))
+    stop("Use only with 'fitdistcens' objects")
+  if(!is.null(x$weights))
+    stop("The plot of the fit is not yet available when using weights")
   plotdistcens(censdata=x$censdata, distr=x$distname, 
-    para=c(as.list(x$estimate), as.list(x$fix.arg)), ...)
-    
+               para=c(as.list(x$estimate), as.list(x$fix.arg)), ...)
+  
 }
 
 summary.fitdistcens <- function(object, ...){
-    if (!inherits(object, "fitdistcens"))
-        stop("Use only with 'fitdistcens' objects")
-    object$ddistname <- paste("d", object$distname, sep="")
-    object$pdistname <- paste("p", object$distname, sep="")
-    
-    class(object) <- c("summary.fitdistcens", class(object))    
-    object
+  if (!inherits(object, "fitdistcens"))
+    stop("Use only with 'fitdistcens' objects")
+  object$ddistname <- paste("d", object$distname, sep="")
+  object$pdistname <- paste("p", object$distname, sep="")
+  
+  class(object) <- c("summary.fitdistcens", class(object))    
+  object
 }
 
 print.summary.fitdistcens <- function(x, ...){
-    if (!inherits(x, "summary.fitdistcens"))
-        stop("Use only with 'fitdistcens' objects")
-    ddistname <- x$ddistname
-    pdistname <- x$pdistname
-    
-    cat("Fitting of the distribution '", x$distname, 
-    "' By maximum likelihood on censored data \n")
-    cat("Parameters\n")
-    print(cbind.data.frame("estimate" = x$estimate, "Std. Error" = x$sd), ...)
-    
-    if(!is.null(x$fix.arg))
+  if (!inherits(x, "summary.fitdistcens"))
+    stop("Use only with 'fitdistcens' objects")
+  ddistname <- x$ddistname
+  pdistname <- x$pdistname
+  
+  cat("Fitting of the distribution '", x$distname, 
+      "' By maximum likelihood on censored data \n")
+  cat("Parameters\n")
+  print(cbind.data.frame("estimate" = x$estimate, "Std. Error" = x$sd), ...)
+  
+  if(!is.null(x$fix.arg))
+  {
+    if(is.null(x$fix.arg.fun))
     {
-      if(is.null(x$fix.arg.fun))
-      {
-        cat("Fixed parameters:\n")
-      }else
-      {
-        cat("Fixed parameters (computed by a user-supplied function):\n")
-      }
-      print(cbind.data.frame("value" = unlist(x$fix.arg)), ...)
-    }
-    
-    cat("Loglikelihood: ", x$loglik, "  ")
-    cat("AIC: ", x$aic, "  ")
-    cat("BIC: ", x$bic, "\n")
-    if (length(x$estimate) > 1) {
-        cat("Correlation matrix:\n")
-        print(x$cor)
-        cat("\n")
+      cat("Fixed parameters:\n")
+    }else
+    {
+      cat("Fixed parameters (computed by a user-supplied function):\n")
     }
-    invisible(x)
+    print(cbind.data.frame("value" = unlist(x$fix.arg)), ...)
+  }
+  
+  cat("Loglikelihood: ", x$loglik, "  ")
+  cat("AIC: ", x$aic, "  ")
+  cat("BIC: ", x$bic, "\n")
+  if (length(x$estimate) > 1) {
+    cat("Correlation matrix:\n")
+    print(x$cor)
+    cat("\n")
+  }
+  invisible(x)
 }
 
 #see quantiles.R for quantile.fitdistcens
diff --git a/R/mgedist.R b/R/mgedist.R
index ce081fe9..00ca0c69 100644
--- a/R/mgedist.R
+++ b/R/mgedist.R
@@ -65,6 +65,7 @@ mgedist <- function (data, distr, gof = "CvM", start=NULL, fix.arg=NULL, optim.m
         if (!(is.numeric(data) & length(data)>1)) 
             stop("data must be a numeric vector of length greater than 1 for non censored data
             or a dataframe with two columns named left and right and more than one line for censored data")
+        checkUncensoredNAInfNan(data)
     }
     else {
         cens <- TRUE
diff --git a/R/mledist.R b/R/mledist.R
index ce856b2a..08efe53d 100644
--- a/R/mledist.R
+++ b/R/mledist.R
@@ -63,6 +63,7 @@ mledist <- function (data, distr, start=NULL, fix.arg=NULL, optim.method="defaul
         cens <- FALSE
         if (!(is.numeric(data) & length(data)>1)) 
             stop(paste(txt1, txt2))
+        checkUncensoredNAInfNan(data)
     }
     else {
         cens <- TRUE
diff --git a/R/mmedist.R b/R/mmedist.R
index 21436b97..334840ab 100644
--- a/R/mmedist.R
+++ b/R/mmedist.R
@@ -63,6 +63,7 @@ mmedist <- function (data, distr, order, memp, start=NULL, fix.arg=NULL,
     }
     if (!(is.numeric(data) & length(data)>1)) 
         stop("data must be a numeric vector of length greater than 1")
+    checkUncensoredNAInfNan(data)
     
     if(is.null(weights))
     {  
diff --git a/R/msedist.R b/R/msedist.R
index 4050a4cb..1bc54c55 100644
--- a/R/msedist.R
+++ b/R/msedist.R
@@ -83,6 +83,7 @@ msedist <- function (data, distr, phidiv="KL", power.phidiv=NULL, start=NULL, fi
         if (!(is.numeric(data) & length(data)>1)) 
             stop("data must be a numeric vector of length greater than 1 for non censored data
             or a dataframe with two columns named left and right and more than one line for censored data")
+        checkUncensoredNAInfNan(data)
     }
     else {
         cens <- TRUE
diff --git a/R/qmedist.R b/R/qmedist.R
index 65ad776f..1cd452ca 100644
--- a/R/qmedist.R
+++ b/R/qmedist.R
@@ -69,6 +69,7 @@ qmedist <- function (data, distr, probs, start=NULL, fix.arg=NULL,
         if (!(is.numeric(data) & length(data)>1)) 
             stop("data must be a numeric vector of length greater than 1 for non censored data
             or a dataframe with two columns named left and right and more than one line for censored data")
+        checkUncensoredNAInfNan(data)
     } else 
     {
         cens <- TRUE
diff --git a/R/util-check-NAInfNan.R b/R/util-check-NAInfNan.R
new file mode 100644
index 00000000..2c71c115
--- /dev/null
+++ b/R/util-check-NAInfNan.R
@@ -0,0 +1,36 @@
+
+
+checkUncensoredNAInfNan <- function(x)
+{
+  if(!is.numeric(x))
+    stop("data is not a numeric vector.")
+  if(any(is.nan(x)))
+    stop("data contain NaN (not a numeric) values.")
+  if(any(is.infinite(x)))
+    stop("data contain Inf (infinite) values.")
+  if(any(is.na(x)))
+    stop("data contain NA values.")
+  invisible(NULL)
+}
+
+
+checkCensoredDataFrameNAInfNan <- function(x)
+{
+  if(!is.data.frame(x))
+    stop("censdata is not a dataframe with two columns named left and right.")
+  if(NCOL(x) != 2)
+    stop("censdata is not a dataframe with two columns named left and right.")
+  if(!"left" %in% colnames(x) || !"right" %in% colnames(x))
+    stop("censdata is not a dataframe with two columns named left and right.")
+  if(any(is.nan(x$left) | is.nan(x$right)))
+    stop("censdata contain NaN (not a numeric) values.")
+  if(any(is.infinite(x$left) | is.infinite(x$right)))
+    stop("censdata contain Inf (infinite) values.")
+  if(any(is.na(x$left) & is.na(x$right)))
+    stop("censdata contain two NA values on the same line.")
+  leftsupright <- x$left > x$right
+  leftsupright <- leftsupright[!is.na(leftsupright)]
+  if (any(leftsupright))
+    stop("each censdata$left value must be less or equal to the corresponding censdata$right value")
+  invisible(NULL)
+}
diff --git a/share/todolist.md b/share/todolist.md
index 8c99fcb0..c07d19e9 100644
--- a/share/todolist.md
+++ b/share/todolist.md
@@ -30,8 +30,8 @@
 26. [x] make introduction for vignette `starting values` and update layout, update README on `github.io` *CD*
 27. [x] update `mledist.Rd` with new starting value list (all distribution in `actuar` except phase-type), link with vignette *CD*
 28. [x] check that all base R distribution have starting values *CD*
-29. [ ] add defensive programming for infinite / NaN / NA values for `fitdist` objects *MLDM* et *CD*
-30. [ ] add defensive programming for NaN values for `fitdistcens` objects, convert `Inf` to `NA` (raise error for inconsistent values), raise error for double `NA` *MLDM* et *CD* 
+29. [x] add defensive programming for infinite / NaN / NA values for `fitdist` objects *MLDM* et *CD*
+30. [x] add defensive programming for NaN values for `fitdistcens` objects, convert `Inf` to `NA` (raise error for inconsistent values), raise error for double `NA` *MLDM* et *CD* 
 31. [x] add a link to `Surv2fitdistcens` in `fitdistcens`, add `Surv2fitdistcens` when infinite values
 32. [x] check link of `Surv2fitdistcens` in the FAQ
 33. [ ] add exploratory tools for assessing the role of covariates and for choosing the appropriate distribution in modelling (survival parametric model), e.g. dose-response
diff --git a/tests/t-fitdist.R b/tests/t-fitdist.R
index 8de8bb7e..cfdb91c4 100644
--- a/tests/t-fitdist.R
+++ b/tests/t-fitdist.R
@@ -26,6 +26,17 @@ if(visualize) {  # check ERROR on aarch64-apple-darwin20.4.0 (64-bit) (2021/05/1
   names(fitdist(serving, "gamma", optim.method="L-BFGS-B", lower=0)$estimate)
 }
 
+# (2) sanity check
+#
+
+
+try(fitdist(c(serving, "a"), "gamma"))
+try(fitdist(c(serving, NA), "gamma"))
+try(fitdist(c(serving, Inf), "gamma"))
+try(fitdist(c(serving, -Inf), "gamma"))
+try(fitdist(c(serving, NaN), "gamma"))
+
+
 # (7) custom optimization function
 #
 
diff --git a/tests/t-fitdistcens.R b/tests/t-fitdistcens.R
index d000a9b7..be632369 100644
--- a/tests/t-fitdistcens.R
+++ b/tests/t-fitdistcens.R
@@ -5,6 +5,25 @@ nsample <- 10
 visualize <- FALSE # TRUE for manual tests with visualization of results
 set.seed(1234)
 
+
+# (1) check input data
+#
+data(smokedfish)
+fitsf  <-  fitdistcens(smokedfish, "lnorm")
+summary(fitsf)
+
+
+smokedfish2 <- tail(smokedfish, 20)
+
+try(fitdistcens(rbind(smokedfish2, c(NA, NA)) , "lnorm"))
+
+try(fitdistcens(rbind(smokedfish2, c(3, Inf)) , "lnorm"))
+
+try(fitdistcens(rbind(smokedfish2, c(NaN, 3)) , "lnorm"))
+
+try(fitdistcens(rbind(smokedfish2, c(5, 3)) , "lnorm"))
+
+
 # (6) custom optimisation function - example with the genetic algorithm
 #
 data(fluazinam)
diff --git a/tests/t-mgedist.R b/tests/t-mgedist.R
index e9ecd95c..99f15bd6 100644
--- a/tests/t-mgedist.R
+++ b/tests/t-mgedist.R
@@ -92,6 +92,13 @@ res <- mgedist(x, "norm2", start=list(a=1))
 attr(try(log("a"), silent=TRUE), "condition")
 
 
+try(mgedist(c(serving, "a"), "gamma"))
+try(mgedist(c(serving, NA), "gamma"))
+try(mgedist(c(serving, Inf), "gamma"))
+try(mgedist(c(serving, -Inf), "gamma"))
+try(mgedist(c(serving, NaN), "gamma"))
+
+
 # (7) test the component optim.message
 x <- rnorm(1000)
 #change parameter to obtain unsuccessful convergence
diff --git a/tests/t-mledist-nocens.R b/tests/t-mledist-nocens.R
index 9653d3fd..d7fd642c 100644
--- a/tests/t-mledist-nocens.R
+++ b/tests/t-mledist-nocens.R
@@ -1,7 +1,8 @@
 library(fitdistrplus)
 
 
-
+data("groundbeef")
+serving <- groundbeef$serving 
 
 # (1) basic fit of a normal distribution with maximum likelihood estimation
 #
@@ -257,6 +258,13 @@ res <- mledist(x, "norm2", start=list(a=1))
 attr(try(log("a"), silent=TRUE), "condition")
 
 
+try(mledist(c(serving, "a"), "gamma"))
+try(mledist(c(serving, NA), "gamma"))
+try(mledist(c(serving, Inf), "gamma"))
+try(mledist(c(serving, -Inf), "gamma"))
+try(mledist(c(serving, NaN), "gamma"))
+
+
 
 # (13) weighted MLE
 #
diff --git a/tests/t-mmedist.R b/tests/t-mmedist.R
index 744ecc51..3aa13ad7 100644
--- a/tests/t-mmedist.R
+++ b/tests/t-mmedist.R
@@ -101,6 +101,13 @@ res <- mmedist(x, "norm3", start=list(a=1), order=1, memp=function(x, order) mea
 attr(try(log("a"), silent=TRUE), "condition")
 
 
+try(mmedist(c(serving, "a"), "gamma"))
+try(mmedist(c(serving, NA), "gamma"))
+try(mmedist(c(serving, Inf), "gamma"))
+try(mmedist(c(serving, -Inf), "gamma"))
+try(mmedist(c(serving, NaN), "gamma"))
+
+
 # (7) fit of a normal distribution with weighted moment matching estimation
 #
 
diff --git a/tests/t-msedist.R b/tests/t-msedist.R
index 5f5d0b5f..448edeab 100644
--- a/tests/t-msedist.R
+++ b/tests/t-msedist.R
@@ -55,6 +55,18 @@ legend("bottomright", lty=1:2, col=c("green", "blue"), leg=c("theoretical value"
 #         upper=c(1,1,0, Inf))
 # pBMT(x, p3=1/2, p4=1/2, p1=-1/2, p2=1/2)
 
+
+
+
+
+try(msedist(c(x1, "a"), "gamma"))
+try(msedist(c(x1, NA), "gamma"))
+try(msedist(c(x1, Inf), "gamma"))
+try(msedist(c(x1, -Inf), "gamma"))
+try(msedist(c(x1, NaN), "gamma"))
+
+
+
 #--------------------------------------------------------
 # lognormal sample
 
diff --git a/tests/t-qmedist.R b/tests/t-qmedist.R
index b96b6215..976a5291 100644
--- a/tests/t-qmedist.R
+++ b/tests/t-qmedist.R
@@ -80,6 +80,16 @@ res <- qmedist(x, "norm3", start=list(a=1), probs=1/2)
 attr(try(log("a"), silent=TRUE), "condition")
 
 
+
+
+
+try(qmedist(c(x1, "a"), "gamma", probs=1/2))
+try(qmedist(c(x1, NA), "gamma", probs=1/2))
+try(qmedist(c(x1, Inf), "gamma", probs=1/2))
+try(qmedist(c(x1, -Inf), "gamma", probs=1/2))
+try(qmedist(c(x1, NaN), "gamma", probs=1/2))
+
+
 # (8) weighted QME
 #
 n <- 1e6