Running 't-CIcdfplot.R' [184s]
  Running 't-Surv2fitdistcens.R' [2s]
  Running 't-bootdist.R' [5s]
  Running 't-bootdistcens.R' [2s]
  Running 't-cdfcomp.R' [166s]
  Running 't-cdfcompcens.R' [173s]
  Running 't-cvg-algo.R' [169s]
  Running 't-denscomp.R' [4s]
  Running 't-descdist.R' [1s]
  Running 't-detectbound.R' [170s]
  Running 't-fitbench.R' [156s]
  Running 't-fitdist-customoptim.R' [156s]
  Running 't-fitdist.R' [3s]
  Running 't-fitdistcens.R' [157s]
  Running 't-gen-max-spacing-estim.R' [173s]
  Running 't-getparam.R' [1s]
  Running 't-gofstat.R' [2s]
  Running 't-init-actuar.R' [0s]
  Running 't-llplot.R' [176s]
  Running 't-lnL-surf.R' [2s]
  Running 't-logLik-vcov-coef.R' [166s]
  Running 't-manageparam.R' [2s]
  Running 't-mgedist.R' [161s]
  Running 't-mledist-cens.R' [166s]
  Running 't-mledist-nocens.R' [2s]
  Running 't-mledist-paramsupport.R' [3s]
  Running 't-mmedist.R' [4s]
  Running 't-msedist.R' [172s]
  Running 't-parallel.R' [156s]
  Running 't-plotdist.R' [2s]
  Running 't-plotdistcens.R' [169s]
  Running 't-ppcomp.R' [3s]
  Running 't-ppcompcens.R' [3s]
  Running 't-prefit.R' [170s]
  Running 't-qme-discrete.R' [177s]
  Running 't-qmedist.R' [179s]
  Running 't-qqcomp.R' [4s]
  Running 't-qqcompcens.R' [4s]
  Running 't-quantiledist.R' [2s]
  Running 't-startfixarg-overall.R' [3s]
  Running 't-startingvalues.R' [1s]
  Running 't-util-npmle.R' [1s]
  Running 't-util-npsurv-mainfunction.R' [164s]
  Running 't-util-testdensity.R' [171s]
  Running 't-weird-ppcomp-cens.R' [156s]
  Running 't-weird-qqcomp-cens.R' [154s]
Running the tests in 'tests/t-CIcdfplot.R' failed.
Complete output:
  > library(fitdistrplus)
  Loading required package: MASS
  Loading required package: survival
  > 
  > nbboot <- 201
  > nbboot <- 10
  > ggplotEx <- requireNamespace("ggplot2", quietly = TRUE)
  > 
  > # (1) Fit of a gamma distribution 
  > #
  > 
  > set.seed(123)
  > s1 <- rgamma(20, 3, 2)
  > f1 <- fitdist(s1, "gamma")
  > b1 <- bootdist(f1, niter=nbboot, silent=TRUE)
  > 
  > plot(b1)
  > quantile(b1)
  (original) estimated quantiles for each specified probability (non-censored data)
               p=0.1    p=0.2     p=0.3    p=0.4   p=0.5    p=0.6    p=0.7
  estimate 0.4742226 0.672897 0.8486169 1.022154 1.20499 1.408707 1.650739
              p=0.8    p=0.9
  estimate 1.966843 2.466269
  Median of bootstrap estimates
               p=0.1     p=0.2     p=0.3    p=0.4    p=0.5    p=0.6    p=0.7
  estimate 0.5083059 0.6866909 0.8522568 1.018742 1.200786 1.384358 1.577513
              p=0.8    p=0.9
  estimate 1.797666 2.213844
  
  two-sided 95 % CI of each quantile
             p=0.1     p=0.2     p=0.3     p=0.4     p=0.5    p=0.6    p=0.7
  2.5 %  0.3930866 0.5545204 0.6866101 0.8156534 0.9504778 1.099636 1.275705
  97.5 % 0.7636623 0.9342928 1.0731198 1.2027355 1.3333688 1.479356 1.762055
            p=0.8    p=0.9
  2.5 %  1.498104 1.827576
  97.5 % 2.151721 2.778154
  > 
  > par(mfrow=c(1,2))
  > CIcdfplot(b1, CI.level=95/100, CI.output = "probability", CI.fill="grey80", CI.col="black")
  > CIcdfplot(b1, CI.level=95/100, CI.output = "quantile", datacol="blue")
  > if(ggplotEx) CIcdfplot(b1, CI.level=95/100, CI.output = "probability", CI.fill="grey80", CI.col="black", plotstyle = "ggplot")
  > if(ggplotEx) CIcdfplot(b1, CI.level=95/100, CI.output = "quantile", datacol="blue", plotstyle = "ggplot")
  >   
  > par(mfrow=c(1,2))
  > CIcdfplot(b1, CI.level=85/100, CI.output = "probability")
  > CIcdfplot(b1, CI.level=85/100, CI.output = "quantile")
  > if(ggplotEx) CIcdfplot(b1, CI.level=85/100, CI.output = "probability", plotstyle = "ggplot")
  > if(ggplotEx) CIcdfplot(b1, CI.level=85/100, CI.output = "quantile", plotstyle = "ggplot")
  > 
  > par(mfrow=c(1,2))
  > CIcdfplot(b1, CI.level=90/100, CI.output = "probability")
  > CIcdfplot(b1, CI.level=90/100, CI.output = "quantile", CI.col="black", CI.type = "less", 
  +           CI.fill="grey85", verticals=TRUE, datacol="blue", do.points=FALSE)
  > if(ggplotEx) CIcdfplot(b1, CI.level=90/100, CI.output = "probability", plotstyle = "ggplot")
  > if(ggplotEx) CIcdfplot(b1, CI.level=90/100, CI.output = "quantile", CI.col="black", CI.type = "less", 
  +           CI.fill="grey85", verticals=TRUE, datacol="blue", do.points=FALSE, plotstyle = "ggplot")
  > 
  > par(mfrow=c(1,2))
  > CIcdfplot(b1, CI.level=90/100, CI.output = "probability", CI.type = "greater")
  > CIcdfplot(b1, CI.level=90/100, CI.output = "quantile", CI.col="black", CI.type = "greater", 
  +           CI.fill="grey90", datacol="blue", datapch=21)
  > if(ggplotEx) CIcdfplot(b1, CI.level=90/100, CI.output = "probability", CI.type = "greater", plotstyle = "ggplot")
  > if(ggplotEx) CIcdfplot(b1, CI.level=90/100, CI.output = "quantile", CI.col="black", CI.type = "greater", 
  +           CI.fill="grey90", datacol="blue", datapch=21, plotstyle = "ggplot")
  > 
  > par(mfrow=c(1,1))
  > CIcdfplot(b1, CI.level=90/100, CI.output = "probability", CI.col="black", CI.type = "less", CI.fill="grey90")
  > CIcdfplot(b1, CI.level=90/100, CI.output = "quantile", CI.col="black", CI.type = "less", CI.fill="grey90", 
  +           verticals=TRUE, datacol="blue", do.points=FALSE)
  > CIcdfplot(b1, CI.level=90/100, CI.output = "quantile", CI.col="grey90", CI.type = "less", CI.fill="grey90", 
  +           verticals=TRUE, datacol="blue", do.points=FALSE, CI.only=TRUE)
  > CIcdfplot(b1, CI.level=90/100, CI.output = "probability", CI.col="grey85", CI.type = "less", CI.fill="grey90", 
  +           CI.only = TRUE)
  > CIcdfplot(b1, CI.output = "probability", fitlty=3, fitlwd=4)
  > 
  > if(ggplotEx) CIcdfplot(b1, CI.level=90/100, CI.output = "probability", CI.col="black", CI.type = "less", CI.fill="grey90", plotstyle = "ggplot")
  > if(ggplotEx) CIcdfplot(b1, CI.level=90/100, CI.output = "quantile", CI.col="black", CI.type = "less", CI.fill="grey90", 
  +           verticals=TRUE, datacol="blue", do.points=FALSE, plotstyle = "ggplot")
  > if(ggplotEx) CIcdfplot(b1, CI.level=90/100, CI.output = "quantile", CI.col="grey90", CI.type = "less", CI.fill="grey90", 
  +           verticals=TRUE, datacol="blue", do.points=FALSE, CI.only=TRUE, plotstyle = "ggplot")
  > if(ggplotEx) CIcdfplot(b1, CI.level=90/100, CI.output = "probability", CI.col="grey85", CI.type = "less", CI.fill="grey90", 
  +           CI.only = TRUE, plotstyle = "ggplot")
  > if(ggplotEx) CIcdfplot(b1, CI.output = "probability", fitlty=3, fitlwd=4, plotstyle = "ggplot")
  > 
  > # (2) an example from ecotoxicology
  > # with censored data
  > #
  > data(salinity)
  > log10LC50 <-log10(salinity)
  > fln <- fitdistcens(log10LC50,"norm")
  > bln <- bootdistcens(fln, niter=nbboot)
  > (HC5ln <- quantile(bln,probs = 0.05))
  (original) estimated quantiles for each specified probability (censored data)
            p=0.05
  estimate 1.11584
  Median of bootstrap estimates
             p=0.05
  estimate 1.110188
  
  two-sided 95 % CI of each quantile
           p=0.05
  2.5 %  1.032713
  97.5 % 1.153806
  > CIcdfplot(bln, CI.output = "quantile", CI.fill = "lightblue", CI.col = "blue",
  +           xlab = "log10(LC50)",xlim=c(0.5,2),lines01 = TRUE)
  > if(ggplotEx) CIcdfplot(bln, CI.output = "quantile", CI.fill = "lightblue", CI.col = "blue",
  +           xlab = "log10(LC50)",xlim=c(0.5,2),lines01 = TRUE, plotstyle = "ggplot")
  > 
  > # zoom around the HC5 with CI on quantiles 
  > CIcdfplot(bln, CI.output = "quantile", CI.fill = "lightblue", CI.col = "blue",
  +           xlab = "log10(LC50)", lines01 = TRUE, xlim = c(0.8, 1.5), ylim = c(0, 0.1))
  > abline(h = 0.05, lty = 1)
  > if(ggplotEx) CIcdfplot(bln, CI.output = "quantile", CI.fill = "lightblue", CI.col = "blue",
  +           xlab = "log10(LC50)", lines01 = TRUE, xlim = c(0.8, 1.5), ylim = c(0, 0.1), plotstyle = "ggplot") + 
  +   ggplot2::geom_hline(yintercept = 0.05)
  > 
  > # zoom around the HC5 with CI on probabilities 
  > CIcdfplot(bln, CI.output = "probability", CI.fill = "lightblue", CI.col = "blue",
  +           xlab = "log10(LC50)", lines01 = TRUE, xlim = c(0.8, 1.5), ylim = c(0, 0.1))
  > abline(h = 0.05, lty = 1)
  > if(ggplotEx) CIcdfplot(bln, CI.output = "probability", CI.fill = "lightblue", CI.col = "blue",
  +           xlab = "log10(LC50)", lines01 = TRUE, xlim = c(0.8, 1.5), ylim = c(0, 0.1), plotstyle = "ggplot") + 
  +   ggplot2::geom_hline(yintercept = 0.05)
  > 
  >   
  > # (3) An example where the difference between "probability"
  > #     and "quantile" is clear on the plot
  > #
  > 
  > set.seed(123)
  > s3 <- rgamma(5, 3, 10)
  > f3 <- fitdist(s3, "norm")
  > b3 <- bootdist(f3, niter=nbboot, silent=TRUE)
  > 
  > par(mfrow=c(1,2))
  > CIcdfplot(b3, CI.level=90/100, CI.output = "probability")
  > CIcdfplot(b3, CI.level=90/100, CI.output = "quantile")
  > 
  > if(ggplotEx) CIcdfplot(b3, CI.level=90/100, CI.output = "probability", plotstyle = "ggplot")
  > if(ggplotEx) CIcdfplot(b3, CI.level=90/100, CI.output = "quantile", plotstyle = "ggplot")
  > 
  > #some ideas from http://edild.github.io/ssd/
  > 
  > proc.time()
     user  system elapsed 
    12.20    0.42   12.60 
Running the tests in 'tests/t-cdfcomp.R' failed.
Complete output:
  > library(fitdistrplus)
  Loading required package: MASS
  Loading required package: survival
  > 
  > # ?cdfcomp
  > visualize <- FALSE # TRUE for manual tests with visualization of results
  > nsample <- 1000
  > nsample <- 10
  > 
  > 
  > # (1) Plot various distributions fitted to serving size data
  > #
  > data(groundbeef)
  > serving <- groundbeef$serving
  > fitW <- fitdist(serving,"weibull")
  > fitln <- fitdist(serving,"lnorm")
  > fitg <- fitdist(serving,"gamma")
  > 
  > #sanity checks
  > try(cdfcomp("list(fitW, fitln, fitg)",horizontals = FALSE), silent=TRUE)
  > try(cdfcomp(list(fitW, fitln, fitg, a=1),horizontals = FALSE), silent=TRUE)
  > 
  > #real call
  > cdfcomp(list(fitW, fitln, fitg), horizontals = FALSE)
  > cdfcomp(list(fitW, fitln, fitg), horizontals = TRUE)
  > cdfcomp(list(fitW, fitln, fitg), horizontals = TRUE, lines01 = TRUE)
  > cdfcomp(list(fitW, fitln, fitg), horizontals = TRUE, verticals = TRUE, datacol = "grey")
  > if (requireNamespace ("ggplot2", quietly = TRUE)) {
  +   cdfcomp(list(fitW, fitln, fitg), horizontals = FALSE, plotstyle = "ggplot")
  + }
  > if (requireNamespace ("ggplot2", quietly = TRUE) & visualize) {
  +   cdfcomp(list(fitW, fitln, fitg), horizontals = TRUE, plotstyle = "ggplot")
  +   cdfcomp(list(fitW, fitln, fitg), horizontals = TRUE, lines01 = TRUE, plotstyle = "ggplot")
  +   cdfcomp(list(fitW, fitln, fitg), horizontals = TRUE, verticals = TRUE, datacol = "grey", plotstyle = "ggplot")
  + }
  > 
  > cdfcomp(list(fitW, fitln, fitg), legendtext = c("Weibull", "lognormal", "gamma"),
  +         main = "ground beef fits", xlab = "serving sizes (g)",
  +         ylab = "F(g)", xlim = c(0, 250), ylim = c(.5, 1))
  > if (requireNamespace ("ggplot2", quietly = TRUE) & visualize) {
  +   cdfcomp(list(fitW, fitln, fitg), legendtext = c("Weibull", "lognormal", "gamma"),
  +           main = "ground beef fits", xlab = "serving sizes (g)",
  +           ylab = "F(g)", xlim = c(0, 250), ylim = c(.5, 1), plotstyle = "ggplot")
  + }
  > 
  > cdfcomp(list(fitW, fitln, fitg), legendtext = c("Weibull", "lognormal", "gamma"),
  +         main = "ground beef fits", xlab = "serving sizes (g)",
  +         ylab = "F(g)", xlogscale = TRUE)
  > if (requireNamespace ("ggplot2", quietly = TRUE) & visualize) {
  +   cdfcomp(list(fitW, fitln, fitg), legendtext = c("Weibull", "lognormal", "gamma"),
  +           main = "ground beef fits", xlab = "serving sizes (g)",
  +           ylab = "F(g)", xlogscale = TRUE, plotstyle = "ggplot")
  + }
  > 
  > cdfcomp(list(fitW,fitln,fitg),legendtext=c("Weibull","lognormal","gamma"),
  +         main="ground beef fits",xlab="serving sizes (g)",
  +         ylab="F(g)", xlogscale=TRUE, ylogscale=TRUE, ylim=c(.005, .99))
  > if (requireNamespace ("ggplot2", quietly = TRUE) & visualize) {
  +   cdfcomp(list(fitW,fitln,fitg),legendtext=c("Weibull","lognormal","gamma"),
  +           main="ground beef fits",xlab="serving sizes (g)",
  +           ylab="F(g)", xlogscale=TRUE, ylogscale=TRUE, ylim=c(.005, .99), plotstyle = "ggplot")
  + }
  > 
  > cdfcomp(list(fitW,fitln,fitg), legendtext=c("Weibull","lognormal","gamma"),
  +         main="ground beef fits",xlab="serving sizes (g)",
  +         ylab="F(g)",xlim = c(0,250), xlegend = "topleft")
  > if (requireNamespace ("ggplot2", quietly = TRUE) & visualize) {
  +   cdfcomp(list(fitW,fitln,fitg), legendtext=c("Weibull","lognormal","gamma"),
  +           main="ground beef fits",xlab="serving sizes (g)",
  +           ylab="F(g)",xlim = c(0,250), xlegend = "topleft", plotstyle = "ggplot")
  + }
  > 
  > 
  > 
  > 
  > # (2) Plot lognormal distributions fitted by 
  > # maximum goodness-of-fit estimation
  > # using various distances (data plotted in log scale)
  > #
  > data(endosulfan)
  > ATV <-subset(endosulfan, group == "NonArthroInvert")$ATV
  > taxaATV <- subset(endosulfan, group == "NonArthroInvert")$taxa
  > flnMGEKS <- fitdist(ATV,"lnorm",method="mge",gof="KS")
  > flnMGEAD <- fitdist(ATV,"lnorm",method="mge",gof="AD")
  > flnMGEADL <- fitdist(ATV,"lnorm",method="mge",gof="ADL")
  > flnMGEAD2L <- fitdist(ATV,"lnorm",method="mge",gof="AD2L")
  > llfit <- list(flnMGEKS, flnMGEAD, flnMGEADL, flnMGEAD2L)
  > 
  > cdfcomp(list(flnMGEKS, flnMGEAD, flnMGEADL, flnMGEAD2L), xlogscale=TRUE, 
  +         main="fits of a lognormal dist. using various GOF dist.")
  > cdfcomp(list(flnMGEKS, flnMGEAD, flnMGEADL, flnMGEAD2L),
  +         xlogscale=TRUE,main="fits of a lognormal dist. using various GOF dist.",
  +         legendtext=c("MGE KS","MGE AD","MGE ADL","MGE AD2L"))
  > cdfcomp(list(flnMGEKS, flnMGEAD, flnMGEADL, flnMGEAD2L),
  +         xlogscale=TRUE,main="fits of a lognormal dist. using various GOF dist.",
  +         legendtext=c("MGE KS","MGE AD","MGE ADL","MGE AD2L"), 
  +         fitcol=c("black", "darkgreen", "yellowgreen", "yellow2"), 
  +         horizontals=FALSE, datapch="+")
  > cdfcomp(llfit, xlogscale=TRUE, main="fits of a lognormal dist. using various GOF dist.",
  +         legendtext=paste("MGE", c("KS","AD","ADL","AD2L")), fitcol="grey35", 	
  +         fitlty="dotted", horizontals=FALSE, datapch=21, datacol="grey30")
  > cdfcomp(list(flnMGEKS, flnMGEAD, flnMGEADL, flnMGEAD2L),
  +         xlogscale=TRUE, verticals=TRUE, xlim=c(10,100000), datapch=21)
  > cdfcomp(flnMGEKS, xlogscale=TRUE, verticals=TRUE, xlim=c(10,100000))
  > cdfcomp(list(flnMGEKS, flnMGEAD, flnMGEADL, flnMGEAD2L),
  +         xlogscale=TRUE, verticals=TRUE, xlim=c(1,100000), datapch=21, name.points=taxaATV)
  > cdfcomp(flnMGEKS, xlogscale=TRUE, verticals=TRUE, xlim=c(1,100000), name.points=taxaATV)
  > 
  > if (requireNamespace ("ggplot2", quietly = TRUE) & visualize) {
  +   cdfcomp(list(flnMGEKS, flnMGEAD, flnMGEADL, flnMGEAD2L), xlogscale=TRUE,
  +           main = "fits of a lognormal dist. using various GOF dist.", plotstyle = "ggplot")
  +   cdfcomp(list(flnMGEKS, flnMGEAD, flnMGEADL, flnMGEAD2L),
  +           xlogscale=TRUE,main="fits of a lognormal dist. using various GOF dist.",
  +           legendtext=c("MGE KS","MGE AD","MGE ADL","MGE AD2L"), plotstyle = "ggplot")
  +   cdfcomp(list(flnMGEKS, flnMGEAD, flnMGEADL, flnMGEAD2L),
  +           xlogscale=TRUE,main="fits of a lognormal dist. using various GOF dist.",
  +           legendtext=c("MGE KS","MGE AD","MGE ADL","MGE AD2L"), 
  +           fitcol=c("black", "darkgreen", "yellowgreen", "yellow2"), 
  +           horizontals=FALSE, datapch="+", plotstyle = "ggplot")
  +   cdfcomp(llfit, xlogscale=TRUE, main="fits of a lognormal dist. using various GOF dist.",
  +           legendtext=paste("MGE", c("KS","AD","ADL","AD2L")), fitcol="grey35", 	
  +           fitlty="dotted", horizontals=FALSE, datapch=21, datacol="grey30", plotstyle = "ggplot")
  +   cdfcomp(list(flnMGEKS, flnMGEAD, flnMGEADL, flnMGEAD2L),
  +           xlogscale=TRUE, verticals=TRUE, xlim=c(10,100000), datapch=21, plotstyle = "ggplot")
  +   cdfcomp(flnMGEKS, xlogscale=TRUE, verticals=TRUE, xlim=c(10,100000), plotstyle = "ggplot")
  +   cdfcomp(list(flnMGEKS, flnMGEAD, flnMGEADL, flnMGEAD2L),
  +           xlogscale=TRUE, verticals=TRUE, xlim=c(1,100000), datapch=21, name.points=taxaATV, 
  +           plotstyle = "ggplot")
  +   cdfcomp(flnMGEKS, xlogscale=TRUE, verticals=TRUE, xlim=c(1,100000), name.points=taxaATV, 
  +           plotstyle = "ggplot")
  + }
  > 
  > 
  > 
  > 
  > # (3) Plot normal and logistic distributions fitted by 
  > # maximum likelihood estimation
  > # using various plotting positions in cdf plots
  > #
  > if (visualize)
  + {
  +   data(endosulfan)
  +   log10ATV <-log10(subset(endosulfan, group == "NonArthroInvert")$ATV)
  +   taxaATV <- subset(endosulfan, group == "NonArthroInvert")$taxa
  +   fln <- fitdist(log10ATV, "norm")
  +   fll <- fitdist(log10ATV, "logis")
  +   
  +   # default plot using Hazen plotting position: (1:n - 0.5)/n
  +   cdfcomp(list(fln, fll), legendtext = c("normal", "logistic"), xlab = "log10ATV", 
  +           name.points=taxaATV, xlim = c(0,5))
  +   if (requireNamespace ("ggplot2", quietly = TRUE)) {
  +     cdfcomp(list(fln, fll), legendtext = c("normal", "logistic"),xlab = "log10ATV", 
  +             name.points=taxaATV, xlim=c(0,5), plotstyle = "ggplot")
  +   }
  +   
  +   # plot using mean plotting position (named also Gumbel plotting position)
  +   # (1:n)/(n + 1)
  +   cdfcomp(list(fln,fll),legendtext=c("normal","logistic"),xlab="log10ATV",
  +           use.ppoints = TRUE, a.ppoints = 0)
  +   if (requireNamespace ("ggplot2", quietly = TRUE)) {
  +     cdfcomp(list(fln,fll),legendtext=c("normal","logistic"),xlab="log10ATV",
  +             use.ppoints = TRUE, a.ppoints = 0, plotstyle = "ggplot")
  +   }
  +   
  +   # plot using basic plotting position: (1:n)/n
  +   cdfcomp(list(fln,fll),legendtext=c("normal","logistic"),xlab="log10ATV",
  +           use.ppoints = FALSE)
  +   if (requireNamespace ("ggplot2", quietly = TRUE)) {
  +     cdfcomp(list(fln,fll),legendtext=c("normal","logistic"),xlab="log10ATV",
  +             use.ppoints = FALSE, plotstyle = "ggplot")
  +   }
  + }
  > 
  > 
  > 
  > 
  > # (4) Plot lognormal distributions fitted by 
  > # maximum goodness-of-fit estimation
  > # using various distances (data plotted in log scale)
  > #
  > if (visualize)
  + {
  +   x1 <- c(6.4,13.3,4.1,1.3,14.1,10.6,9.9,9.6,15.3,22.1,13.4,13.2,8.4,6.3,8.9,5.2,10.9,14.4)
  +   x <- seq(0, 1.1*max(x1), length=100)
  +   
  +   dgumbel <- function(x,a,b) 1/b*exp((a-x)/b)*exp(-exp((a-x)/b))
  +   pgumbel <- function(q,a,b) exp(-exp((a-q)/b))
  +   
  +   f1 <- mledist(x1, "norm")
  +   f2 <- mledist(x1, "gumbel", start = list(a = 10, b = 5))
  +   f3 <- mledist(x1, "exp")
  +   
  +   # graph 1
  +   plot(ecdf(x1))
  +   lines(x, pnorm(x, f1$estimate[1], f1$estimate[2]), col = "red")
  +   lines(x, pgumbel(x, f2$estimate[1], f2$estimate[2]), col = "green")
  +   lines(x, pexp(x, f3$estimate[1]), col = "blue")
  +   legend("bottomright", lty = 1, leg = c("Normal", "Gumbel", "Exp"), col = c("red", "green", "blue"))
  +   
  +   # graph 2
  +   f1 <- fitdist(x1, "norm")
  +   f2 <- fitdist(x1, "gumbel", start = list(a = 10, b = 5))
  +   f3 <- fitdist(x1, "exp")
  +   cdfcomp(list(f1, f2, f3), xlim=range(x), fitcol = c("red", "green", "blue"), fitlty = 1, legendtext = c("Normal", "Gumbel", "Exp"))
  +   cdfcomp(list(f1, f2, f3), xlim=range(x), fitcol = c("red", "green", "blue"), fitlty = 1, fitlwd = (1:3)*2, legendtext = c("Normal", "Gumbel", "Exp"))
  +   
  +   # graph 3
  +   if (requireNamespace ("ggplot2", quietly = TRUE)) {
  +     cdfcomp(list(f1, f2, f3), xlim=range(x), fitcol=c("red","green","blue"), fitlty = 1, legendtext = c("Normal", "Gumbel", "Exp"), plotstyle = "ggplot")
  +     cdfcomp(list(f1, f2, f3), xlim=range(x), fitcol=c("red","green","blue"), fitlty = 1, fitlwd = 1:3, legendtext = c("Normal", "Gumbel", "Exp"), plotstyle = "ggplot")
  +   }
  + }
  > 
  > # (5) normal mixture
  > #
  > 
  > # mixture of two normal distributions
  > # density
  > dnorm2 <- function(x, poid, m1, s1, m2, s2)
  +   poid*dnorm(x, m1, s1) + (1-poid)*dnorm(x, m2, s2)
  > # numerically approximate quantile function
  > qnorm2 <- function(p, poid, m1, s1, m2, s2)
  + {
  +   L2 <- function(x, prob)
  +     (prob - pnorm2(x, poid, m1, s1, m2, s2))^2	
  +   sapply(p, function(pr) optimize(L2, c(-1000, 1000), prob=pr)$minimum)
  + }	
  > # distribution function		
  > pnorm2 <- function(q, poid, m1, s1, m2, s2)
  +   poid*pnorm(q, m1, s1) + (1-poid)*pnorm(q, m2, s2)		
  > 
  > 
  > # basic normal distribution
  > set.seed(1234)
  > x2 <- c(rnorm(nsample, 5),  rnorm(nsample, 10))
  > # MLE fit
  > fit1 <- fitdist(x2, "norm2", "mle", start=list(poid=1/3, m1=4, s1=2, m2=8, s2=2), 
  +                 lower=c(0, 0, 0, 0, 0))
  > fit2 <- fitdist(x2, "norm2", "qme", probs=c(1/6, 1/4, 1/3, 1/2, 2/3), 
  +                 start=list(poid=1/3, m1=4, s1=2, m2=8, s2=2), 
  +                 lower=c(0, 0, 0, 0, 0), upper=c(1/2, Inf, Inf, Inf, Inf))
  > fit3 <- fitdist(x2, "norm2", "mge", gof="AD", 
  +                 start=list(poid=1/3, m1=4, s1=2, m2=8, s2=2), 
  +                 lower=c(0, 0, 0, 0, 0), upper=c(1/2, Inf, Inf, Inf, Inf))
  > 
  > cdfcomp(list(fit1, fit2, fit3), datapch=".")
  > cdfcomp(list(fit1, fit2, fit3), datapch=".", xlim=c(6, 8), ylim=c(.4, .55))
  > 
  > if (requireNamespace ("ggplot2", quietly = TRUE) & visualize) {
  +   cdfcomp(list(fit1, fit2, fit3), datapch=".", plotstyle = "ggplot")
  +   cdfcomp(list(fit1, fit2, fit3), datapch=".", xlim=c(6, 8), ylim=c(.4, .55), plotstyle = "ggplot")
  + }
  > 
  > 
  > 
  > 
  > # (6) discrete example
  > #
  > set.seed(1234)
  > x3 <- rpois(20, 10)
  > 
  > fit1 <- fitdist(x3, "pois", "mle")
  > fit2 <- fitdist(x3, "nbinom", "qme", probs=c(1/3, 2/3))
  > 
  > cdfcomp(list(fit1, fit2), datapch=21, horizontals=FALSE)
  > if (requireNamespace ("ggplot2", quietly = TRUE) & visualize) {
  +   cdfcomp(list(fit1, fit2), datapch=21, horizontals=FALSE, plotstyle = "ggplot")
  + }
  > 
  > 
  > # (7) large dataset
  > #
  > if (visualize)
  + {
  +   n <- 1e4
  +   f1 <- fitdist(rlnorm(n), "lnorm")
  +   
  +   cdfcomp(f1, do.points=TRUE)
  +   cdfcomp(f1, do.points=FALSE)
  +   cdfcomp(f1, horizontals = FALSE, verticals = FALSE, do.points = FALSE)
  +   if (requireNamespace ("ggplot2", quietly = TRUE) & visualize) {
  +     cdfcomp(f1, do.points = TRUE, plotstyle = "ggplot")
  +     cdfcomp(f1, do.points = FALSE, plotstyle = "ggplot")
  +     cdfcomp(f1, horizontals = FALSE, verticals = FALSE, do.points = FALSE, plotstyle = "ggplot")
  +   }
  +   
  + }
  > 
  > # (8) argument add (must give the same plot (except colors) as ex. 6)
  > #
  > set.seed(1234)
  > x3 <- rpois(nsample, 10)
  > 
  > fit1 <- fitdist(x3, "pois", "mle")
  > cdfcomp(fit1, fitcol = "red", horizontals=FALSE, addlegend = FALSE)
  > fit2 <- fitdist(x3, "nbinom", "qme", probs=c(1/3, 2/3))
  > cdfcomp(fit2, fitcol = "blue", horizontals=FALSE, addlegend = FALSE, add = TRUE)
  > cdfcomp(list(fit1, fit2), horizontals=FALSE, addlegend = FALSE, 
  +         fitcol=c("red", "blue"))
  > 
  > 
  > if (requireNamespace ("ggplot2", quietly = TRUE) & visualize) {
  +   # the argument add is not available when plotstyle = "ggplot"
  +   cdfcomp(list(fit1, fit2), fitcol = c("red", "blue"), fitlty = 1, horizontals = FALSE, addlegend = FALSE, plotstyle = "ggplot")
  + }
  > 
  > 
  > # (9) test legend labels
  > #
  > if (visualize)
  + {
  +   serving <- groundbeef$serving
  +   fitW <- fitdist(serving,"weibull")
  +   fitW2 <- fitdist(serving,"weibull", method="qme", probs=c(1/3,2/3))
  +   fitW3 <- fitdist(serving,"weibull", method="qme", probs=c(1/2,2/3))
  +   fitln <- fitdist(serving,"lnorm")
  +   fitg <- fitdist(serving,"gamma")
  +   
  +   cdfcomp(list(fitW, fitln, fitg)) #distrib
  +   cdfcomp(list(fitW, fitW2, fitln, fitg)) #distrib+method
  +   cdfcomp(list(fitW, fitW2, fitW3, fitln, fitg)) #distrib+method+num
  +   if (requireNamespace ("ggplot2", quietly = TRUE))
  +     cdfcomp(list(fitW, fitW2, fitW3, fitln, fitg), plotstyle = "ggplot") #distrib+method+num
  +   
  + }
  > 
  > proc.time()
     user  system elapsed 
     2.67    0.37    3.03 
Running the tests in 'tests/t-cdfcompcens.R' failed.
Complete output:
  > library(fitdistrplus)
  Loading required package: MASS
  Loading required package: survival
  > 
  > visualize <- FALSE # TRUE for manual tests with visualization of results
  > 
  > # (1) Plot various distributions fitted to bacterial contamination data
  > #
  > data(smokedfish)
  > Clog10 <- log10(smokedfish)
  > 
  > fitsfn <- fitdistcens(Clog10,"norm")
  > fitsfl <- fitdistcens(Clog10,"logis")
  > 
  > dgumbel <- function(x,a,b) 1/b*exp((a-x)/b)*exp(-exp((a-x)/b))
  > pgumbel <- function(q,a,b) exp(-exp((a-q)/b))
  > qgumbel <- function(p,a,b) a-b*log(-log(p))
  > fitsfg <- fitdistcens(Clog10, "gumbel", start=list(a=-3,b=3))
  > 
  > cdfcompcens(list(fitsfn,fitsfl,fitsfg))
  > cdfcompcens(list(fitsfn,fitsfl,fitsfg), fitlty=1, fitlwd=3)
  > 
  > # Same plot in y logscale
  > cdfcompcens(list(fitsfn, fitsfl, fitsfg), NPMLE.method = "Turnbull",
  +             ylogscale = TRUE, ylim=c(.5, .99))
  Warning message:
  In cdfcompcens(list(fitsfn, fitsfl, fitsfg), NPMLE.method = "Turnbull",  :
    Turnbull is now a deprecated option for NPMLE.method. You should use Turnbull.middlepoints
              of Turnbull.intervals. It was here fixed as Turnbull.middlepoints, equivalent to former Turnbull.
  > cdfcompcens(list(fitsfn, fitsfl, fitsfg), NPMLE.method = "Wang",
  +             ylogscale = TRUE, ylim=c(.5, .99))
  > 
  > if (requireNamespace ("ggplot2", quietly = TRUE)) {
  +   cdfcompcens(list(fitsfn,fitsfl,fitsfg), plotstyle = "ggplot")
  +   cdfcompcens(list(fitsfn,fitsfl,fitsfg), plotstyle = "ggplot", fitlty=1, fitlwd=3)
  +   cdfcompcens(list(fitsfn,fitsfl,fitsfg), datacol="grey", 
  +               legendtext=c("normal","logistic","Gumbel"),
  +               main="bacterial contamination fits",
  +               xlab="bacterial concentration (CFU/g)", ylab="F",
  +               xlegend = "center", lines01 = TRUE, plotstyle = "ggplot")
  +   
  +   # Same plot in y logscale
  +   cdfcompcens(list(fitsfn, fitsfl, fitsfg), NPMLE.method = "Wang",
  +               ylogscale = TRUE, ylim=c(.5, .99), plotstyle = "ggplot")
  + }
  Warning messages:
  1: In ggplot2::scale_y_continuous(trans = "log10") :
    log-10 transformation introduced infinite values.
  2: In ggplot2::scale_y_continuous(trans = "log10") :
    log-10 transformation introduced infinite values.
  3: In ggplot2::scale_y_continuous(trans = "log10") :
    log-10 transformation introduced infinite values.
  4: In ggplot2::scale_y_continuous(trans = "log10") :
    log-10 transformation introduced infinite values.
  > 
  > # Use of x logscale
  > if (visualize)
  + {
  +   if(any(installed.packages()[,"Package"] == "actuar"))
  +   {
  +     require(actuar)
  +     data(smokedfish)
  +     fln <- fitdistcens(smokedfish,"lnorm")
  +     fll <- fitdistcens(smokedfish,"llogis")
  +     cdfcompcens(list(fln, fll))
  +     cdfcompcens(list(fln, fll), xlogscale = TRUE)
  +     cdfcompcens(list(fln, fll), xlogscale = TRUE, xlim = c(0.01, 1000))
  +     cdfcompcens(list(fln, fll), NPMLE.method = "Turnbull",
  +                 xlogscale = TRUE, xlim = c(0.01, 1000))
  +     
  +     if (requireNamespace ("ggplot2", quietly = TRUE)) {
  +       cdfcompcens(list(fln, fll), plotstyle = "ggplot")
  +       cdfcompcens(list(fln, fll), xlogscale = TRUE, plotstyle = "ggplot")
  +       cdfcompcens(list(fln, fll), xlogscale = TRUE, xlim = c(0.01, 1000), plotstyle = "ggplot")
  +     }
  +   }
  +   # same plot using argument add
  +   cdfcompcens(fitsfn, addlegend = FALSE, fitcol = "red")
  +   cdfcompcens(fitsfl, addlegend = FALSE, fillrect = NA, fitcol = "green", add = TRUE)
  +   cdfcompcens(fitsfg, addlegend = FALSE, fillrect = NA, fitcol = "blue", add = TRUE)
  +   
  +   cdfcompcens(list(fitsfn, fitsfl, fitsfg), addlegend = FALSE, fitcol = 2:4, fitlty = 1, plotstyle = "ggplot")
  +   cdfcompcens(list(fitsfn, fitsfl, fitsfg), addlegend = FALSE, fitcol = 2:4, fitlty = 1, plotstyle = "ggplot")
  + }
  > 
  > # Test on the salinity data set
  > #
  > data(salinity)
  > log10LC50 <-log10(salinity)
  > plotdistcens(log10LC50)
  > plotdistcens(log10LC50, NPMLE = FALSE)
  > fn <- fitdistcens(log10LC50,"norm")
  > fl <- fitdistcens(log10LC50,"logis")
  > plot(fn)
  > plot(fl)
  > cdfcompcens(list(fn, fl))
  > 
  > if (requireNamespace("ggplot2", quietly = TRUE)) {
  +   cdfcompcens(list(fn, fl), plotstyle = "ggplot")
  + }
  > 
  > 
  > proc.time()
     user  system elapsed 
     3.46    0.50    3.98 
Running the tests in 'tests/t-cvg-algo.R' failed.
Complete output:
  > 
  > if(FALSE)
  + {  
  +   library(fitdistrplus)
  +   
  +   
  +   
  +   #(1) beta distribution
  +   #
  +   
  +   n <- 100
  +   set.seed(12345)
  +   x <- rbeta(n, 3, 3/4)
  +   
  +   psi <- function(x) digamma(x)
  +   grbetalnl <- function(x, a, b)
  +     c(log(x)-psi(a)+psi(a+b), log(1-x)-psi(b)+psi(a+b))
  +   #grbetalnl(x, 3, 4)
  +   grlnL <- function(par, obs, ...) 
  +     -rowSums(sapply(obs, function(x) grbetalnl(x, a=par[1], b=par[2])))
  +   #rowSums(sapply(x, function(x) grbetalnl(x, 3, 4)))
  +   #grlnL(c(3, 4), x)
  +   #grlnL(c(3, 3/4), x)
  +   
  +   
  +   constrOptim2 <- function(par, fn, gr=NULL, ui, ci, ...)
  +     constrOptim(theta=unlist(par), f=fn, grad=gr, ui=ui, ci=ci, ...)
  +   
  +   
  +   
  +   #control parameters
  +   ctr <- list(trace=3, REPORT=1, maxit=1000)
  +   ctr <- list(trace=0, REPORT=1, maxit=1000)
  +   
  +   bfgs_gr$time <- system.time(bfgs_gr <- mledist(x, dist="beta", optim.method="BFGS", gr=grlnL, control=ctr))[3]
  +   bfgs <- mledist(x, dist="beta", optim.method="BFGS", control=ctr)
  +   
  +   lbfgs_gr <- mledist(x, dist="beta", optim.method="L-BFGS-B", gr=grlnL, control=ctr, lower=c(0,0))
  +   lbfgs <- mledist(x, dist="beta", optim.method="L-BFGS-B", control=ctr, lower=c(0,0))
  +   
  +   cg_gr <- mledist(x, dist="beta", optim.method="CG", gr=grlnL, control=ctr)
  +   cg <- mledist(x, dist="beta", optim.method="CG", control=ctr)
  +   
  +   nm_gr <- mledist(x, dist="beta", optim.method="Nelder", gr=grlnL, control=ctr)
  +   nm <- mledist(x, dist="beta", optim.method="Nelder", control=ctr)
  +   
  +   constr_nm_gr <- mledist(x, dist="beta", custom.optim=constrOptim2,
  +                           ui = diag(2), ci = c(0, 0), optim.method="Nelder", gr=grlnL, control=ctr)
  +   constr_nm <- mledist(x, dist="beta", custom.optim=constrOptim2,
  +                        ui = diag(2), ci = c(0, 0), optim.method="Nelder", control=ctr)
  +   
  +   constr_bfgs_gr <- mledist(x, dist="beta", custom.optim=constrOptim2,
  +                             ui = diag(2), ci = c(0, 0), optim.method="BFGS", gr=grlnL, control=ctr)
  +   constr_bfgs <- mledist(x, dist="beta", custom.optim=constrOptim2,
  +                          ui = diag(2), ci = c(0, 0), optim.method="BFGS", control=ctr)
  +   
  +   constr_cg_gr <- mledist(x, dist="beta", custom.optim=constrOptim2,
  +                           ui = diag(2), ci = c(0, 0), optim.method="CG", gr=grlnL, control=ctr)
  +   constr_cg <- mledist(x, dist="beta", custom.optim=constrOptim2,
  +                        ui = diag(2), ci = c(0, 0), optim.method="CG", control=ctr)
  +   
  +   lnL <- function(par, fix.arg, obs, ddistnam) 
  +   {
  +     fitdistrplus:::loglikelihood(par, fix.arg, obs, ddistnam, weights = rep(1, NROW(obs))) 
  +   }
  +   
  +   
  +   constrOptim2(c(shape1=1, shape2=1), lnL, obs=x, fix.arg=NULL, ddistnam="dbeta",
  +                ui = diag(2), ci = c(0, 0))
  +   
  +   
  +   #no log
  +   dbeta3 <- function(x, shape1, shape2)
  +     dbeta(x, shape1, shape2)
  +   
  +   #Ripley trick : param transform
  +   dbeta2 <- function(x, shape1, shape2, log)
  +     dbeta(x, exp(shape1), exp(shape2), log=log)
  +   pbeta2 <- function(q, shape1, shape2, log.p)
  +     pbeta(q, exp(shape1), exp(shape2), log.p=log.p)
  +   
  +   bfgs2 <- mledist(x, dist="beta2", optim.method="BFGS", control=ctr, 
  +                    start=list(shape1=0, shape2=0))
  +   bfgs3 <- mledist(x, dist="beta3", optim.method="BFGS", control=ctr, 
  +                    start=list(shape1=1, shape2=1))
  +   
  +   
  +   
  +   
  +   getval <- function(x)
  +     c(x$estimate, loglik=x$loglik, x$counts)
  +   getval2 <- function(x)
  +     c(exp(x$estimate), loglik=x$loglik, x$counts)
  +   
  +   cbind(trueval=c(3, 3/4, lnL(c(3, 3/4), NULL, x, "dbeta"), NA, NA),
  +         NM=getval(nm), NMgrad=getval(nm_gr), 
  +         constr_NM=getval(constr_nm), constr_NMgrad=getval(constr_nm_gr),
  +         CG=getval(cg), CGgrad=getval(cg_gr), 
  +         constr_CG=getval(constr_cg), constr_CGgrad=getval(constr_cg_gr),
  +         BFGS=getval(bfgs), BFGSgrad=getval(bfgs_gr),
  +         constr_BFGS=getval(constr_bfgs), constr_BFGSgrad=getval(constr_bfgs_gr),
  +         BFGS_exp=getval2(bfgs2), BFGS_nolog=getval(bfgs3))
  +   
  +   
  +   
  +   llsurface(min.arg=c(0.1, 0.1), max.arg=c(7, 3), plot.arg=c("shape1", "shape2"),
  +             lseq=50, data=x, distr="beta")
  +   points(bfgs$estimate[1], bfgs$estimate[2], pch="+", col="red")
  +   points(3, 3/4, pch="x", col="green")
  +   
  + }
  > 
  > proc.time()
     user  system elapsed 
     0.26    0.01    0.26 
Running the tests in 'tests/t-detectbound.R' failed.
Complete output:
  > require(fitdistrplus)
  Loading required package: fitdistrplus
  Loading required package: MASS
  Loading required package: survival
  > 
  > #case where the density returns a Not-an-Numeric value.
  > detectbound("gamma", c(shape=3, scale=3), 1:10, echo=TRUE)
  [1] "shape" "scale"
  shape scale 
      3     3 
  shape 
  1 	tested value -1 3 
   [1] NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
   [1] NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
  2 	tested value 1.490116e-08 3 
   [1] NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
   [1] 1.067715e-08 3.825256e-09 1.827277e-09 9.819759e-10 5.628932e-10
   [6] 3.361088e-10 2.064279e-10 1.294230e-10 8.243168e-11 5.315839e-11
  3 	tested value 1 3 
   [1] 0.23884377 0.17113904 0.12262648 0.08786571 0.06295853 0.04511176
   [7] 0.03232399 0.02316115 0.01659569 0.01189133
   [1] 0.23884377 0.17113904 0.12262648 0.08786571 0.06295854 0.04511176
   [7] 0.03232399 0.02316115 0.01659569 0.01189133
  scale 
  1 	tested value 3 -1 
   [1] NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
   [1] NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
  2 	tested value 3 1.490116e-08 
   [1] NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
   [1] 0 0 0 0 0 0 0 0 0 0
  3 	tested value 3 1 
   [1] 0.183939726 0.270670571 0.224041808 0.146525109 0.084224335 0.044617537
   [7] 0.022341107 0.010734803 0.004998097 0.002269996
   [1] 0.183939715 0.270670562 0.224041808 0.146525113 0.084224340 0.044617541
   [7] 0.022341109 0.010734805 0.004998098 0.002269997
       shape scale
  lowb     0     0
  uppb   Inf   Inf
  > detectbound("logis", c(location=3, scale=3), 1:10)
       location scale
  lowb     -Inf     0
  uppb      Inf   Inf
  > detectbound("geom", c(prob=1/2), 1:10)
       prob
  lowb    0
  uppb    1
  > 
  > #test rate-scale arg
  > detectbound("gamma", c(shape=1, scale=3), 1:10)
       shape scale
  lowb     0     0
  uppb   Inf   Inf
  > detectbound("gamma", c(shape=1, rate=1/3), 1:10)
       shape scale
  lowb     0     0
  uppb   Inf   Inf
  > 
  > 
  > #case where the density returns a Not-an-Numeric value and one parameter is fixed.
  > detectbound("gamma", c(shape=3), 1:10, fix.arg=c(scale=3))
       shape
  lowb     0
  uppb   Inf
  > 
  > 
  > #case where the density returns an error rather than a Not-an-Numeric value.
  > dgeom2 <- function(x, prob, log=FALSE)
  + {
  +   stopifnot(prob >= 0 && prob <= 1)
  +   dgeom(x, prob, log)
  + }
  > detectbound("geom2", c(prob=1/2), 1:10)
       prob
  lowb    0
  uppb    1
  > 
  > #case where the density returns a Not-an-Numeric value for actuar package
  > require(actuar)
  Loading required package: actuar
  
  Attaching package: 'actuar'
  
  The following objects are masked from 'package:stats':
  
      sd, var
  
  The following object is masked from 'package:grDevices':
  
      cm
  
  > 
  > detectbound("burr", c(shape1=3, shape2=3, rate=1), 1:10)
       shape1 shape2 scale
  lowb      0      0     0
  uppb    Inf    Inf   Inf
  > detectbound("llogis", c(shape=3, rate=1), 1:10)
       shape scale
  lowb     0     0
  uppb   Inf   Inf
  > 
  > proc.time()
     user  system elapsed 
     1.35    0.17    1.51 
Running the tests in 'tests/t-fitbench.R' failed.
Complete output:
  > if(FALSE)
  + {  
  +   require(fitdistrplus)
  +   
  +   x <- rgamma(1e3, shape=3/2, rate= 1/2)
  +   
  +   initval <- unlist(fitdistrplus:::startargdefault(x, "gamma"))
  +   
  +   fitdistrplus:::fitbench(x, "gamma", "mle")
  +   
  +   fitdistrplus:::fitbench(x, "gamma", "mle", lower=0)
  +   
  +   
  +   # grgam(1:2, x)
  +   # grgamlnl(x[1], shape=1, rate=2)
  +   # grgamlnl(x[1:2], shape=1, rate=2)
  +   
  +   fitdistrplus:::fitbench(x, "gamma", "mle", grad=fitdistrplus:::grlnlgamma)
  +   
  +   fitdistrplus:::fitbench(x, "gamma", "mle", grad=fitdistrplus:::grlnlgamma, lower=0)
  +   
  +   #mledist(x, "gamma", grad=grgam, lower=0, optim.method = "CG")
  +   #mledist(x, "gamma", grad=grgam, lower=0, optim.method = "BFGS")
  +   
  + }
  > 
  > proc.time()
     user  system elapsed 
     0.14    0.07    0.21 
Running the tests in 'tests/t-fitdist-customoptim.R' failed.
Complete output:
  > visualize <- FALSE # TRUE for manual tests with visualization of results
  > require(fitdistrplus)
  Loading required package: fitdistrplus
  Loading required package: MASS
  Loading required package: survival
  > 
  > myoptimize <- function(fn,par,ui,ci,...){
  +   res <- constrOptim(f=fn,theta=par,method="Nelder-Mead",ui=ui,ci=ci, ...)
  +   standardres <- c(res,convergence=0,value=res$objective,par=res$minimum,hessian=NA)
  +   return(standardres)
  + }
  > 
  > 
  > #one parameter example
  > x <- rexp(100)
  > 
  > #binding example 
  > fitdist(x, "exp", custom.optim=myoptimize, ui=1, ci=2, start=list(rate=10))
  Fitting of the distribution ' exp ' by maximum likelihood 
  Parameters:
       estimate Std. Error
  rate        2         NA
  > fitdist(x, "exp", lower= 2, optim.method="L-BFGS-B")
  Fitting of the distribution ' exp ' by maximum likelihood 
  Parameters:
       estimate Std. Error
  rate        2        0.2
  > 
  > #two parameter example
  > 
  > if(visualize) {  # check ERROR on aarch64-apple-darwin20.4.0 (64-bit) (2021/05/12)
  +   x <- rbeta(100, pi, 1/pi)
  + 
  +     set.seed(1234)
  +   fitdist(x, "beta")
  +   
  +   #binding example 
  +   fitdist(x, "beta", custom.optim=myoptimize, ui=rbind(1,1), ci=c(1/2,1/2), start=list(shape1=5, shape2=5))
  +   fitdist(x, "beta", lower= c(1/2,1/2), optim.method="L-BFGS-B")
  + }
  > 
  > 
  > #true example
  > library(GeneralizedHyperbolic)
  > args(dnig)
  function (x, mu = 0, delta = 1, alpha = 1, beta = 0, param = c(mu, 
      delta, alpha, beta)) 
  NULL
  > 
  > 
  > x <- rnig(100, 3, 1/2, 1/2, 1/4)
  > ui<-rbind(c(0,1,0,0),c(0,0,1,0),c(0,0,1,-1),c(0,0,1,1))
  > ci<-c(0,0,0,0)
  > fitdist(x, "nig", custom.optim=myoptimize, ui=ui, ci=ci, start=list(mu = 0, delta = 1, alpha = 1, beta = 0))
  Fitting of the distribution ' nig ' by maximum likelihood 
  Parameters:
         estimate Std. Error
  mu    3.0965290         NA
  delta 0.5761845         NA
  alpha 0.5773815         NA
  beta  0.1584031         NA
  Warning messages:
  1: In fitdist(x, "nig", custom.optim = myoptimize, ui = ui, ci = ci,  :
    The dnig function should return a vector of with NaN values when input has inconsistent parameters and not raise an error
  2: In fitdist(x, "nig", custom.optim = myoptimize, ui = ui, ci = ci,  :
    The pnig function should return a vector of with NaN values when input has inconsistent values and not raise an error
  > 
  > 
  > proc.time()
     user  system elapsed 
     1.14    0.20    1.32 
Running the tests in 'tests/t-fitdistcens.R' failed.
Complete output:
  > library(fitdistrplus)
  Loading required package: MASS
  Loading required package: survival
  > 
  > nsample <- 500
  > nsample <- 10
  > visualize <- FALSE # TRUE for manual tests with visualization of results
  > set.seed(1234)
  > 
  > # (6) custom optimisation function - example with the genetic algorithm
  > #
  > data(fluazinam)
  > log10EC50 <-log10(fluazinam)
  > 
  >   #wrap genoud function rgenoud package
  >   mygenoud  <-  function(fn, par, ...) 
  +   {
  +     require(rgenoud)
  +     res  <-  genoud(fn, starting.values=par, ...)        
  +     standardres  <-  c(res, convergence=0)
  +     
  +     return(standardres)
  +   }
  >   
  >   # call fitdistcens with a 'custom' optimization function
  >   fit.with.genoud <- fitdistcens(log10EC50, "logis", custom.optim=mygenoud, nvars=2, 
  +                                  start=list(location=0, scale=1),    
  +                                  Domains=cbind(c(0,0), c(5, 5)), boundary.enforcement=1, 
  +                                  print.level=1, hessian=TRUE)
  Loading required package: rgenoud
  ##  rgenoud (Version 5.9-0.10, Build Date: 2023-12-13)
  ##  See http://sekhon.berkeley.edu/rgenoud for additional documentation.
  ##  Please cite software as:
  ##   Walter Mebane, Jr. and Jasjeet S. Sekhon. 2011.
  ##   ``Genetic Optimization Using Derivatives: The rgenoud package for R.''
  ##   Journal of Statistical Software, 42(11): 1-26. 
  ##
  
  
  
  Sat Mar 23 14:28:33 2024
  Domains:
   0.000000e+00   <=  X1   <=    5.000000e+00 
   0.000000e+00   <=  X2   <=    5.000000e+00 
  
  Data Type: Floating Point
  Operators (code number, name, population) 
  	(1) Cloning........................... 	122
  	(2) Uniform Mutation.................. 	125
  	(3) Boundary Mutation................. 	125
  	(4) Non-Uniform Mutation.............. 	125
  	(5) Polytope Crossover................ 	125
  	(6) Simple Crossover.................. 	126
  	(7) Whole Non-Uniform Mutation........ 	125
  	(8) Heuristic Crossover............... 	126
  	(9) Local-Minimum Crossover........... 	0
  
  HARD Maximum Number of Generations: 100
  Maximum Nonchanging Generations: 10
  Population size       : 1000
  Convergence Tolerance: 1.000000e-03
  
  Using the BFGS Derivative Based Optimizer on the Best Individual Each Generation.
  Checking Gradients before Stopping.
  Not Using Out of Bounds Individuals But Allowing Trespassing.
  
  Minimization Problem.
  
  
  Generation#	    Solution Value
  
        0 	2.060254e+01
        1 	2.055391e+01
  
  'wait.generations' limit reached.
  No significant improvement in 10 generations.
  
  Solution Fitness Value: 2.055391e+01
  
  Parameters at the Solution (parameter, gradient):
  
   X[ 1] :	2.151909e+00	G[ 1] :	-9.910020e-06
   X[ 2] :	6.909662e-01	G[ 2] :	3.172033e-06
  
  Solution Found Generation 1
  Number of Generations Run 12
  
  Sat Mar 23 14:28:34 2024
  Total run time : 0 hours 0 minutes and 1 seconds
  >   
  >   summary(fit.with.genoud)
  Fitting of the distribution ' logis ' By maximum likelihood on censored data 
  Parameters
            estimate Std. Error
  location 2.1519094  0.3222566
  scale    0.6909662  0.1744835
  Loglikelihood:  -20.55391   AIC:  45.10781   BIC:  46.38593 
  Correlation matrix:
             [,1]       [,2]
  [1,] 1.00000000 0.05106398
  [2,] 0.05106398 1.00000000
  
  > 
  > 
  > # (9) check keepdata
  > #
  > if (visualize) # LONG TO RUN ON CRAN AND NEEDS VISALIZATION OF RESULTS
  + {
  +   set.seed(1234)
  +   x <- rexp(1e3, 5)
  +   # x <- data.frame(left=x, right=x+rexp(x, 1/2))
  +   x <- data.frame(left=x, right=x)
  +   f1 <- fitdistcens(x, "exp", keepdata=FALSE)
  +   f2 <- fitdistcens(x, "exp", keepdata=TRUE)
  +   f1$censdata
  +   f2$censdata
  +   plot(f1)
  +   plot(f2)
  +   
  + }
  > 
  > 
  > # (9) fixing parameters
  > #
  > x <- rexp(nsample, 5)
  > x <- data.frame(left=x, right=x+.1)
  > 
  > f1 <- fitdistcens(x, "gamma", fix.arg=list(shape=1.5))
  > f1
  Fitting of the distribution ' gamma ' on censored data by maximum likelihood 
  Parameters:
       estimate
  rate 9.840874
  Fixed parameters:
        value
  shape   1.5
  > f1$fix.arg
  $shape
  [1] 1.5
  
  > 
  > f1 <- fitdistcens(x, "gamma", fix.arg=function(x) list(shape=1.5))
  > f1
  Fitting of the distribution ' gamma ' on censored data by maximum likelihood 
  Parameters:
       estimate
  rate 9.840874
  Fixed parameters (computed by a user-supplied function):
        value
  shape   1.5
  > f1$fix.arg.fun
  function (x) 
  list(shape = 1.5)
  > 
  > 
  > # (10) weights
  > #
  > data(salinity)
  > salinity.unique <- unique(salinity)
  > string.unique <- paste(salinity.unique$left, salinity.unique$right)
  > string.salinity <- paste(salinity$left, salinity$right)
  > nobs <- nrow(salinity.unique)
  > salinity.weights <- numeric(nobs)
  > for (i in 1:nobs)
  + {
  +   salinity.weights[i] <- length(which(string.salinity == string.unique[i]))
  + }
  > cbind(salinity.unique, salinity.weights)
       left right salinity.weights
  1   20.00    NA                8
  6   21.50  21.5                1
  7   15.00  30.0                1
  8   20.00  25.0                1
  9   23.70  23.7                1
  10  25.00    NA                5
  13  20.00  30.0                3
  18   3.20  47.0                2
  20  26.10  26.1                1
  21  26.20  26.2                1
  22  25.00  30.0                1
  23  29.10  29.1                2
  25  30.00    NA               10
  30  30.00  30.0                1
  36  30.00  35.0                2
  38  35.00  35.0                1
  39  35.00    NA                2
  41  25.00  47.0                1
  42  30.00  47.0                1
  43  30.00  50.0                1
  44  35.00  47.0                1
  45  43.90  43.9                1
  46  35.00  55.0                1
  47  45.70  45.7                1
  48  47.00    NA                3
  50  47.00  47.0                2
  53  49.00  49.0                1
  54   0.10    NA                5
  59   0.23    NA                1
  60   3.20    NA                3
  63   6.40    NA                6
  69   3.20  12.8                1
  70   6.40  12.8                1
  71   0.10  20.0                1
  72   3.20  20.0                1
  73  12.60  12.6                1
  74  12.80  12.8                2
  75  12.80    NA               10
  86  12.80  15.0                2
  88  15.00    NA                7
  90  15.00  15.0                3
  98   6.40  25.0                1
  99  15.00  20.0                1
  100  6.40  30.0                1
  101 12.80  25.0                2
  103  3.20  35.0                1
  104 15.00  25.0                1
  105 12.80  20.0                1
  > 
  > (fa <- fitdistcens(salinity, "lnorm"))
  Fitting of the distribution ' lnorm ' on censored data by maximum likelihood 
  Parameters:
           estimate
  meanlog 3.3854230
  sdlog   0.4961333
  > (fb <- fitdistcens(salinity.unique, "lnorm", weights = salinity.weights)) # should give the same results
  Fitting of the distribution ' lnorm ' on censored data by maximum likelihood 
  Parameters:
           estimate
  meanlog 3.3854550
  sdlog   0.4961777
  Warning message:
  In mledist(censdata, distname, start = arg_startfix$start.arg, fix.arg = arg_startfix$fix.arg,  :
    weights are not taken into account in the default initial values
  > 
  > # (11) check the warning messages when using weights in the fit followed by functions
  > # that do not yet take weights into account
  > # with an example to be used later to see if weights are well taken into account
  > #
  > x <- rexp(100, 5)
  > x <- sort(x)
  > x <- data.frame(left=x, right=x+.1)
  > (f <- fitdistcens(x, "gamma", weights=c(rep(10, 50), rep(1, 50))))
  Fitting of the distribution ' gamma ' on censored data by maximum likelihood 
  Parameters:
         estimate
  shape  2.413252
  rate  19.147661
  Warning message:
  In mledist(censdata, distname, start = arg_startfix$start.arg, fix.arg = arg_startfix$fix.arg,  :
    weights are not taken into account in the default initial values
  > try(plot(f))
  Error in plot.fitdistcens(f) : 
    The plot of the fit is not yet available when using weights
  > try(cdfcompcens(f))
  Error in cdfcompcens(f) : 
    cdfcompcens is not yet available when using weights
  > (f2 <- fitdistcens(x, "weibull", weights=c(rep(10, 50), rep(1, 50))))
  Fitting of the distribution ' weibull ' on censored data by maximum likelihood 
  Parameters:
         estimate
  shape 1.3770720
  scale 0.1372821
  Warning message:
  In mledist(censdata, distname, start = arg_startfix$start.arg, fix.arg = arg_startfix$fix.arg,  :
    weights are not taken into account in the default initial values
  > try(cdfcompcens(list(f, f2)))
  Error in cdfcompcens(list(f, f2)) : 
    cdfcompcens is not yet available when using weights
  > try(bootdistcens(f))
  Error in bootdistcens(f) : 
    Bootstrap is not yet available when using weights
  > 
  > proc.time()
     user  system elapsed 
     2.98    0.29    3.21 
Running the tests in 'tests/t-gen-max-spacing-estim.R' failed.
Complete output:
  > library(fitdistrplus)
  Loading required package: MASS
  Loading required package: survival
  > 
  > set.seed(123)
  > 
  > nsample <- 10
  > 
  > #--------------------------------------------------------
  > # exponential sample
  > 
  > x1 <- rexp(nsample)
  > 
  > mseKL_exp <- msedist(x1, "exp")
  > mseJ_exp <- msedist(x1, "exp", phidiv = "J")
  > mseR2_exp <- msedist(x1, "exp", phidiv = "R", power.phidiv=2)
  > mseR1o2_exp <- msedist(x1, "exp", phidiv = "R", power.phidiv=1/2)
  > mseH3o2_exp <- msedist(x1, "exp", phidiv = "H", power.phidiv=3/2)
  > mseV3o2_exp <- msedist(x1, "exp", phidiv = "V", power.phidiv=3/2)
  > 
  > c(true=1, mseKL_exp$estimate, mseJ_exp$estimate, mseR2_exp$estimate, 
  +   mseR1o2_exp$estimate, mseH3o2_exp$estimate, mseV3o2_exp$estimate)
      true     rate     rate     rate     rate     rate     rate 
  1.000000 1.336924 1.329121 1.164714 1.195323 1.167898 1.171918 
  > 
  > 
  > 
  > mseKL_exp <- fitdist(x1, "exp", method="mse", phidiv="KL")
  > mseJ_exp <- fitdist(x1, "exp", method="mse", phidiv = "J")
  > mseR2_exp <- fitdist(x1, "exp", method="mse", phidiv = "R", power.phidiv=2)
  > mseR1o2_exp <- fitdist(x1, "exp", method="mse", phidiv = "R", power.phidiv=1/2)
  > mseH3o2_exp <- fitdist(x1, "exp", method="mse", phidiv = "H", power.phidiv=3/2)
  > mseV3o2_exp <- fitdist(x1, "exp", method="mse", phidiv = "V", power.phidiv=3/2)
  > 
  > 
  > 
  > gofstat(list(mseKL_exp, mseJ_exp, mseR2_exp, mseR1o2_exp, mseH3o2_exp, mseV3o2_exp))
  Goodness-of-fit statistics
                               1-mse-exp 2-mse-exp 3-mse-exp 4-mse-exp 5-mse-exp
  Kolmogorov-Smirnov statistic 0.2549881 0.2566077 0.2916788 0.2850105 0.2909822
  Cramer-von Mises statistic   0.1544480 0.1568271 0.2185483 0.2052244 0.2171190
  Anderson-Darling statistic   1.0571407 1.0703179 1.4200650 1.3435948 1.4118434
                               6-mse-exp
  Kolmogorov-Smirnov statistic 0.2901037
  Cramer-von Mises statistic   0.2153293
  Anderson-Darling statistic   1.4015539
  
  Goodness-of-fit criteria
                                 1-mse-exp 2-mse-exp 3-mse-exp 4-mse-exp
  Akaike's Information Criterion  13.22043  13.23811  13.78497  13.65601
  Bayesian Information Criterion  13.52301  13.54070  14.08756  13.95860
                                 5-mse-exp 6-mse-exp
  Akaike's Information Criterion  13.77093  13.75341
  Bayesian Information Criterion  14.07351  14.05599
  > 
  > 
  > cdfcomp(list(mseKL_exp, mseJ_exp, mseR2_exp, mseR1o2_exp, mseH3o2_exp, mseV3o2_exp), do.points=FALSE,
  +         legendtext = c("Kullback-Leibler", "Jeffreys", "Renyi alpha=2", "Renyi alpha=1/2", "Hellinger p=3/2",
  +                        "Vajda beta=3/2"))
  > 
  > qqcomp(list(mseKL_exp, mseJ_exp, mseR2_exp, mseR1o2_exp, mseH3o2_exp, mseV3o2_exp), 
  +         legendtext = c("Kullback-Leibler", "Jeffreys", "Renyi alpha=2", "Renyi alpha=1/2", "Hellinger p=3/2",
  +                        "Vajda beta=3/2"))
  > 
  > denscomp(list(mseKL_exp, mseJ_exp, mseR2_exp, mseR1o2_exp, mseH3o2_exp, mseV3o2_exp), demp = TRUE,
  +        legendtext = c("Kullback-Leibler", "Jeffreys", "Renyi alpha=2", "Renyi alpha=1/2", "Hellinger p=3/2",
  +                       "Vajda beta=3/2"))
  > 
  > 
  > #defensive test 
  > 
  > try(msedist(x1, "exp", phidiv="ABC"))
  Error in match.arg(phidiv, c("KL", "J", "R", "H", "V")) : 
    'arg' should be one of "KL", "J", "R", "H", "V"
  > try(msedist(x1, "exp", phidiv="K", power.phidiv = "a"))
  Error in msedist(x1, "exp", phidiv = "K", power.phidiv = "a") : 
    is.null(power.phidiv) is not TRUE
  > try(msedist(x1, "exp", phidiv="J", power.phidiv = "a"))
  Error in msedist(x1, "exp", phidiv = "J", power.phidiv = "a") : 
    is.null(power.phidiv) is not TRUE
  > try(msedist(x1, "exp", phidiv="R", power.phidiv = 0))
  Error in msedist(x1, "exp", phidiv = "R", power.phidiv = 0) : 
    power.phidiv > 0 && power.phidiv != 1 is not TRUE
  > try(msedist(x1, "exp", phidiv="R", power.phidiv = 1:10))
  Error in msedist(x1, "exp", phidiv = "R", power.phidiv = 1:10) : 
    length(power.phidiv) == 1 is not TRUE
  > try(msedist(x1, "exp", phidiv="H", power.phidiv = 0))
  Error in msedist(x1, "exp", phidiv = "H", power.phidiv = 0) : 
    power.phidiv >= 1 is not TRUE
  > try(msedist(x1, "exp", phidiv="H", power.phidiv = 1:10))
  Error in msedist(x1, "exp", phidiv = "H", power.phidiv = 1:10) : 
    length(power.phidiv) == 1 is not TRUE
  > 
  > 
  > 
  > #--------------------------------------------------------
  > # Poisson sample
  > 
  > x1 <- rpois(nsample, lambda=15)
  > 
  > #no weight
  > mseKL_pois1 <- fitdist(x1, "pois", method="mse", phidiv="KL")
  > mseJ_pois1 <- fitdist(x1, "pois", method="mse", phidiv = "J")
  > mseR2_pois1 <- fitdist(x1, "pois", method="mse", phidiv = "R", power.phidiv=2)
  > mseR1o2_pois1 <- fitdist(x1, "pois", method="mse", phidiv = "R", power.phidiv=1/2)
  > mseH3o2_pois1 <- fitdist(x1, "pois", method="mse", phidiv = "H", power.phidiv=3/2)
  > mseV3o2_pois1 <- fitdist(x1, "pois",method="mse",  phidiv = "V", power.phidiv=3/2)
  > 
  > 
  > #with weight
  > mseKL_pois2 <- fitdist(unique(sort(x1)), "pois", method="mse", phidiv="KL", weights=as.numeric(table(x1)))
  Warning message:
  In msedist(data, distname, start = arg_startfix$start.arg, fix.arg = arg_startfix$fix.arg,  :
    weights are not taken into account in the default initial values
  > mseJ_pois2 <- fitdist(unique(sort(x1)), "pois", method="mse", phidiv = "J", weights=as.numeric(table(x1)))
  Warning message:
  In msedist(data, distname, start = arg_startfix$start.arg, fix.arg = arg_startfix$fix.arg,  :
    weights are not taken into account in the default initial values
  > mseR2_pois2 <- fitdist(unique(sort(x1)), "pois", method="mse", phidiv = "R", power.phidiv=2, weights=as.numeric(table(x1)))
  Warning message:
  In msedist(data, distname, start = arg_startfix$start.arg, fix.arg = arg_startfix$fix.arg,  :
    weights are not taken into account in the default initial values
  > mseR1o2_pois2 <- fitdist(unique(sort(x1)), "pois", method="mse", phidiv = "R", power.phidiv=1/2, weights=as.numeric(table(x1)))
  Warning message:
  In msedist(data, distname, start = arg_startfix$start.arg, fix.arg = arg_startfix$fix.arg,  :
    weights are not taken into account in the default initial values
  > mseH3o2_pois2 <- fitdist(unique(sort(x1)), "pois", method="mse", phidiv = "H", power.phidiv=3/2, weights=as.numeric(table(x1)))
  Warning message:
  In msedist(data, distname, start = arg_startfix$start.arg, fix.arg = arg_startfix$fix.arg,  :
    weights are not taken into account in the default initial values
  > mseV3o2_pois2 <- fitdist(unique(sort(x1)), "pois",method="mse",  phidiv = "V", power.phidiv=3/2, weights=as.numeric(table(x1)))
  Warning message:
  In msedist(data, distname, start = arg_startfix$start.arg, fix.arg = arg_startfix$fix.arg,  :
    weights are not taken into account in the default initial values
  > 
  > #taking into account partially unbias the estimation of the true cdf
  > cdfcomp(list(mseKL_pois1, mseJ_pois1, mseR2_pois1), do.points=FALSE, fitlty=1,
  +         legendtext = c("Kullback-Leibler", "Jeffreys", "Renyi alpha=2"))
  > 
  > cdfcomp(list(mseKL_pois2, mseJ_pois2, mseR2_pois2), do.points=FALSE, add=TRUE, 
  +         fitlty=2, addlegend = FALSE, datacol="white")
  > 
  > 
  > cdfcomp(list(mseR1o2_pois1, mseH3o2_pois1, mseV3o2_pois1), do.points=FALSE, fitlty=1,
  +         legendtext = c("Renyi alpha=1/2", "Hellinger p=3/2", "Vajda beta=3/2"))
  > 
  > cdfcomp(list(mseR1o2_pois2, mseH3o2_pois2, mseV3o2_pois2), do.points=FALSE, add=TRUE, 
  +         fitlty=2, addlegend = FALSE, datacol="white")
  > 
  > 
  > 
  > proc.time()
     user  system elapsed 
     1.48    0.23    1.70 
Running the tests in 'tests/t-llplot.R' failed.
Complete output:
  > require(fitdistrplus)
  Loading required package: fitdistrplus
  Loading required package: MASS
  Loading required package: survival
  > visualize <- FALSE # TRUE for manual tests with visualization of results
  > nsample <- 10000
  > nsample <- 10
  > 
  > # (1) tests with the Burr distribution (three parameters)
  > #
  > if(any(installed.packages()[, "Package"] == "actuar"))
  + {
  +   require(actuar)
  +   
  +   data(endosulfan)
  +   ATV <-endosulfan$ATV
  +   library("actuar")
  +   fBurr <- fitdist(ATV, "burr", start = list(shape1 = 0.3, shape2 = 1, rate = 1))
  +   llplot(fBurr)
  +   
  +   fBurr2 <- fitdist(ATV, "burr", start = list(shape1 = 0.3, shape2 = 1),
  +                     fix.arg = list(rate = 1.5))
  +   llplot(fBurr2)
  +   
  +   fBurr3 <- fitdist(ATV, "burr", start = list(shape1 = 0.3, rate = 1),
  +                     fix.arg = list(shape2 = 1.5))
  +   llplot(fBurr3)
  + }
  Loading required package: actuar
  
  Attaching package: 'actuar'
  
  The following objects are masked from 'package:stats':
  
      sd, var
  
  The following object is masked from 'package:grDevices':
  
      cm
  
  > 
  > 
  > # (2) An example on discrete data with or without weights
  > #
  > set.seed(1234)
  > x <- rpois(nsample, 10)
  > xtab <- table(x)
  > xval <- sort(unique(x))
  > f1 <- fitdist(x, "pois")
  > f2 <- fitdist(xval, "pois", weights = xtab)
  Warning message:
  In mledist(data, distname, start = arg_startfix$start.arg, fix.arg = arg_startfix$fix.arg,  :
    weights are not taken into account in the default initial values
  > 
  > f1$estimate
  lambda 
     8.4 
  > f2$estimate # should give the same
  lambda 
     8.4 
  > llplot(f1, fit.show = TRUE)
  > llplot(f2, fit.show = TRUE) # should give the same
  > llplot(f1, loglik = FALSE, fit.show = TRUE)
  > llplot(f2, loglik = FALSE,fit.show = TRUE) # should give the same
  > 
  > # (3) An example on censored data with or without weights
  > #
  > if(visualize)
  + {
  +   data(salinity)
  +   salinity.unique <- unique(salinity)
  +   string.unique <- paste(salinity.unique$left, salinity.unique$right)
  +   string.salinity <- paste(salinity$left, salinity$right)
  +   nobs <- nrow(salinity.unique)
  +   salinity.weights <- numeric(nobs)
  +   for (i in 1:nobs)
  +   {
  +     salinity.weights[i] <- length(which(string.salinity == string.unique[i]))
  +   }
  +   cbind(salinity.unique, salinity.weights)
  +   
  +   (fa <- fitdistcens(salinity, "lnorm"))
  +   (fb <- fitdistcens(salinity.unique, "lnorm", weights = salinity.weights))
  +   llplot(fa, fit.show = TRUE)
  +   llplot(fb, fit.show = TRUE) # should give the same
  +   llplot(fa, fit.show = TRUE, loglik = FALSE)
  +   llplot(fb, fit.show = TRUE, loglik = FALSE) # should give the same
  +   
  + }
  > 
  > proc.time()
     user  system elapsed 
     2.78    0.25    3.03 
Running the tests in 'tests/t-logLik-vcov-coef.R' failed.
Complete output:
  > library(fitdistrplus)
  Loading required package: MASS
  Loading required package: survival
  > 
  > # (1) basic fit of a gamma distribution by maximum likelihood estimation
  > #
  > data(groundbeef)
  > serving <- groundbeef$serving
  > fitg <- fitdist(serving, "gamma")
  > 
  > logLik(fitg)
  [1] -1253.625
  > vcov(fitg)
              shape         rate
  shape 0.116510296 1.580943e-03
  rate  0.001580943 2.435884e-05
  > coef(fitg)
       shape       rate 
  4.00825257 0.05441911 
  > 
  > fitg <- fitdist(serving, "gamma", method="mme")
  > 
  > logLik(fitg)
  [1] -1253.825
  > vcov(fitg) 
  NULL
  Warning message:
  In vcov.fitdist(fitg) :
    The variance-covariance matrix can only be calculated for fits using the mle method
  > coef(fitg)
       shape       rate 
  4.22848617 0.05741663 
  > 
  > # (2) Fit of a lognormal distribution to bacterial contamination data
  > #
  > data(smokedfish)
  > fitsf  <-  fitdistcens(smokedfish,"lnorm")
  > logLik(fitsf)
  [1] -90.65154
  > vcov(fitsf) 
              meanlog       sdlog
  meanlog  0.21502901 -0.09782285
  sdlog   -0.09782285  0.23781325
  > coef(fitsf)
    meanlog     sdlog 
  -3.627606  3.544570 
  > 
  > 
  > proc.time()
     user  system elapsed 
     1.37    0.15    1.51 
Running the tests in 'tests/t-mgedist.R' failed.
Complete output:
  > library(fitdistrplus)
  Loading required package: MASS
  Loading required package: survival
  > 
  > nsample <- 10
  > 
  > # (1) Fit of a Weibull distribution to serving size data by maximum 
  > # goodness-of-fit estimation using all the distances available
  > # 
  > 
  > data(groundbeef)
  > serving <- groundbeef$serving
  > mgedist(serving, "weibull", gof="CvM")
  $estimate
      shape     scale 
   2.093204 82.660014 
  
  $convergence
  [1] 0
  
  $value
  [1] 0.6556672
  
  $hessian
             shape      scale
  shape 4.05295367 0.09244476
  scale 0.09244476 0.02418777
  
  $optim.function
  [1] "optim"
  
  $optim.method
  [1] "Nelder-Mead"
  
  $fix.arg
  NULL
  
  $fix.arg.fun
  NULL
  
  $weights
  NULL
  
  $counts
  function gradient 
        65       NA 
  
  $optim.message
  NULL
  
  $loglik
  [1] -1255.623
  
  $gof
  [1] "CvM"
  
  > mgedist(serving, "weibull", gof="CvM", silent = FALSE)
  $estimate
      shape     scale 
   2.093204 82.660014 
  
  $convergence
  [1] 0
  
  $value
  [1] 0.6556672
  
  $hessian
             shape      scale
  shape 4.05295367 0.09244476
  scale 0.09244476 0.02418777
  
  $optim.function
  [1] "optim"
  
  $optim.method
  [1] "Nelder-Mead"
  
  $fix.arg
  NULL
  
  $fix.arg.fun
  NULL
  
  $weights
  NULL
  
  $counts
  function gradient 
        65       NA 
  
  $optim.message
  NULL
  
  $loglik
  [1] -1255.623
  
  $gof
  [1] "CvM"
  
  Warning messages:
  1: In pweibull(c(10, 11.5, 17, 20, 20, 20, 20, 20, 20, 20, 20, 24,  :
    NaNs produced
  2: In pweibull(c(10, 11.5, 17, 20, 20, 20, 20, 20, 20, 20, 20, 24,  :
    NaNs produced
  > mgedist(serving, "weibull", gof="KS")
  $estimate
      shape     scale 
   2.065634 81.450487 
  
  $convergence
  [1] 0
  
  $value
  [1] 0.112861
  
  $hessian
             shape    scale
  shape 122.668263 6.509057
  scale   6.509057 7.599584
  
  $optim.function
  [1] "optim"
  
  $optim.method
  [1] "Nelder-Mead"
  
  $fix.arg
  NULL
  
  $fix.arg.fun
  NULL
  
  $weights
  NULL
  
  $counts
  function gradient 
       127       NA 
  
  $optim.message
  NULL
  
  $loglik
  [1] -1255.975
  
  $gof
  [1] "KS"
  
  > mgedist(serving, "weibull", gof="AD")
  $estimate
      shape     scale 
   2.125473 82.890260 
  
  $convergence
  [1] 0
  
  $value
  [1] 3.501035
  
  $hessian
             shape     scale
  shape 29.4165108 0.1823375
  scale  0.1823375 0.1354409
  
  $optim.function
  [1] "optim"
  
  $optim.method
  [1] "Nelder-Mead"
  
  $fix.arg
  NULL
  
  $fix.arg.fun
  NULL
  
  $weights
  NULL
  
  $counts
  function gradient 
        63       NA 
  
  $optim.message
  NULL
  
  $loglik
  [1] -1255.392
  
  $gof
  [1] "AD"
  
  > mgedist(serving, "weibull", gof="ADR")
  $estimate
      shape     scale 
   2.072087 82.761868 
  
  $convergence
  [1] 0
  
  $value
  [1] 1.610479
  
  $hessian
             shape       scale
  shape 13.5240921 -0.33242262
  scale -0.3324226  0.07977375
  
  $optim.function
  [1] "optim"
  
  $optim.method
  [1] "Nelder-Mead"
  
  $fix.arg
  NULL
  
  $fix.arg.fun
  NULL
  
  $weights
  NULL
  
  $counts
  function gradient 
        65       NA 
  
  $optim.message
  NULL
  
  $loglik
  [1] -1255.836
  
  $gof
  [1] "ADR"
  
  > mgedist(serving, "weibull", gof="ADL")
  $estimate
      shape     scale 
   2.197498 82.016005 
  
  $convergence
  [1] 0
  
  $value
  [1] 1.845939
  
  $hessian
             shape      scale
  shape 15.3272022 0.54407117
  scale  0.5440712 0.05549883
  
  $optim.function
  [1] "optim"
  
  $optim.method
  [1] "Nelder-Mead"
  
  $fix.arg
  NULL
  
  $fix.arg.fun
  NULL
  
  $weights
  NULL
  
  $counts
  function gradient 
        65       NA 
  
  $optim.message
  NULL
  
  $loglik
  [1] -1255.415
  
  $gof
  [1] "ADL"
  
  > mgedist(serving, "weibull", gof="AD2R")
  $estimate
     shape    scale 
   1.90328 81.33464 
  
  $convergence
  [1] 0
  
  $value
  [1] 11.56415
  
  $hessian
            shape       scale
  shape 334.61081 -10.4227495
  scale -10.42275   0.5223167
  
  $optim.function
  [1] "optim"
  
  $optim.method
  [1] "Nelder-Mead"
  
  $fix.arg
  NULL
  
  $fix.arg.fun
  NULL
  
  $weights
  NULL
  
  $counts
  function gradient 
        67       NA 
  
  $optim.message
  NULL
  
  $loglik
  [1] -1259.112
  
  $gof
  [1] "AD2R"
  
  > mgedist(serving, "weibull", gof="AD2L")
  $estimate
      shape     scale 
   2.483836 78.252113 
  
  $convergence
  [1] 0
  
  $value
  [1] 9.786977
  
  $hessian
             shape     scale
  shape 113.511932 4.1108355
  scale   4.110836 0.2341312
  
  $optim.function
  [1] "optim"
  
  $optim.method
  [1] "Nelder-Mead"
  
  $fix.arg
  NULL
  
  $fix.arg.fun
  NULL
  
  $weights
  NULL
  
  $counts
  function gradient 
        69       NA 
  
  $optim.message
  NULL
  
  $loglik
  [1] -1265.933
  
  $gof
  [1] "AD2L"
  
  > mgedist(serving, "weibull", gof="AD2")
  $estimate
      shape     scale 
   2.081168 85.281194 
  
  $convergence
  [1] 0
  
  $value
  [1] 26.95166
  
  $hessian
           shape       scale
  shape 534.9606 -10.5940982
  scale -10.5941   0.7606462
  
  $optim.function
  [1] "optim"
  
  $optim.method
  [1] "Nelder-Mead"
  
  $fix.arg
  NULL
  
  $fix.arg.fun
  NULL
  
  $weights
  NULL
  
  $counts
  function gradient 
        67       NA 
  
  $optim.message
  NULL
  
  $loglik
  [1] -1256.313
  
  $gof
  [1] "AD2"
  
  > 
  > # (2) Fit of a uniform distribution using Cramer-von Mises or
  > # Kolmogorov-Smirnov distance
  > # 
  > 
  > set.seed(1234)
  > u <- runif(nsample, min=5, max=10)
  > mgedist(u, "unif", gof="CvM")
  $estimate
  min max 
    0   1 
  
  $convergence
  [1] 0
  
  $value
  [1] 3.333333
  
  $hessian
      min max
  min   0   0
  max   0   0
  
  $optim.function
  [1] "optim"
  
  $optim.method
  [1] "Nelder-Mead"
  
  $fix.arg
  NULL
  
  $fix.arg.fun
  NULL
  
  $weights
  NULL
  
  $counts
  function gradient 
         3       NA 
  
  $optim.message
  NULL
  
  $loglik
  [1] -Inf
  
  $gof
  [1] "CvM"
  
  > mgedist(u, "unif", gof="KS")
  $estimate
  min max 
    0   1 
  
  $convergence
  [1] 0
  
  $value
  [1] 1
  
  $hessian
      min max
  min   0   0
  max   0   0
  
  $optim.function
  [1] "optim"
  
  $optim.method
  [1] "Nelder-Mead"
  
  $fix.arg
  NULL
  
  $fix.arg.fun
  NULL
  
  $weights
  NULL
  
  $counts
  function gradient 
         3       NA 
  
  $optim.message
  NULL
  
  $loglik
  [1] -Inf
  
  $gof
  [1] "KS"
  
  > 
  > # (3) Fit of a triangular distribution using Cramer-von Mises or
  > # Kolmogorov-Smirnov distance
  > # 
  > require(mc2d)
  Loading required package: mc2d
  Loading required package: mvtnorm
  
  Attaching package: 'mc2d'
  
  The following objects are masked from 'package:base':
  
      pmax, pmin
  
  > set.seed(1234)
  > t <- rtriang(nsample, min=5, mode=6, max=10)
  > mgedist(t, "triang", start = list(min=4, mode=6,max=9), gof="CvM")
  $estimate
       min     mode      max 
  4.602048 7.367066 8.682229 
  
  $convergence
  [1] 0
  
  $value
  [1] 0.08579645
  
  $hessian
             min      mode       max
  min  0.1914605 0.2739648 0.1892873
  mode 0.2739648 0.7323464 0.4060730
  max  0.1892873 0.4060730 0.2910040
  
  $optim.function
  [1] "optim"
  
  $optim.method
  [1] "Nelder-Mead"
  
  $fix.arg
  NULL
  
  $fix.arg.fun
  NULL
  
  $weights
  NULL
  
  $counts
  function gradient 
       180       NA 
  
  $optim.message
  NULL
  
  $loglik
  [1] -11.84833
  
  $gof
  [1] "CvM"
  
  Warning message:
  In checkparamlist(arg_startfix$start.arg, arg_startfix$fix.arg,  :
    Some parameter names have no starting/fixed value but have a default value: mean.
  > mgedist(t, "triang", start = list(min=4, mode=6,max=9), gof="KS")
  $estimate
       min     mode      max 
  5.191300 7.414283 8.306100 
  
  $convergence
  [1] 0
  
  $value
  [1] 0.1853719
  
  $hessian
             min     mode       max
  min  148.15115 110.9594  93.74732
  mode 110.95941 292.0984 208.67195
  max   93.74732 208.6720 208.66467
  
  $optim.function
  [1] "optim"
  
  $optim.method
  [1] "Nelder-Mead"
  
  $fix.arg
  NULL
  
  $fix.arg.fun
  NULL
  
  $weights
  NULL
  
  $counts
  function gradient 
       298       NA 
  
  $optim.message
  NULL
  
  $loglik
  [1] -Inf
  
  $gof
  [1] "KS"
  
  Warning message:
  In checkparamlist(arg_startfix$start.arg, arg_startfix$fix.arg,  :
    Some parameter names have no starting/fixed value but have a default value: mean.
  > 
  > 
  > # (4) scaling problem
  > # the simulated dataset (below) has particularly small values, hence without scaling (10^0),
  > # the optimization raises an error. The for loop shows how scaling by 10^i
  > # for i=1,...,6 makes the fitting procedure work correctly.
  > if(FALSE)
  + {
  +   set.seed(1234)
  +   x2 <- rnorm(nsample, 1e-4, 2e-4)
  +   for(i in 6:0)
  +     cat(i, try(mgedist(x2*10^i,"cauchy")$estimate, silent=TRUE), "\n")
  + }
  > 
  > # (5) scaling problem
  > #
  > if(FALSE)
  + {
  +   x <- c(-0.00707717, -0.000947418, -0.00189753, 
  +          -0.000474947, -0.00190205, -0.000476077, 0.00237812, 0.000949668, 
  +          0.000474496, 0.00284226, -0.000473149, -0.000473373, 0, 0, 0.00283688, 
  +          -0.0037843, -0.0047506, -0.00238379, -0.00286807, 0.000478583, 
  +          0.000478354, -0.00143575, 0.00143575, 0.00238835, 0.0042847, 
  +          0.00237248, -0.00142281, -0.00142484, 0, 0.00142484, 0.000948767, 
  +          0.00378609, -0.000472478, 0.000472478, -0.0014181, 0, -0.000946522, 
  +          -0.00284495, 0, 0.00331832, 0.00283554, 0.00141476, -0.00141476, 
  +          -0.00188947, 0.00141743, -0.00236351, 0.00236351, 0.00235794, 
  +          0.00235239, -0.000940292, -0.0014121, -0.00283019, 0.000472255, 
  +          0.000472032, 0.000471809, -0.0014161, 0.0014161, -0.000943842, 
  +          0.000472032, -0.000944287, -0.00094518, -0.00189304, -0.000473821, 
  +          -0.000474046, 0.00331361, -0.000472701, -0.000946074, 0.00141878, 
  +          -0.000945627, -0.00189394, -0.00189753, -0.0057143, -0.00143369, 
  +          -0.00383326, 0.00143919, 0.000479272, -0.00191847, -0.000480192, 
  +          0.000960154, 0.000479731, 0, 0.000479501, 0.000958313, -0.00383878, 
  +          -0.00240674, 0.000963391, 0.000962464, -0.00192586, 0.000481812, 
  +          -0.00241138, -0.00144963)
  +   
  +   
  +   #only i == 0, no scaling, should not converge.
  +   for(i in 6:0)
  +     cat(i, try(mgedist(x*10^i,"cauchy")$estimate, silent=TRUE), "\n")
  + }
  > 
  > 
  > # (6) test error messages
  > #
  > 
  > dnorm2 <- pnorm2 <- function(x, a)
  +   "NA"
  > x <- rexp(10)
  > 
  > #should get a one-line error 
  > res <- mgedist(x, "norm2", start=list(a=1))
  <simpleError in theop - (2 * 1:n - 1)/(2 * n): non-numeric argument to binary operator>
  > #as in 
  > attr(try(log("a"), silent=TRUE), "condition")
  <simpleError in log("a"): non-numeric argument to mathematical function>
  > 
  > 
  > # (7) test the component optim.message
  > x <- rnorm(1000)
  > #change parameter to obtain unsuccessful convergence
  > mgedist(x, "norm", control=list(maxit=2), start=list(mean=1e5, sd=1), optim.method="L-BFGS-B", lower=0)
  $estimate
   mean    sd 
  1e+05 1e+00 
  
  $convergence
  [1] 0
  
  $value
  [1] 333.3333
  
  $hessian
       mean sd
  mean    0  0
  sd      0  0
  
  $optim.function
  [1] "optim"
  
  $optim.method
  [1] "L-BFGS-B"
  
  $fix.arg
  NULL
  
  $fix.arg.fun
  NULL
  
  $weights
  NULL
  
  $counts
  function gradient 
         1        1 
  
  $optim.message
  [1] "CONVERGENCE: NORM OF PROJECTED GRADIENT <= PGTOL"
  
  $loglik
  [1] -Inf
  
  $gof
  [1] "CvM"
  
  > 
  > # (8) test bounds
  > x <- rnorm(1000)
  > mgedist(x, "norm", optim.method="L-BFGS-B", lower=c(-Inf, 0)) #optim and L-BFGS-B
  $estimate
        mean         sd 
  0.01134416 0.97945793 
  
  $convergence
  [1] 0
  
  $value
  [1] 0.02003669
  
  $hessian
               mean          sd
  mean 191.57400735 -0.05237777
  sd    -0.05237777 62.19267021
  
  $optim.function
  [1] "optim"
  
  $optim.method
  [1] "L-BFGS-B"
  
  $fix.arg
  NULL
  
  $fix.arg.fun
  NULL
  
  $weights
  NULL
  
  $counts
  function gradient 
         8        8 
  
  $optim.message
  [1] "CONVERGENCE: REL_REDUCTION_OF_F <= FACTR*EPSMCH"
  
  $loglik
  [1] -1399.894
  
  $gof
  [1] "CvM"
  
  > mgedist(x, "norm", optim.method="Nelder", lower=c(-Inf, 0))
  $estimate
        mean         sd 
  0.01134315 0.97945779 
  
  $convergence
  [1] 0
  
  $value
  [1] 0.02003669
  
  $hessian
  NULL
  
  $optim.function
  [1] "constrOptim"
  
  $optim.method
  [1] "Nelder-Mead"
  
  $fix.arg
  NULL
  
  $fix.arg.fun
  NULL
  
  $weights
  NULL
  
  $counts
  function gradient 
        59       NA 
  
  $optim.message
  NULL
  
  $loglik
  [1] -1399.894
  
  $gof
  [1] "CvM"
  
  > 
  > 
  > 
  > proc.time()
     user  system elapsed 
     2.43    0.29    2.70 
Running the tests in 'tests/t-mledist-cens.R' failed.
Complete output:
  > library(fitdistrplus)
  Loading required package: MASS
  Loading required package: survival
  > 
  > 
  > 
  > 
  > # (1) fit a user-defined continuous distribution (Gumbel) to censored data.
  > #
  > 
  > data(fluazinam)
  > log10EC50 <-log10(fluazinam)
  > # definition of the Gumbel distribution
  > dgumbel  <-  function(x,a,b) 1/b*exp((a-x)/b)*exp(-exp((a-x)/b))
  > pgumbel  <-  function(q,a,b) exp(-exp((a-q)/b))
  > qgumbel  <-  function(p,a,b) a-b*log(-log(p))
  > 
  > mledist(log10EC50, "gumbel", start=list(a=0,b=2), optim.method="Nelder-Mead")
  $estimate
         a        b 
  1.632726 1.144113 
  
  $convergence
  [1] 0
  
  $value
  [1] 20.27771
  
  $hessian
            a         b
  a 10.683955 -4.764691
  b -4.764691 16.810611
  
  $optim.function
  [1] "optim"
  
  $optim.method
  [1] "Nelder-Mead"
  
  $fix.arg
  NULL
  
  $fix.arg.fun
  NULL
  
  $weights
  NULL
  
  $counts
  function gradient 
        65       NA 
  
  $optim.message
  NULL
  
  $loglik
  [1] -20.27771
  
  > mledist(log10EC50, "gumbel", start=list(a=0,b=2)) #default NM
  $estimate
         a        b 
  1.632726 1.144113 
  
  $convergence
  [1] 0
  
  $value
  [1] 20.27771
  
  $hessian
            a         b
  a 10.683955 -4.764691
  b -4.764691 16.810611
  
  $optim.function
  [1] "optim"
  
  $optim.method
  [1] "Nelder-Mead"
  
  $fix.arg
  NULL
  
  $fix.arg.fun
  NULL
  
  $weights
  NULL
  
  $counts
  function gradient 
        65       NA 
  
  $optim.message
  NULL
  
  $loglik
  [1] -20.27771
  
  > 
  > 
  > 
  > # (2) test optimization arguments to censored data MLE.
  > #
  > 
  > mledist(log10EC50, "lnorm", optim.method="BFGS")
  $estimate
    meanlog     sdlog 
  0.6229737 0.9281694 
  
  $convergence
  [1] 0
  
  $value
  [1] 21.06638
  
  $hessian
            meanlog     sdlog
  meanlog 15.358336 -2.676329
  sdlog   -2.676329 23.992338
  
  $optim.function
  [1] "optim"
  
  $optim.method
  [1] "BFGS"
  
  $fix.arg
  NULL
  
  $fix.arg.fun
  NULL
  
  $weights
  NULL
  
  $counts
  function gradient 
        13        6 
  
  $optim.message
  NULL
  
  $loglik
  [1] -21.06638
  
  > mledist(log10EC50, "lnorm", optim.method="Nelder")
  $estimate
    meanlog     sdlog 
  0.6230097 0.9282263 
  
  $convergence
  [1] 0
  
  $value
  [1] 21.06638
  
  $hessian
            meanlog     sdlog
  meanlog 15.356415 -2.676889
  sdlog   -2.676889 23.985490
  
  $optim.function
  [1] "optim"
  
  $optim.method
  [1] "Nelder-Mead"
  
  $fix.arg
  NULL
  
  $fix.arg.fun
  NULL
  
  $weights
  NULL
  
  $counts
  function gradient 
        45       NA 
  
  $optim.message
  NULL
  
  $loglik
  [1] -21.06638
  
  > 
  > #optim() is used
  > mledist(log10EC50, "lnorm", optim.method="L-BFGS-B", lower=0)
  $estimate
    meanlog     sdlog 
  0.6229875 0.9281720 
  
  $convergence
  [1] 0
  
  $value
  [1] 21.06638
  
  $hessian
            meanlog     sdlog
  meanlog 15.358242 -2.676762
  sdlog   -2.676762 23.992168
  
  $optim.function
  [1] "optim"
  
  $optim.method
  [1] "L-BFGS-B"
  
  $fix.arg
  NULL
  
  $fix.arg.fun
  NULL
  
  $weights
  NULL
  
  $counts
  function gradient 
         9        9 
  
  $optim.message
  [1] "CONVERGENCE: REL_REDUCTION_OF_F <= FACTR*EPSMCH"
  
  $loglik
  [1] -21.06638
  
  > mledist(log10EC50, "lnorm", optim.method="BFGS", lower=0) #a simple warning
  $estimate
    meanlog     sdlog 
  0.6229875 0.9281720 
  
  $convergence
  [1] 0
  
  $value
  [1] 21.06638
  
  $hessian
            meanlog     sdlog
  meanlog 15.358242 -2.676762
  sdlog   -2.676762 23.992168
  
  $optim.function
  [1] "optim"
  
  $optim.method
  [1] "L-BFGS-B"
  
  $fix.arg
  NULL
  
  $fix.arg.fun
  NULL
  
  $weights
  NULL
  
  $counts
  function gradient 
         9        9 
  
  $optim.message
  [1] "CONVERGENCE: REL_REDUCTION_OF_F <= FACTR*EPSMCH"
  
  $loglik
  [1] -21.06638
  
  Warning message:
  In mledist(log10EC50, "lnorm", optim.method = "BFGS", lower = 0) :
    The BFGS method cannot be used with bounds without provided the gradient. The method is changed to L-BFGS-B.
  > 
  > #constrOptim() is used
  > mledist(log10EC50, "lnorm", optim.method="Nelder", lower=0)
  $estimate
    meanlog     sdlog 
  0.6230097 0.9282263 
  
  $convergence
  [1] 0
  
  $value
  [1] 21.06638
  
  $hessian
  NULL
  
  $optim.function
  [1] "constrOptim"
  
  $optim.method
  [1] "Nelder-Mead"
  
  $fix.arg
  NULL
  
  $fix.arg.fun
  NULL
  
  $weights
  NULL
  
  $counts
  function gradient 
        45       NA 
  
  $optim.message
  NULL
  
  $loglik
  [1] -21.06638
  
  > mledist(log10EC50, "lnorm", lower=0) #error
  $estimate
    meanlog     sdlog 
  0.6230097 0.9282263 
  
  $convergence
  [1] 0
  
  $value
  [1] 21.06638
  
  $hessian
  NULL
  
  $optim.function
  [1] "constrOptim"
  
  $optim.method
  [1] "Nelder-Mead"
  
  $fix.arg
  NULL
  
  $fix.arg.fun
  NULL
  
  $weights
  NULL
  
  $counts
  function gradient 
        45       NA 
  
  $optim.message
  NULL
  
  $loglik
  [1] -21.06638
  
  > 
  > # (3) weighted MLE
  > #
  > xleft <- c(-1.8, -0.6, -0.1, 0.07, 0.14, 1, 1.2, 1.2, 1.2)
  > xright <- c(-1.8, -0.6, -0.1, 0.07, 0.14, 1, NA, NA, NA)
  > d <- data.frame(left = xleft, right = xright)
  > f1 <- mledist(d, "norm")
  > 
  > dbis <- data.frame(left = c(-1.8, -0.6, -0.1, 0.07, 0.14, 1, 1.2), 
  +                    right = c(-1.8, -0.6, -0.1, 0.07, 0.14, 1, NA))
  > f2 <- mledist(dbis, "norm", weights = c(rep(1,6),3))
  Warning message:
  In mledist(dbis, "norm", weights = c(rep(1, 6), 3)) :
    weights are not taken into account in the default initial values
  > 
  > # f1 and f2 must give quite the same results (only starting values differ)
  > f1$estimate
       mean        sd 
  0.5447786 1.3422448 
  > f2$estimate
       mean        sd 
  0.5448978 1.3423076 
  > 
  > # (4) test the definition of fix.arg/start.arg as functions
  > #     if defined as functions, start.arg and fix.arg must be
  > #     functions of pseudo data (output pseudo of cens2pseudo())
  > 
  > mledist(d, "norm", start = function(x) list(mean = 0, sd = 1))$estimate
       mean        sd 
  0.5448718 1.3420260 
  > mledist(d, "norm", start = function(x) list(mean = mean(x), sd = 1))$estimate
       mean        sd 
  0.5447851 1.3422520 
  > mledist(d, "norm", fix.arg = function(x) list(mean = 0))$estimate
        sd 
  1.399337 
  > mledist(d, "norm", fix.arg = function(x) list(mean = 0.544))$estimate
        sd 
  1.341995 
  > mledist(d, "norm", fix.arg = function(x) list(mean = mean(x)))$estimate
        sd 
  1.327269 
  > 
  > 
  > proc.time()
     user  system elapsed 
     1.51    0.15    1.65 
Running the tests in 'tests/t-msedist.R' failed.
Complete output:
  > library(fitdistrplus)
  Loading required package: MASS
  Loading required package: survival
  > nboot <- 1000
  > nboot <- 10
  > nsample <- 10
  > 
  > set.seed(123)
  > 
  > #manual implementation
  > geomeanspacing <- function(pdist, obs, ...)
  + {
  +   sx <- c(-Inf, sort(obs), Inf)
  +   n <- length(sx)
  +   Di <- pdist(sx[-1], ...) - pdist(sx[-n], ...)
  +   mean(log(Di))
  + }
  > 
  > #--------------------------------------------------------
  > # exponential sample
  > 
  > x1 <- rexp(nsample)
  > lam <- seq(0.1, 2, length=51)
  > Sn <- sapply(lam, function(x) geomeanspacing(pexp, obs=x1, rate=x))
  > 
  > Dn <- function(theta)
  +   -geomeanspacing(pexp, obs=x1, rate=theta[1])
  > #theoretical optimum
  > Dn(1/mean(x1))
  [1] 3.090208
  > 
  > #check curve behavior
  > par(mar=c(4,4,2,1))
  > plot(lam, Sn, type="l", xlab="theta", main="Average spacing logarithm")
  > abline(v=1, col="green")
  > abline(v=msedist(x1, "exp")$estimate, lty=2, col="blue")
  > legend("bottomright", lty=1:2, col=c("green", "blue"), leg=c("theoretical value", "fitted value"))
  > 
  > msedist(x1, "exp", control=list(trace=0, REPORT=1))
  $estimate
      rate 
  1.336924 
  
  $convergence
  [1] 0
  
  $value
  [1] 3.078204
  
  $hessian
            rate
  rate 0.4930124
  
  $optim.function
  [1] "optim"
  
  $optim.method
  [1] "BFGS"
  
  $fix.arg
  NULL
  
  $fix.arg.fun
  NULL
  
  $weights
  NULL
  
  $counts
  function gradient 
        11        6 
  
  $optim.message
  NULL
  
  $loglik
  [1] -5.610213
  
  $phidiv
  [1] "KL"
  
  $power.phidiv
  NULL
  
  > 
  > mse_exp <- fitdist(x1, "exp", method="mse")
  > plot(mse_exp)
  > summary(mse_exp)
  Parameters : 
       estimate
  rate 1.336924
  Loglikelihood:  -5.610213   AIC:  13.22043   BIC:  13.52301 
  > gofstat(mse_exp)
  Goodness-of-fit statistics
                               1-mse-exp
  Kolmogorov-Smirnov statistic 0.2549881
  Cramer-von Mises statistic   0.1544480
  Anderson-Darling statistic   1.0571407
  
  Goodness-of-fit criteria
                                 1-mse-exp
  Akaike's Information Criterion  13.22043
  Bayesian Information Criterion  13.52301
  > 
  > mse_exp_boot <- bootdist(mse_exp, niter = nboot)
  > plot(mse_exp_boot)
  > abline(v=1, col="green")
  > abline(v=msedist(x1, "exp")$estimate, lty=2, col="blue")
  > legend("bottomright", lty=1:2, col=c("green", "blue"), leg=c("theoretical value", "fitted value"))
  > 
  > 
  > # library(BMT)
  > # x <- rBMT(1e3, 1/4, 3/4)
  > # 
  > # BMTfit.mpse(x)
  > # fitdist(x, "BMT", method="mse", start=list(p3=1/2, p4=1/2, p1=-1/2, p2=1/2), lower=c(0, 0, -Inf, 0), 
  > #         upper=c(1,1,0, Inf))
  > # pBMT(x, p3=1/2, p4=1/2, p1=-1/2, p2=1/2)
  > 
  > #--------------------------------------------------------
  > # lognormal sample
  > 
  > 
  > x1 <- rlnorm(nsample, 0, 1)
  > mu <- seq(-1, 1, length=51)
  > Sn <- sapply(mu, 
  +              function(x) geomeanspacing(plnorm, obs=x1, mean=x, sd=1))
  > 
  > 
  > Dn <- function(theta)
  +   -geomeanspacing(plnorm, obs=x1, mean=theta[1], sd=theta[2])
  > 
  > plot(mu, Sn, type="l")
  > abline(v=0)
  > 
  > 
  > optim(c(2,2), Dn)
  $par
  [1] -0.5939607  0.6723916
  
  $value
  [1] 2.614255
  
  $counts
  function gradient 
        81       NA 
  
  $convergence
  [1] 0
  
  $message
  NULL
  
  > msedist(x1, "lnorm", control=list(trace=0, REPORT=1))
  $estimate
     meanlog      sdlog 
  -0.5939281  0.6723368 
  
  $convergence
  [1] 0
  
  $value
  [1] 2.614255
  
  $hessian
             meanlog      sdlog
  meanlog 2.12465975 0.01259491
  sdlog   0.01259491 3.14939812
  
  $optim.function
  [1] "optim"
  
  $optim.method
  [1] "Nelder-Mead"
  
  $fix.arg
  NULL
  
  $fix.arg.fun
  NULL
  
  $weights
  NULL
  
  $counts
  function gradient 
        47       NA 
  
  $optim.message
  NULL
  
  $loglik
  [1] -2.380166
  
  $phidiv
  [1] "KL"
  
  $power.phidiv
  NULL
  
  > 
  > 
  > mse_lnorm <- fitdist(x1, "lnorm", method="mse")
  > mle_lnorm <- fitdist(x1, "lnorm", method="mle")
  > plot(mse_lnorm)
  > summary(mse_lnorm)
  Parameters : 
            estimate
  meanlog -0.5939281
  sdlog    0.6723368
  Loglikelihood:  -2.380166   AIC:  8.760331   BIC:  9.365502 
  > cdfcomp(list(mse_lnorm, mle_lnorm))
  > gofstat(list(mse_lnorm, mle_lnorm))
  Goodness-of-fit statistics
                               1-mse-lnorm 2-mle-lnorm
  Kolmogorov-Smirnov statistic  0.16002509  0.15515574
  Cramer-von Mises statistic    0.03898949  0.02565167
  Anderson-Darling statistic    0.26173571  0.18525761
  
  Goodness-of-fit criteria
                                 1-mse-lnorm 2-mle-lnorm
  Akaike's Information Criterion    8.760331    7.735518
  Bayesian Information Criterion    9.365502    8.340688
  > mse_lnorm_boot <- bootdist(mse_lnorm, niter = nboot)
  > par(mar=c(4,4,2,1))
  > plot(mse_lnorm_boot, enhance = TRUE, trueval=c(0,1))
  > 
  > #--------------------------------------------------------
  > # Pareto sample
  > 
  > library(actuar)
  
  Attaching package: 'actuar'
  
  The following objects are masked from 'package:stats':
  
      sd, var
  
  The following object is masked from 'package:grDevices':
  
      cm
  
  > 
  > x1 <- rburr(nsample, 2,2,2)
  > 
  > Dn <- function(theta)
  +   -geomeanspacing(pburr, obs=x1, shape1=theta[1], shape2=theta[2], rate=theta[3])
  > Dn(c(1,1,10))
  [1] 3.408594
  > 
  > optim(c(1,1,10), Dn)
  $par
  [1] 0.609447 3.343173 3.811287
  
  $value
  [1] 2.738615
  
  $counts
  function gradient 
       126       NA 
  
  $convergence
  [1] 0
  
  $message
  NULL
  
  Warning messages:
  1: In pdist(sx[-1], ...) : NaNs produced
  2: In pdist(sx[-n], ...) : NaNs produced
  3: In pdist(sx[-1], ...) : NaNs produced
  4: In pdist(sx[-n], ...) : NaNs produced
  > 
  > msedist(x1, "burr", start=list(shape1=1, shape2=1, rate=10), control=list(trace=0, REPORT=1))
  $estimate
    shape1   shape2     rate 
  0.609447 3.343173 3.811287 
  
  $convergence
  [1] 0
  
  $value
  [1] 2.738615
  
  $hessian
            shape1     shape2       rate
  shape1 2.3885660 0.29619164 0.53774428
  shape2 0.2961916 0.09251466 0.02910551
  rate   0.5377443 0.02910551 0.19023261
  
  $optim.function
  [1] "optim"
  
  $optim.method
  [1] "Nelder-Mead"
  
  $fix.arg
  NULL
  
  $fix.arg.fun
  NULL
  
  $weights
  NULL
  
  $counts
  function gradient 
       126       NA 
  
  $optim.message
  NULL
  
  $loglik
  [1] 3.023078
  
  $phidiv
  [1] "KL"
  
  $power.phidiv
  NULL
  
  > 
  > 
  > mse_burr <- fitdist(x1, "burr", method="mse", start=list(shape1=1, shape2=1, rate=10))
  > mle_burr <- fitdist(x1, "burr", method="mle", start=list(shape1=1, shape2=1, rate=10))
  > plot(mse_burr)
  > summary(mse_burr)
  Parameters : 
         estimate
  shape1 0.609447
  shape2 3.343173
  rate   3.811287
  Loglikelihood:  3.023078   AIC:  -0.046155   BIC:  0.8616003 
  > cdfcomp(list(mse_burr, mle_burr))
  > gofstat(list(mse_burr, mle_burr))
  Goodness-of-fit statistics
                               1-mse-burr 2-mle-burr
  Kolmogorov-Smirnov statistic 0.20751258 0.15833478
  Cramer-von Mises statistic   0.07113128 0.03458676
  Anderson-Darling statistic   0.40398658 0.28801421
  
  Goodness-of-fit criteria
                                 1-mse-burr  2-mle-burr
  Akaike's Information Criterion -0.0461550 -0.99734543
  Bayesian Information Criterion  0.8616003 -0.08959015
  > mse_burr_boot <- bootdist(mse_burr, niter = pmin(nboot,100))
  > plot(mse_burr_boot, enhance = TRUE, trueval=c(2,2,2))
  > 
  > 
  > 
  > #--------------------------------------------------------
  > # Poisson sample
  > 
  > x1 <- rpois(nsample, 15)
  > 
  > geomeanSpacingUnique <- function(pdist, obs, ...)
  + {
  +   sx <- c(-Inf, unique(sort(obs)), Inf)
  +   n <- length(sx)
  +   Di <- pdist(sx[-1], ...) - pdist(sx[-n], ...)
  +   mean(log(Di))
  + }
  > 
  > geomeanSpacingWeight <- function(pdist, obs, weights, ...)
  + {
  +   sx <- c(-Inf, unique(sort(obs)), Inf)
  +   weights <- c(1, weights)
  +   n <- length(sx)
  +   Di <- pdist(sx[-1], ...) - pdist(sx[-n], ...)
  +   mean(weights*log(Di))
  + }
  > 
  > DnUnique <- function(theta)
  +   -geomeanSpacingUnique(ppois, obs=x1, lambda=theta[1])
  > 
  > DnWeight <- function(theta, weights)
  +   -geomeanSpacingWeight(ppois, obs=x1, lambda=theta[1], weights=weights)
  > 
  > 
  > optimize(DnWeight, c(1, 30), weights=as.numeric(table(x1)))
  $minimum
  [1] 15.86236
  
  $objective
  [1] 2.796532
  
  > optimize(DnUnique, c(1, 30))
  $minimum
  [1] 15.49055
  
  $objective
  [1] 2.581978
  
  > optimize(Dn, c(1, 30)) #does not converge
  $minimum
  [1] 29.99993
  
  $objective
  [1] NA
  
  There were 27 warnings (use warnings() to see them)
  > 
  > mle_pois1 <- fitdist(x1, "pois", method="mle")
  > #no weight
  > mse_pois1 <- fitdist(x1, "pois", method="mse")
  > #with weight
  > mse_pois2 <- fitdist(unique(sort(x1)), "pois", method="mse", weights=as.numeric(table(x1)))
  Warning message:
  In msedist(data, distname, start = arg_startfix$start.arg, fix.arg = arg_startfix$fix.arg,  :
    weights are not taken into account in the default initial values
  > plot(mse_pois1)
  > plot(mse_pois2)
  > summary(mse_pois1)
  Parameters : 
         estimate
  lambda 15.49046
  Loglikelihood:  -29.76579   AIC:  61.53157   BIC:  61.83416 
  > gofstat(mse_pois1)
  Chi-squared statistic:  6.798192 
  Degree of freedom of the Chi-squared distribution:  3 
  Chi-squared p-value:  0.07861594 
     the p-value may be wrong with some theoretical counts < 5  
  Chi-squared table:
        obscounts theocounts
  <= 9  2.0000000  0.5544272
  <= 13 2.0000000  2.6248528
  <= 17 2.0000000  3.8811167
  <= 18 2.0000000  0.7716062
  > 18  2.0000000  2.1679971
  
  Goodness-of-fit criteria
                                 1-mse-pois
  Akaike's Information Criterion   61.53157
  Bayesian Information Criterion   61.83416
  > gofstat(mse_pois2)
  Chi-squared statistic:  6.650898 
  Degree of freedom of the Chi-squared distribution:  2 
  Chi-squared p-value:  0.03595636 
     the p-value may be wrong with some theoretical counts < 5  
  Chi-squared table:
        obscounts theocounts
  <= 9  2.0000000  0.4168917
  <= 13 2.0000000  2.1558228
  <= 17 2.0000000  3.4765332
  > 17  3.0000000  2.9507523
  
  Goodness-of-fit criteria
                                 1-mse-pois
  Akaike's Information Criterion   56.74383
  Bayesian Information Criterion   56.94105
  > 
  > par(mfrow=c(1,1))
  > cdfcomp(list(mle_pois1, mse_pois1), addlegend = FALSE, fitlty = 1)
  > curve(ppois(x, lambda=mse_pois2$estimate), type="s", col="blue", add=TRUE)
  > legend("bottomright", lty=1, col=c("red", "green", "blue"), leg=c("MLE", "MSE no weight", "MSE with weight"))
  > 
  > 
  > 
  > #--------------------------------------------------------
  > # real dataset
  > # library(CASdatasets)
  > # data("ushustormloss")
  > # x <- ushustormloss$Normalized.CL05
  > # 
  > # plot(Normalized.CL05 ~ Year, data=ushustormloss, type="h", main="Normalized Hurricane Damages in United States")
  > # 
  > # mse_burr <- fitdist(x, "burr", method="mse", start=list(shape1=1, shape2=1, rate=10), lower=0)
  > # mle_burr0 <- fitdist(x, "burr", method="mle", start=list(shape1=1, shape2=1, rate=10), lower=0)
  > # 
  > # cbind(MSE=coef(mse_burr), MLE=coef(mle_burr0))
  > # 
  > # 
  > # setwd("~/Desktop/")
  > # par(mar=c(4,4,2,1))
  > # pdf("Ushustorm-cdfcomp.pdf", 6, 6)
  > # cdfcomp(list(mse_burr, mle_burr0), xlogscale = TRUE, do.points = FALSE)
  > # dev.off()
  > # pdf("Ushustorm-qqcomp.pdf", 6, 6)
  > # qqcomp(list(mse_burr, mle_burr0), xlogscale=TRUE, ylogscale=TRUE)
  > # dev.off()
  > # pdf("Ushustorm-ppcomp.pdf", 6, 6)
  > # ppcomp(list(mse_burr, mle_burr0))
  > # dev.off()
  > # 
  > # gofstat(list(mse_burr, mle_burr0))
  > # 
  > # mse_iburr <- fitdist(x, "invburr", method="mse", start=list(shape1=1, shape2=1, rate=10), lower=0)
  > # mle_iburr0 <- fitdist(x, "invburr", method="mle", start=list(shape1=1, shape2=1, rate=10), lower=0)
  > # 
  > # gofstat(list(mse_iburr, mle_iburr0))
  > # cdfcomp(list(mse_iburr, mle_iburr0))
  > 
  > 
  > proc.time()
     user  system elapsed 
     2.42    0.25    2.64 
Running the tests in 'tests/t-parallel.R' failed.
Complete output:
  > visualize <- FALSE # TRUE for manual tests with visualization of results
  > 
  > if(visualize)
  + {
  +   
  +   library(parallel)
  +   
  +   
  +   #fonction basique evaluant la moyenne empirique d'un echantillon gaussien
  +   f <- function(i) mean(rnorm(1e4, mean=i))
  +   
  +   #exemple simple
  +   glist <- 1:20 	
  +   cl <- parallel::makeCluster(2)
  +   system.time(
  +     res <- parallel::parLapply(cl, glist, f)
  +   )
  +   parallel::stopCluster(cl)
  +   
  +   
  +   #exemple en faisant varier le nombre de coeurs et le nombre de simulations
  +   nbsimu <- 10^(1:2)
  +   cores <- 1:4
  +   cores <- 1:getOption("cl.cores", 2)
  +   
  +   partime <- matrix(NA, length(nbsimu), length(cores)+1)
  +   colnames(partime) <- c("R", paste("core",cores))
  +   rownames(partime) <- paste("n", nbsimu, sep="=")
  +   
  +   partime[, 1] <- sapply(1:length(nbsimu), function(i) 
  +     system.time(lapply(1:nbsimu[i], f))[3])
  +   for(j in 1:length(cores))
  +   {
  +     print(cores[j])
  +     cl <- parallel::makeCluster(cores[j])
  +     partime[, j+1] <- sapply(1:length(nbsimu), function(i) 
  +       system.time(parallel::parLapply(cl, 1:nbsimu[i], f))[3])
  +     parallel::stopCluster(cl)
  +   }		
  +   
  +   
  +   partime
  +   
  + }
  > 
  > proc.time()
     user  system elapsed 
     0.17    0.06    0.20 
Running the tests in 'tests/t-plotdistcens.R' failed.
Complete output:
  > library(fitdistrplus)
  Loading required package: MASS
  Loading required package: survival
  > 
  > # (1) Plot of an empirical censored distribution (censored data) as a CDF
  > # using the default Turnbull method
  > #
  > data(smokedfish)
  > plotdistcens(smokedfish)
  > d1 <- as.data.frame(log10(smokedfish))
  > plotdistcens(d1)
  > 
  > 
  > # (2) Add the CDF of a normal distribution and QQ and PP plots
  > #
  > plotdistcens(smokedfish,"lnorm", para=list(meanlog=-3.6,sdlog=3.5))
  > plotdistcens(d1,"norm", para=list(mean=-1.6,sd=1.5))
  > 
  > # (3) Various plots of the same empirical distribution 
  > #
  > # default Wang plot
  > plotdistcens(d1, NPMLE = TRUE, NPMLE.method = "Wang")
  > plotdistcens(d1, NPMLE = TRUE, NPMLE.method = "Wang", lwd = 3, main = "Wang ECDF plot")
  > # Turnbull plot
  > plotdistcens(d1, NPMLE = TRUE, NPMLE.method = "Turnbull", col = "red", 
  +              main = "Turnbull ECDF plot")
  Warning message:
  In plotdistcens(d1, NPMLE = TRUE, NPMLE.method = "Turnbull", col = "red",  :
    Turnbull is now a deprecated option for NPMLE.method. You should use Turnbull.middlepoints
              of Turnbull.intervals. It was here fixed as Turnbull.middlepoints, equivalent to former Turnbull.
  > plotdistcens(d1,Turnbull = TRUE) # deprecated way to do it
  Warning message:
  In plotdistcens(d1, Turnbull = TRUE) :
    When Turnbull.confint is TRUE NPMLE.method is forced to Turnbull.middlepoints.
  > # Turnbull plot with confidence intervals
  > plotdistcens(d1,NPMLE = TRUE, NPMLE.method = "Turnbull", Turnbull.confint = TRUE)
  Warning message:
  In plotdistcens(d1, NPMLE = TRUE, NPMLE.method = "Turnbull", Turnbull.confint = TRUE) :
    Turnbull is now a deprecated option for NPMLE.method. You should use Turnbull.middlepoints
              of Turnbull.intervals. It was here fixed as Turnbull.middlepoints, equivalent to former Turnbull.
  > plotdistcens(d1,Turnbull = TRUE,Turnbull.confint = TRUE) # deprecated way to do it
  Warning messages:
  1: In plotdistcens(d1, Turnbull = TRUE, Turnbull.confint = TRUE) :
    When Turnbull.confint is TRUE NPMLE.method is forced to Turnbull.middlepoints.
  2: In plot.window(...) : "Turnbull" is not a graphical parameter
  3: In plot.xy(xy, type, ...) : "Turnbull" is not a graphical parameter
  4: In axis(side = side, at = at, labels = labels, ...) :
    "Turnbull" is not a graphical parameter
  5: In axis(side = side, at = at, labels = labels, ...) :
    "Turnbull" is not a graphical parameter
  6: In box(...) : "Turnbull" is not a graphical parameter
  7: In title(...) : "Turnbull" is not a graphical parameter
  > # with intervals and points
  > plotdistcens(d1,NPMLE = FALSE)
  > plotdistcens(d1,NPMLE = FALSE, col = "red", lwd = 2)
  > plotdistcens(d1,rightNA=3, NPMLE = FALSE)
  > plotdistcens(d1,rightNA=3, Turnbull = FALSE) # deprecated way to do it
  > # with intervals and points
  > # defining a minimum value for left censored values
  > plotdistcens(d1,leftNA=-3, NPMLE = FALSE)
  > 
  > # (4) Goodness-of-fit plots for the same dataset after logarithmic transformation
  > #   with a lognormal distribution, successively using the three proposed methods
  > #
  > d3 <- smokedfish
  > plotdistcens(d3,"lnorm",para=list(meanlog=-3.6,sdlog=3.5), main = "Wang plot")
  > plotdistcens(d3,"lnorm",para=list(meanlog=-3.6,sdlog=3.5), 
  +              NPMLE.method = "Turnbull", main = "Turnbull plot")
  Warning messages:
  1: In plotdistcens(d3, "lnorm", para = list(meanlog = -3.6, sdlog = 3.5),  :
    Turnbull is now a deprecated option for NPMLE.method. You should use Turnbull.middlepoints
              of Turnbull.intervals. It was here fixed as Turnbull.middlepoints, equivalent to former Turnbull.
  2: In plotdistcens(d3, "lnorm", para = list(meanlog = -3.6, sdlog = 3.5),  :
    Q-Q plot and P-P plot are available only  
              with the arguments NPMLE.method at Wang (default value) or Turnbull.intervals.
  > plotdistcens(d3,"lnorm",para=list(meanlog=-3.6,sdlog=3.5),
  +              NPMLE = FALSE, leftNA=0, main = "Plot of ordered intervals")
  Warning message:
  In plotdistcens(d3, "lnorm", para = list(meanlog = -3.6, sdlog = 3.5),  :
    When NPMLE is FALSE the nonparametric maximum likelihood estimation 
              of the cumulative distribution function is not computed.
              Q-Q plot and P-P plot are available only with the arguments NPMLE.method at Wang 
              (default value) or Turnbull.intervals.
  > 
  > # Test with the salinity data set
  > #
  > data(salinity)
  > log10LC50 <-log10(salinity)
  > plotdistcens(log10LC50)
  > plotdistcens(log10LC50, NPMLE.method = "Turnbull")
  Warning message:
  In plotdistcens(log10LC50, NPMLE.method = "Turnbull") :
    Turnbull is now a deprecated option for NPMLE.method. You should use Turnbull.middlepoints
              of Turnbull.intervals. It was here fixed as Turnbull.middlepoints, equivalent to former Turnbull.
  > plotdistcens(log10LC50, NPMLE = FALSE)
  > fn <- fitdistcens(log10LC50,"norm")
  > fl <- fitdistcens(log10LC50,"logis")
  > plot(fn)
  > plot(fl)
  > 
  > proc.time()
     user  system elapsed 
     1.89    0.17    2.04 
Running the tests in 'tests/t-prefit.R' failed.
Complete output:
  > require(fitdistrplus)
  Loading required package: fitdistrplus
  Loading required package: MASS
  Loading required package: survival
  > nsample <- 10
  > 
  > # (1) gamma
  > x <- rgamma(nsample, 5/2, 7/2)
  > 
  > prefit(x, "gamma", "mle", list(shape=3, scale=3), lower=-Inf, upper=Inf, silent=TRUE, control=list(trace=1, REPORT=1))
  initial  value 6.560080 
  iter   2 value 4.236045
  iter   3 value 3.779057
  iter   4 value 3.375129
  iter   5 value 3.344779
  iter   6 value 3.341388
  iter   7 value 3.341369
  iter   8 value 3.341361
  iter   9 value 3.341353
  iter  10 value 3.341353
  iter  10 value 3.341353
  iter  10 value 3.341353
  final  value 3.341353 
  converged
  $shape
  [1] 2.422232
  
  $scale
  [1] 3.954371
  
  > prefit(x, "gamma", "mle", list(shape=1, scale=1), lower=-Inf, upper=Inf, silent=TRUE)
  $shape
  [1] 2.422231
  
  $scale
  [1] 3.954373
  
  > 
  > prefit(x, "gamma", "mle", list(shape=3), fix.arg=list(scale=7/2), lower=-Inf, upper=Inf, silent=TRUE)
  $shape
  [1] 0.5000425
  
  > 
  > prefit(x, "gamma", "qme", list(shape=1, scale=1), probs=1:2/3, lower=-Inf, upper=Inf, silent=TRUE)
  $shape
  [1] 10.70754
  
  $scale
  [1] 19.73535
  
  > 
  > prefit(x, "gamma", "mge", list(shape=1, scale=1), gof="CvM", lower=-Inf, upper=Inf, silent=TRUE)
  $shape
  [1] 2.881339
  
  $scale
  [1] 4.976772
  
  > prefit(x, "gamma", "mge", list(shape=1, scale=1), gof="AD", lower=-Inf, upper=Inf, silent=TRUE)
  $shape
  [1] 2.145594
  
  $scale
  [1] 3.463989
  
  > 
  > # (2) geometric
  > x <- rgeom(nsample, 1/7)
  > prefit(x, "geom", "mle", list(prob=1/2), lower=-Inf, upper=Inf, silent=TRUE)
  $prob
  [1] 0.1587296
  
  > tbx <- table(x)
  > prefit(as.numeric(names(tbx)), "geom", "mle", list(prob=1/2), lower=-Inf, upper=Inf, silent=TRUE, weights=tbx)
  $prob
  [1] 0.1587296
  
  > prefit(x, "geom", "qme", list(prob=1/2), probs=1/2, lower=-Inf, upper=Inf)
  $prob
  [1] 0.3100255
  
  > 
  > 
  > # (3) Pareto
  > require(actuar)
  Loading required package: actuar
  
  Attaching package: 'actuar'
  
  The following objects are masked from 'package:stats':
  
      sd, var
  
  The following object is masked from 'package:grDevices':
  
      cm
  
  > x  <-  rpareto(nsample, 6, 2)
  > 
  > prefit(x, "pareto", "mme", list(shape=10, scale=10), order=1:2, memp=function(x, order) mean(x^order), lower=-Inf, upper=Inf)
  $shape
  [1] 32.52952
  
  $scale
  [1] 12.40677
  
  > 
  > 
  > proc.time()
     user  system elapsed 
     1.46    0.26    1.70 
Running the tests in 'tests/t-qme-discrete.R' failed.
Complete output:
  > library(fitdistrplus)
  Loading required package: MASS
  Loading required package: survival
  > set.seed(1234)
  > nsample <- 10
  > 
  > #Poisson
  > 
  > x <- rpois(nsample, lambda=7.5)
  > L2 <- function(lam)
  +   (qpois(1/2, lambda = lam) - median(x))^2
  > curve(L2(x), 5, 9, xlab=expression(lambda), ylab=expression(L2(lambda)), main="squared differences", n=201)
  > 
  > fitdist(x, "pois", method="qme", probs=1/2, start=list(lambda=2), control=list(trace=1, REPORT=1))
  initial  value 36.000000 
  iter   1 value 36.000000
  final  value 36.000000 
  converged
  Fitting of the distribution ' pois ' by matching quantiles 
  Parameters:
         estimate
  lambda        2
  > fitdist(x, "pois", method="qme", probs=1/2, start=list(lambda=6.8), control=list(trace=1, REPORT=1))
  initial  value 1.000000 
  iter   1 value 1.000000
  final  value 1.000000 
  converged
  Fitting of the distribution ' pois ' by matching quantiles 
  Parameters:
         estimate
  lambda      6.8
  > fitdist(x, "pois", method="qme", probs=1/2, start=list(lambda=15), control=list(trace=1, REPORT=1))
  initial  value 49.000000 
  iter   1 value 49.000000
  final  value 49.000000 
  converged
  Fitting of the distribution ' pois ' by matching quantiles 
  Parameters:
         estimate
  lambda       15
  > 
  > fitdist(x, "pois", method="qme", optim.method="SANN", probs=1/2, start=list(lambda=2), control=list(trace=1, REPORT=100))
  sann objective function values
  initial       value 36.000000
  iter     1000 value 0.000000
  iter     2000 value 0.000000
  iter     3000 value 0.000000
  iter     4000 value 0.000000
  iter     5000 value 0.000000
  iter     6000 value 0.000000
  iter     7000 value 0.000000
  iter     8000 value 0.000000
  iter     9000 value 0.000000
  iter     9999 value 0.000000
  final         value 0.000000
  sann stopped after 9999 iterations
  Fitting of the distribution ' pois ' by matching quantiles 
  Parameters:
         estimate
  lambda 8.425269
  > fitdist(x, "pois", method="qme", optim.method="SANN", probs=1/2, start=list(lambda=17), control=list(trace=1, REPORT=100))
  sann objective function values
  initial       value 81.000000
  iter     1000 value 0.000000
  iter     2000 value 0.000000
  iter     3000 value 0.000000
  iter     4000 value 0.000000
  iter     5000 value 0.000000
  iter     6000 value 0.000000
  iter     7000 value 0.000000
  iter     8000 value 0.000000
  iter     9000 value 0.000000
  iter     9999 value 0.000000
  final         value 0.000000
  sann stopped after 9999 iterations
  Fitting of the distribution ' pois ' by matching quantiles 
  Parameters:
         estimate
  lambda 8.183321
  > 
  > 
  > 
  > #Geometric
  > 
  > x <- rgeom(nsample, 1/3)
  > L2 <- function(p)
  +   (qgeom(1/2, p) - median(x))^2
  > curve(L2(x), 0.10, 0.95, xlab=expression(p), ylab=expression(L2(p)), main="squared differences", n=301)
  > 
  > fitdist(x, "geom", method="qme", probs=1/2, start=list(prob=1/2), control=list(trace=1, REPORT=1))
  initial  value 16.000000 
  iter   2 value 4.000000
  iter   2 value 4.000000
  iter   2 value 4.000000
  final  value 4.000000 
  converged
  Fitting of the distribution ' geom ' by matching quantiles 
  Parameters:
       estimate
  prob    0.276
  > fitdist(x, "geom", method="qme", probs=1/2, start=list(prob=1/20), control=list(trace=1, REPORT=1))
  initial  value 81.000000 
  iter   1 value 81.000000
  final  value 81.000000 
  converged
  Fitting of the distribution ' geom ' by matching quantiles 
  Parameters:
       estimate
  prob     0.05
  > 
  > proc.time()
     user  system elapsed 
     6.59    0.21    6.81 
Running the tests in 'tests/t-qmedist.R' failed.
Complete output:
  > library(fitdistrplus)
  Loading required package: MASS
  Loading required package: survival
  > nsample <- 10
  > 
  > # (1) basic fit of a normal distribution 
  > #
  > 
  > set.seed(1234)
  > x1 <- rnorm(n=nsample)
  > qmedist(x1, "norm", probs=c(1/3, 2/3))
  $estimate
        mean         sd 
  -0.1486435  0.9892013 
  
  $convergence
  [1] 0
  
  $value
  [1] 4.312801e-10
  
  $hessian
               mean           sd
  mean 4.000000e+00 5.548066e-14
  sd   5.548066e-14 7.421040e-01
  
  $optim.function
  [1] "optim"
  
  $optim.method
  [1] "Nelder-Mead"
  
  $fix.arg
  NULL
  
  $fix.arg.fun
  NULL
  
  $weights
  NULL
  
  $counts
  function gradient 
        55       NA 
  
  $optim.message
  NULL
  
  $loglik
  [1] -13.92195
  
  $probs
  [1] 0.3333333 0.6666667
  
  > 
  > 
  > # (2) defining your own distribution functions, here for the Gumbel 
  > # distribution for other distributions, see the CRAN task view dedicated 
  > # to probability distributions
  > 
  > dgumbel <- function(x, a, b) 1/b*exp((a-x)/b)*exp(-exp((a-x)/b))
  > qgumbel <- function(p, a, b) a - b*log(-log(p))
  > qmedist(x1, "gumbel", probs=c(1/3, 2/3), start=list(a=10,b=5))
  $estimate
           a          b 
  -0.4935571  0.8557281 
  
  $convergence
  [1] 0
  
  $value
  [1] 2.735865e-06
  
  $hessian
           a        b
  a 4.000000 1.617345
  b 1.617345 1.647498
  
  $optim.function
  [1] "optim"
  
  $optim.method
  [1] "Nelder-Mead"
  
  $fix.arg
  NULL
  
  $fix.arg.fun
  NULL
  
  $weights
  NULL
  
  $counts
  function gradient 
        61       NA 
  
  $optim.message
  NULL
  
  $loglik
  [1] -16.79828
  
  $probs
  [1] 0.3333333 0.6666667
  
  > 
  > # (3) fit a discrete distribution (Poisson)
  > #
  > 
  > set.seed(1234)
  > x2 <- rpois(n=nsample,lambda = 2)
  > qmedist(x2, "pois", probs=1/2)
  $estimate
  lambda 
     1.7 
  
  $convergence
  [1] 0
  
  $value
  [1] 0
  
  $hessian
         lambda
  lambda      0
  
  $optim.function
  [1] "optim"
  
  $optim.method
  [1] "BFGS"
  
  $fix.arg
  NULL
  
  $fix.arg.fun
  NULL
  
  $weights
  NULL
  
  $counts
  function gradient 
         1        1 
  
  $optim.message
  NULL
  
  $loglik
  [1] -15.31626
  
  $probs
  [1] 0.5
  
  > 
  > # (4) fit a finite-support distribution (beta)
  > #
  > 
  > set.seed(1234)
  > x3 <- rbeta(n=nsample, shape1=5, shape2=10)
  > qmedist(x3, "beta", probs=c(1/3, 2/3))
  $estimate
    shape1   shape2 
  4.010999 8.397853 
  
  $convergence
  [1] 0
  
  $value
  [1] 4.99688e-11
  
  $hessian
              shape1       shape2
  shape1  0.01325492 -0.006123230
  shape2 -0.00612323  0.002877022
  
  $optim.function
  [1] "optim"
  
  $optim.method
  [1] "Nelder-Mead"
  
  $fix.arg
  NULL
  
  $fix.arg.fun
  NULL
  
  $weights
  NULL
  
  $counts
  function gradient 
        61       NA 
  
  $optim.message
  NULL
  
  $loglik
  [1] 6.642124
  
  $probs
  [1] 0.3333333 0.6666667
  
  > 
  > 
  > # (5) fit frequency distributions on USArrests dataset.
  > #
  > 
  > x4 <- USArrests$Assault
  > qmedist(x4, "pois", probs=1/2)
  $estimate
  lambda 
  170.76 
  
  $convergence
  [1] 0
  
  $value
  [1] 144
  
  $hessian
         lambda
  lambda      0
  
  $optim.function
  [1] "optim"
  
  $optim.method
  [1] "BFGS"
  
  $fix.arg
  NULL
  
  $fix.arg.fun
  NULL
  
  $weights
  NULL
  
  $counts
  function gradient 
         1        1 
  
  $optim.message
  NULL
  
  $loglik
  [1] -1211.705
  
  $probs
  [1] 0.5
  
  > qmedist(x4, "nbinom", probs=c(1/3, 2/3))
  $estimate
        size         mu 
    2.518966 182.313344 
  
  $convergence
  [1] 0
  
  $value
  [1] 0.2222222
  
  $hessian
       size mu
  size    0  0
  mu      0  0
  
  $optim.function
  [1] "optim"
  
  $optim.method
  [1] "Nelder-Mead"
  
  $fix.arg
  NULL
  
  $fix.arg.fun
  NULL
  
  $weights
  NULL
  
  $counts
  function gradient 
        37       NA 
  
  $optim.message
  NULL
  
  $loglik
  [1] -292.5969
  
  $probs
  [1] 0.3333333 0.6666667
  
  > 
  > # (6) normal mixture
  > #
  > 
  > #mixture of two normal distributions
  > #density
  > dnorm2 <- function(x, poid, m1, s1, m2, s2)
  + 	poid*dnorm(x, m1, s1) + (1-poid)*dnorm(x, m2, s2)
  > #numerically approximate quantile function
  > qnorm2 <- function(p, poid, m1, s1, m2, s2)
  + {
  + 	L2 <- function(x, prob)
  + 		(prob - pnorm2(x, poid, m1, s1, m2, s2))^2	
  + 	sapply(p, function(pr) optimize(L2, c(-20, 30), prob=pr)$minimum)
  + }	
  > #distribution function		
  > pnorm2 <- function(q, poid, m1, s1, m2, s2)
  + 	poid*pnorm(q, m1, s1) + (1-poid)*pnorm(q, m2, s2)		
  > 
  > 
  > #basic normal distribution
  > x <- c(rnorm(nsample, 5),  rnorm(nsample, 10))
  > #QME
  > fit2 <- qmedist(x, "norm2", probs=c(1/6, 1/4, 1/3, 1/2, 2/3), 
  + 	start=list(poid=1/3, m1=4, s1=2, m2=8, s2=2), 
  + 	lower=c(0, 0, 0, 0, 0), upper=c(1/2, Inf, Inf, Inf, Inf))
  > 
  > 
  > 
  > # (7) test error messages
  > #
  > 
  > dnorm3 <- qnorm3 <- function(x, a)
  +   "NA"
  > x <- rexp(nsample)
  > 
  > #should get a one-line error 
  > res <- qmedist(x, "norm3", start=list(a=1), probs=1/2)
  <simpleError in qemp - qtheo: non-numeric argument to binary operator>
  > #as in 
  > attr(try(log("a"), silent=TRUE), "condition")
  <simpleError in log("a"): non-numeric argument to mathematical function>
  > 
  > 
  > # (8) weighted QME
  > #
  > n <- 1e6
  > x <- rpois(nsample, 10)
  > xtab <- table(x)
  > xval <- sort(unique(x))
  > f1 <- qmedist(x, "pois", start=list(lambda=mean(x)), lower=0, upper=100, probs=1/2) #, control=list(trace=1, REPORT=1)
  Warning message:
  In qmedist(x, "pois", start = list(lambda = mean(x)), lower = 0,  :
    The BFGS method cannot be used with bounds without provided the gradient. The method is changed to L-BFGS-B.
  > f2 <- qmedist(xval, "pois", weights=xtab, start=list(lambda=mean(x)), probs=1/2)
  Warning message:
  In qmedist(xval, "pois", weights = xtab, start = list(lambda = mean(x)),  :
    weights are not taken into account in the default initial values
  > 
  > f1$estimate
  lambda 
     9.6 
  > f2$estimate #should be identical
  lambda 
     9.6 
  > 
  > x <- rexp(nsample)
  > f3 <- qmedist(x, "exp", probs=1/2)
  > w4 <- c(rep(1, nsample/2), round(sqrt(1:(nsample/2))))
  > f4 <- qmedist(x, "exp", weights=w4, probs=1/2)
  Warning message:
  In qmedist(x, "exp", weights = w4, probs = 1/2) :
    weights are not taken into account in the default initial values
  > f3$estimate
       rate 
  0.4809858 
  > f4$estimate
       rate 
  0.4860816 
  > 
  > f3$loglik
  [1] -14.43196
  > f4$loglik
  [1] -18.91757
  > 
  > median(x)
  [1] 1.441099
  > median(tail(x, 50))
  [1] 1.441099
  > 
  > 
  > #try non integer weights
  > try(qmedist(x, "exp", weights=c(rep(1, n/2), sqrt(1:(n/2))), probs=1/2))
  Error in qmedist(x, "exp", weights = c(rep(1, n/2), sqrt(1:(n/2))), probs = 1/2) : 
    weights should be a vector of (strictly) positive integers
  > 
  > # (9) test the component optim.message
  > x <- rnorm(nsample)
  > #change parameter to obtain unsuccessful convergence
  > qmedist(x, "norm", probs=1:2/3, control=list(maxit=2), start=list(mean=1e5, sd=1), optim.method="L-BFGS-B", lower=0)
  $estimate
       mean        sd 
  0.7855744 1.7499993 
  
  $convergence
  [1] 1
  
  $value
  [1] 0.008861144
  
  $hessian
                mean            sd
  mean  4.000000e+00 -3.903128e-12
  sd   -3.903128e-12  7.421040e-01
  
  $optim.function
  [1] "optim"
  
  $optim.method
  [1] "L-BFGS-B"
  
  $fix.arg
  NULL
  
  $fix.arg.fun
  NULL
  
  $weights
  NULL
  
  $counts
  function gradient 
        11       11 
  
  $optim.message
  [1] "NEW_X"
  
  $loglik
  [1] -17.25942
  
  $probs
  [1] 0.3333333 0.6666667
  
  > 
  > # (10) test bounds
  > x <- rnorm(nsample)
  > qmedist(x, "norm", probs=1:2/3, optim.method="L-BFGS-B", lower=c(-Inf, 0)) #via optim
  $estimate
       mean        sd 
  0.2474441 0.9306178 
  
  $convergence
  [1] 0
  
  $value
  [1] 3.081488e-33
  
  $hessian
                mean            sd
  mean  4.000000e+00 -1.110672e-13
  sd   -1.110672e-13  7.421040e-01
  
  $optim.function
  [1] "optim"
  
  $optim.method
  [1] "L-BFGS-B"
  
  $fix.arg
  NULL
  
  $fix.arg.fun
  NULL
  
  $weights
  NULL
  
  $counts
  function gradient 
         6        6 
  
  $optim.message
  [1] "CONVERGENCE: REL_REDUCTION_OF_F <= FACTR*EPSMCH"
  
  $loglik
  [1] -13.62491
  
  $probs
  [1] 0.3333333 0.6666667
  
  > qmedist(x, "norm", probs=1:2/3, optim.method="Nelder", lower=c(-Inf, 0)) #via constrOptim
  $estimate
       mean        sd 
  0.2474440 0.9306188 
  
  $convergence
  [1] 0
  
  $value
  [1] 4.042279e-13
  
  $hessian
  NULL
  
  $optim.function
  [1] "constrOptim"
  
  $optim.method
  [1] "Nelder-Mead"
  
  $fix.arg
  NULL
  
  $fix.arg.fun
  NULL
  
  $weights
  NULL
  
  $counts
  function gradient 
       124       NA 
  
  $optim.message
  NULL
  
  $loglik
  [1] -13.62491
  
  $probs
  [1] 0.3333333 0.6666667
  
  > 
  > proc.time()
     user  system elapsed 
     4.64    0.32    4.95 
Running the tests in 'tests/t-util-npsurv-mainfunction.R' failed.
Complete output:
  > library(fitdistrplus)
  Loading required package: MASS
  Loading required package: survival
  > vizualise <- FALSE
  > mytrace <- 1
  > mytrace <- 0
  > 
  > #-------------------------------------------------
  > # example with right censoring from package npsurv 
  > ap <- cbind(L=c(1:15, 1:15), 
  +             R=c(1:15, rep(Inf, 15)),
  +             count=c(456, 226, 152, 171, 135, 125,  83,  74,  51,  42,  
  +                     43,  34,  18,   9,   6,  39,  22,  23,  24,  107, 
  +                     133, 102,  68,  64,  45,  53,  33,  27,  23,  30))
  > dim(ap)
  [1] 30  3
  > 
  > resap <- fitdistrplus:::npsurvminimal(ap, w=1, maxit=100, tol=1e-6, verb=mytrace, pkg="stats") 
  > 
  > #------------------------------------------------
  > # example with left censoring from package npsurv 
  > ap <- cbind(L=c(1:15, rep(-Inf, 15)), 
  +             R=c(1:15, 1:15),
  +             count=c(456, 226, 152, 171, 135, 125,  83,  74,  51,  42,  
  +                     43,  34,  18,   9,   6,  39,  22,  23,  24,  107, 
  +                     133, 102,  68,  64,  45,  53,  33,  27,  23,  30))
  > dim(ap)
  [1] 30  3
  > 
  > resap <- fitdistrplus:::npsurvminimal(ap, w=1, maxit=100, tol=1e-6, verb=mytrace, pkg="stats") 
  > 
  > if(vizualise)
  + {
  +   #-------------------------------------------------------------------------------------------------
  +   # example with interval censoring from package npsurv leading to 16 maximal intersection intervals  
  +   ap <- cbind(L=c(0:14,1:15), 
  +               R=c(1:15, rep(Inf, 15)),
  +               count=c(456, 226, 152, 171, 135, 125,  83,  74,  51,  42,  
  +                       43,  34,  18,   9,   6,  39,  22,  23,  24,  107, 
  +                       133, 102,  68,  64,  45,  53,  33,  27,  23,  30))
  +   dim(ap)
  +   
  +   ap.x2 <- fitdistrplus:::icendata(ap, w=1)
  +   ap.x <- rbind(cbind(ap.x2$t, ap.x2$t), ap.x2$o)
  +   ap.Delta <- fitdistrplus:::Deltamatrix(ap.x)
  +   #cbind(ap.Delta$left, unique(ap.x[,"L"]))
  +   
  +   resap <- fitdistrplus:::npsurvminimal(ap, w=1, maxit=100, tol=1e-6, verb=mytrace, pkg="stats") 
  +   
  +   # if(FALSE)
  +   # {
  +   #   require(npsurv)
  +   #   rescheck <- npsurv::npsurv(ap, w=1, maxit=100, tol=1e-6, verb=3) 
  +   #   c(resap$ll, rescheck$ll)
  +   #   
  +   #   cbind(resap$f$left, rescheck$f$left)
  +   #   cbind(resap$f$right, rescheck$f$right)
  +   #   cbind(resap$f$p, rescheck$f$p)
  +   #   sum(abs(resap$f$p- rescheck$f$p))
  +   # }
  +   
  +   
  +   
  +   #----------------------------------------------------------------------------------
  +   # example with interval censoring leading to a single maximal intersection interval
  +   
  +   LR <- matrix(1:100, ncol=2)
  +   cnt <- round(1000*pgeom(1:NROW(LR)-1, prob=1/10, lower=FALSE))
  +   fakedata <- cbind(L=LR[,1], R=LR[,2], count=cnt)
  +   
  +   fakedata.x2 <- fitdistrplus:::icendata(fakedata, w=1)
  +   fakedata.x2.x <- rbind(cbind(fakedata.x2$t, fakedata.x2$t), fakedata.x2$o)
  +   fakedata.x2.Delta <- fitdistrplus:::Deltamatrix(fakedata.x2.x) # a single vector
  +   
  +   
  +   resfk <- fitdistrplus:::npsurvminimal(fakedata, w=1, maxit=100, tol=1e-6, verb=mytrace, pkg="stats") 
  +   
  +   #---------------------------------------------------------------------------------------
  +   # geometric example with interval censoring leading to 50 maximal intersection intervals
  +   
  +   LR <- matrix(1:100, ncol=2, byrow = TRUE)
  +   theop <- dgeom(1:NROW(LR)-1, prob=1/10)
  +   cnt <- round(1000*pgeom(1:NROW(LR)-1, prob=1/10, lower=FALSE))
  +   fakedata <- cbind(L=LR[,1], R=LR[,2], count=cnt)
  +   
  +   fakedata.x2 <- fitdistrplus:::icendata(fakedata, w=1)
  +   fakedata.x2.x <- rbind(cbind(fakedata.x2$t, fakedata.x2$t), fakedata.x2$o)
  +   fakedata.x2.Delta <- fitdistrplus:::Deltamatrix(fakedata.x2.x) #is diagonal
  +   
  +   
  +   resfk <- fitdistrplus:::npsurvminimal(fakedata, w=1, maxit=100, tol=1e-6, verb=mytrace, pkg="stats") 
  +   
  +   head(cbind(optimized=resfk$f$p, theoretical=theop))
  +   
  +   #----------------------------------------------------------------------------------------
  +   # multimodal example with interval censoring leading to 50 maximal intersection intervals
  +   
  +   LR <- matrix(1:100, ncol=2, byrow = TRUE)
  +   cnt <- round(450*dgeom(1:NROW(LR)-1, prob=1/10)
  +                +550*dbinom(1:NROW(LR)-1, size=NROW(LR), prob=4/5))
  +   theop <- cnt/sum(cnt)
  +   fakedata <- cbind(L=LR[,1], R=LR[,2], count=cnt)
  +   
  +   fakedata.x2 <- fitdistrplus:::icendata(fakedata, w=1)
  +   fakedata.x2.x <- rbind(cbind(fakedata.x2$t, fakedata.x2$t), fakedata.x2$o)
  +   fakedata.x2.Delta <- fitdistrplus:::Deltamatrix(fakedata.x2.x) #is diagonal
  +   
  +   
  +   resfk <- fitdistrplus:::npsurvminimal(fakedata, w=1, maxit=100, tol=1e-6, verb=mytrace, pkg="stats") 
  +   
  +   head(cbind(optimized=resfk$f$p, theoretical=theop[1:resfk$m]))
  +   
  +   
  +   #----------------------------------------------------------------------------------------
  +   # multimodal example with interval censoring leading to 43 maximal intersection intervals
  +   
  +   n <- 100
  +   set.seed(123)
  +   LR <- sample.int(n)
  +   LR <- cbind(LR, LR+sample.int(n)/10)
  +   cnt <- 1+round(450*dgeom(1:NROW(LR)-1, prob=1/10)
  +                  +550*dbinom(1:NROW(LR)-1, size=NROW(LR), prob=4/5))
  +   theop <- cnt/sum(cnt)
  +   fakedata <- cbind(L=LR[,1], R=LR[,2], count=cnt)
  +   
  +   fakedata.x2 <- fitdistrplus:::icendata(fakedata, w=1)
  +   fakedata.x2.x <- rbind(cbind(fakedata.x2$t, fakedata.x2$t), fakedata.x2$o)
  +   fakedata.x2.Delta <- fitdistrplus:::Deltamatrix(fakedata.x2.x)
  +   str(fakedata.x2.Delta)
  +   
  +   resfk <- fitdistrplus:::npsurvminimal(fakedata, w=1, maxit=100, tol=1e-6, verb=mytrace, pkg="stats") 
  +   
  +   
  +   # if(FALSE)
  +   # {
  +   #   require(npsurv)
  +   #   rescheck <- npsurv::npsurv(fakedata, w=1, maxit=100, tol=1e-6, verb=3)
  +   #   c(resfk$ll, rescheck$ll)
  +   # 
  +   #   cbind(resfk$f$left, rescheck$f$left)
  +   #   cbind(resfk$f$right, rescheck$f$right)
  +   #   cbind(resfk$f$p, rescheck$f$p)
  +   #   sum(abs(resfk$f$p- rescheck$f$p))
  +   # }
  +   
  +   
  +   #----------------------------------------------------------------------------------------
  +   # simulated example with interval censoring leading to 258 maximal intersection intervals
  +   
  +   set.seed(1232)
  +   ns <- 500
  +   ns <- 100
  +   r <- rnorm(ns)
  +   d8 <- data.frame(left = r, right = r)
  +   delta <- rlnorm(ns)
  +   icensored <- rbinom(ns, size = 1, prob = 0.2) 
  +   Lcensored <- rbinom(ns, size = 1, prob = 0.2*(1 - icensored))
  +   Rcensored <- rbinom(ns, size = 1, prob = 0.3*(1 - icensored)*(1 - Lcensored))
  +   # icensored +  Lcensored + Rcensored
  +   d8$left <- d8$left * (1 - Lcensored) + (-1000) * Lcensored
  +   d8$right <- d8$right * (1 - Rcensored) + (1000) * Rcensored
  +   d8$right <- d8$right + delta * icensored
  +   d8$right[d8$right == 1000] <- +Inf
  +   d8$left[d8$left == -1000] <- -Inf
  +   
  +   
  +   d8.x2 <- fitdistrplus:::icendata(d8, w=1)
  +   d8.x2.x <- rbind(cbind(d8.x2$t, d8.x2$t), d8.x2$o)
  +   d8.x2.Delta <- fitdistrplus:::Deltamatrix(d8.x2.x)
  +   str(d8.x2.Delta)
  +   
  +   system.time(
  +     resd8 <- fitdistrplus:::npsurvminimal(d8, w=1, maxit=100, tol=1e-6, verb=mytrace, pkg="stats") 
  +   )
  +   
  +   # if(FALSE)
  +   # {
  +   #   require(npsurv)
  +   #   d8bis <- d8
  +   #   d8bis$left[is.na(d8bis$left)] <- -Inf
  +   #   d8bis$right[is.na(d8bis$right)] <- Inf
  +   #   
  +   #   rescheck <- npsurv::npsurv(d8bis, w=1, maxit=100, tol=1e-6, verb=3)
  +   #   c(resd8$ll, rescheck$ll)
  +   #   
  +   #   cbind(resd8$f$left, rescheck$f$left)
  +   #   cbind(resd8$f$right, rescheck$f$right)
  +   #   head(cbind(resd8$f$p, rescheck$f$p))
  +   #   sum(abs(resd8$f$p- rescheck$f$p))
  +   # }
  +   
  +   
  +   #------------------------------------------------------
  +   # crash example with wrong interval censoring intervals
  +   
  +   set.seed(1232)
  +   ns <- 100
  +   r <- rnorm(ns)
  +   d8 <- data.frame(left = r, right = r)
  +   delta <- rlnorm(ns)
  +   icensored <- rbinom(ns, size = 1, prob = 0.2) 
  +   Lcensored <- rbinom(ns, size = 1, prob = 0.2*(1 - icensored))
  +   Rcensored <- rbinom(ns, size = 1, prob = 0.3*(1 - icensored)*(1 - Lcensored))
  +   # icensored +  Lcensored + Rcensored
  +   d8$left <- d8$left * (1 - Lcensored) + (-1000) * Lcensored
  +   d8$right <- d8$right * (1 - Rcensored) + (1000) * Rcensored
  +   d8$right <- d8$right + delta * icensored
  +   d8$right[d8$right == 1000] <- -Inf
  +   d8$left[d8$left == -1000] <- +Inf
  +   
  +   try(resd8 <- fitdistrplus:::npsurvminimal(d8, w=1, maxit=100, tol=1e-6, verb=2, pkg="stats"))
  +   
  + }
  > 
  > proc.time()
     user  system elapsed 
     1.46    0.09    1.54 
Running the tests in 'tests/t-util-testdensity.R' failed.
Complete output:
  > library(fitdistrplus)
  Loading required package: MASS
  Loading required package: survival
  > 
  > 
  > testdpqfun <- fitdistrplus:::testdpqfun
  > 
  > 
  > 
  > ##### first argument ##### 
  > #a data.frame of TRUE and ""
  > testdpqfun("exp", start=c(rate=1))
      ok txt
  1 TRUE    
  2 TRUE    
  3 TRUE    
  > #a data.frame with error messages
  > dEXP <- function(y, rate) dexp(x, rate)
  > pEXP <- function(y, rate) pexp(x, rate)
  > qEXP <- function(y, rate) qexp(x, rate)
  > testdpqfun("EXP", start=c(rate=1))
       ok                                                                    txt
  1 FALSE The dEXP function should have its first argument named: x as in base R
  2 FALSE The pEXP function should have its first argument named: q as in base R
  3 FALSE The qEXP function should have its first argument named: p as in base R
  > 
  > 
  > ##### existence ##### 
  > #a data.frame of TRUE and ""
  > testdpqfun("exp", start=c(rate=1))
      ok txt
  1 TRUE    
  2 TRUE    
  3 TRUE    
  > #a data.frame with error messages
  > testdpqfun("exp2", start=c(rate=1))
       ok                                txt
  1 FALSE The dexp2 function must be defined
  2 FALSE The pexp2 function must be defined
  3 FALSE The qexp2 function must be defined
  > 
  > ##### void vector ##### 
  > dexp2 <- function(x, rate)
  +   ifelse(length(x)==0, stop("zero input"), dexp(x,rate))
  > dexp3 <- function(x, rate)
  +   ifelse(length(x)==0, NA, dexp(x,rate))
  > #TRUE
  > testdpqfun("exp", "d", c(rate=1))
      ok txt
  1 TRUE    
  > #error message
  > testdpqfun("exp2", "d", c(rate=1))
       ok
  1 FALSE
                                                                                                        txt
  1 The dexp2 function should return a zero-length vector when input has length zero and not raise an error
  > #error message
  > testdpqfun("exp3", "d", c(rate=1))
       ok
  1 FALSE
                                                                                 txt
  1 The dexp3 function should return a zero-length vector when input has length zero
  > 
  > ##### inconsistent value ##### 
  > pexp2 <- function(q, rate)
  + {
  +   res <- pexp(q, rate)
  +   if(any(is.nan(res)))
  +     stop("NaN values")
  +   res
  + }
  > pexp3 <- function(q, rate)
  + {
  +   res <- pexp(q, rate)
  +   if(any(is.infinite(q)))
  +     stop("Inf values")
  +   res
  + }
  > 
  > #TRUE
  > testdpqfun("exp", "p", c(rate=1))
      ok txt
  1 TRUE    
  > #error message
  > testdpqfun("exp2", "p", c(rate=1))
       ok
  1 FALSE
                                                                                                                       txt
  1 The pexp2 function should return a vector of with NaN values when input has inconsistent values and not raise an error
  > #error message
  > testdpqfun("exp3", "p", c(rate=1))
       ok
  1 FALSE
                                                                                                                       txt
  1 The pexp3 function should return a vector of with NaN values when input has inconsistent values and not raise an error
  > 
  > ##### missing value ##### 
  > qexp2 <- function(p, rate)
  + {
  +   res <- qexp(p, rate)
  +   if(any(is.na(res)))
  +     stop("NA values")
  +   res
  + }
  > qexp3 <- function(p, rate)
  + {
  +   res <- qexp(p, rate)
  +   res[!is.na(res)]
  + }
  > 
  > #TRUE
  > testdpqfun("exp", "q", c(rate=1))
      ok txt
  1 TRUE    
  > #error message
  > testdpqfun("exp2", "q", c(rate=1))
       ok
  1 FALSE
                                                                                                                       txt
  1 The qexp2 function should return a vector of with NaN values when input has inconsistent values and not raise an error
  > #error message
  > testdpqfun("exp3", "q", c(rate=1))
       ok
  1 FALSE
                                                                                                                        txt
  1 The qexp3 function should return a vector of with NA values when input has missing values and not remove missing values
  > 
  > ##### inconsistent parameter ##### 
  > dnorm2 <- function(x, mean, sd)
  + {
  +   if(sd < 0)
  +     stop("negative param")
  +   else
  +     dnorm(x,mean,sd)
  + }
  > #TRUE
  > testdpqfun("norm", "d", c(mean=1, sd=1))
      ok txt
  1 TRUE    
  > #error message
  > testdpqfun("norm2", "d", c(mean=1, sd=1))
       ok
  1 FALSE
                                                                                                                            txt
  1 The dnorm2 function should return a vector of with NaN values when input has inconsistent parameters and not raise an error
  > 
  > ##### inconsistent name ##### 
  > dnorm2 <- function(x, mean=0, sd=1, ...)
  +     dnorm(x, mean, sd)
  > 
  > dnorm3 <- dnorm2
  > pnorm3 <- pnorm
  > qnorm3 <- qnorm
  > 
  > #TRUE
  > testdpqfun("norm", "d", c(mean=1, sd=1))
      ok txt
  1 TRUE    
  > #error message
  > testdpqfun("norm2", "d", c(mean=1, sd=1))
       ok
  1 FALSE
                                                                           txt
  1 The dnorm2 function should raise an error when names are incorrectly named
  > 
  > 
  > #a data.frame with error messages
  > testdpqfun("norm", c("d", "p", "q"), c(mean=1, sd=1))
      ok txt
  1 TRUE    
  2 TRUE    
  3 TRUE    
  > testdpqfun("norm2", c("d", "p", "q"), c(mean=1, sd=1))
       ok
  1 FALSE
  2 FALSE
  3 FALSE
                                                                           txt
  1 The dnorm2 function should raise an error when names are incorrectly named
  2                                        The pnorm2 function must be defined
  3                                        The qnorm2 function must be defined
  > testdpqfun("norm3", c("d", "p", "q"), c(mean=1, sd=1))
       ok
  1 FALSE
  2  TRUE
  3  TRUE
                                                                           txt
  1 The dnorm3 function should raise an error when names are incorrectly named
  2                                                                           
  3                                                                           
  > 
  > x <- rnorm(100)
  > fitdist(x, "norm") #ok
  Fitting of the distribution ' norm ' by maximum likelihood 
  Parameters:
         estimate Std. Error
  mean 0.06737797 0.08815511
  sd   0.88155109 0.06233471
  > fitdist(x, "norm2", start=list(mean=1, sd=1)) #pnorm2 not defined
  Fitting of the distribution ' norm2 ' by maximum likelihood 
  Parameters:
         estimate Std. Error
  mean 0.06729994 0.08814857
  sd   0.88148565 0.06232315
  Warning messages:
  1: In fitdist(x, "norm2", start = list(mean = 1, sd = 1)) :
    The dnorm2 function should raise an error when names are incorrectly named
  2: In fitdist(x, "norm2", start = list(mean = 1, sd = 1)) :
    The pnorm2 function must be defined
  > fitdist(x, "norm3", start=list(mean=1, sd=1)) #The dnorm3 function should return raise an error when names are incorrectly named
  Fitting of the distribution ' norm3 ' by maximum likelihood 
  Parameters:
         estimate Std. Error
  mean 0.06729994 0.08814857
  sd   0.88148565 0.06232315
  Warning message:
  In fitdist(x, "norm3", start = list(mean = 1, sd = 1)) :
    The dnorm3 function should raise an error when names are incorrectly named
  > 
  > 
  > proc.time()
     user  system elapsed 
     1.39    0.17    1.59 
Running the tests in 'tests/t-weird-ppcomp-cens.R' failed.
Complete output:
  > library(fitdistrplus)
  Loading required package: MASS
  Loading required package: survival
  > 
  > visualize <- FALSE # TRUE for manual tests with visualization of results
  > set.seed(123)
  > n <- 20
  > if (visualize)
  + {
  +   # (1) test ppcomp/ppcompcens on a good example
  +   #
  +   
  +   
  +   x <- rlnorm(n, 0, 1)
  +   dx <- data.frame(left=x, right=x)
  +   dx$right[1:(n/2)*2] <- NA
  +   dx$left[2:(n/4)*4-1] <- NA
  +   
  +   f1 <- fitdist(x, "lnorm")
  +   f1c <- fitdistcens(dx, "lnorm")
  +   
  +   f3 <- fitdist(x, "lnorm", fix.arg=list(sdlog=1))
  +   f3c <- fitdistcens(dx, "lnorm", fix.arg=list(sdlog=1))
  +   
  +   
  +   
  +   par(mfrow=1:2, mar=c(4,4,2,1))
  +   ppcomp(f1)
  +   ppcompcens(f1c)
  +   par(mfrow = c(1,1))
  +   if (requireNamespace ("ggplot2", quietly = TRUE)) {
  +     ppcomp(f1, plotstyle = "ggplot")
  +     ppcompcens(f1c, plotstyle = "ggplot")
  +   }
  +   
  +   #test log-scale
  +   par(mfrow=1:2, mar=c(4,4,2,1))
  +   ppcomp(f1, xlogscale = TRUE, ylogscale = TRUE)
  +   ppcompcens(f1c, xlogscale = TRUE, ylogscale = TRUE, xlim=c(.1, 1), ylim=c(.1, 1))
  +   par(mfrow = c(1,1))
  +   if (requireNamespace ("ggplot2", quietly = TRUE)) {
  +     ppcomp(f1, xlogscale = TRUE, ylogscale = TRUE, plotstyle = "ggplot")
  +     ppcompcens(f1c, xlogscale = TRUE, ylogscale = TRUE, xlim=c(.1, 1), ylim=c(.1, 1), plotstyle = "ggplot")
  +   }
  +   
  +   
  +   # (2) test ppcomp/ppcompcens on a weird example
  +   #
  +   
  +   f2 <- fitdist(x, "unif")
  +   f2c <- fitdistcens(dx, "unif")
  +   
  +   par(mfrow=1:2, mar=c(4,4,2,1))
  +   ppcomp(list(f1, f2, f3))
  +   ppcompcens(list(f1c, f2c, f3c))
  +   par(mfrow = c(1,1))
  +   if (requireNamespace ("ggplot2", quietly = TRUE)) {
  +     ppcomp(list(f1, f2, f3), plotstyle = "ggplot")
  +     ppcompcens(list(f1c, f2c, f3c), plotstyle = "ggplot")
  +   }
  +   
  +   #test log-scale
  +   par(mfrow=1:2, mar=c(4,4,2,1))
  +   ppcomp(list(f1, f2, f3), xlogscale = TRUE, ylogscale = TRUE, xlim=c(.1, 1), ylim=c(.1, 1))
  +   ppcompcens(list(f1c, f2c, f3c), xlogscale = TRUE, ylogscale = TRUE, xlim=c(.1, 1), ylim=c(.1, 1))
  +   par(mfrow = c(1,1))
  +   if (requireNamespace ("ggplot2", quietly = TRUE)) {
  +     ppcomp(list(f1, f2, f3), xlogscale = TRUE, ylogscale = TRUE, xlim=c(.1, 1), ylim=c(.1, 1), plotstyle = "ggplot")
  +     ppcompcens(list(f1c, f2c, f3c), xlogscale = TRUE, ylogscale = TRUE, xlim=c(.1, 1), ylim=c(.1, 1), plotstyle = "ggplot")
  +   }
  +   
  +   #test y noise
  +   par(mfrow=1:2, mar=c(4,4,2,1))
  +   ppcomp(list(f1, f2, f3))
  +   ppcomp(list(f1, f2, f3), ynoise=FALSE)
  +   par(mfrow = c(1,1))
  +   if (requireNamespace ("ggplot2", quietly = TRUE)) {
  +     ppcomp(list(f1, f2, f3), plotstyle = "ggplot")
  +     ppcomp(list(f1, f2, f3), ynoise=FALSE, plotstyle = "ggplot")
  +   }
  +   
  +   par(mfrow=1:2, mar=c(4,4,2,1))
  +   ppcompcens(list(f1c, f2c, f3c))
  +   ppcompcens(list(f1c, f2c, f3c), ynoise=FALSE)
  +   par(mfrow = c(1,1))
  +   if (requireNamespace ("ggplot2", quietly = TRUE)) {
  +     ppcompcens(list(f1c, f2c, f3c), ynoise=FALSE, plotstyle = "ggplot")
  +   }
  +   
  +   #test log-scale y-noise
  +   par(mfrow=1:2, mar=c(4,4,2,1))
  +   ppcomp(list(f1, f2, f3), xlogscale = TRUE, ylogscale = TRUE, xlim=c(.1, 1), ylim=c(.1, 1))
  +   ppcomp(list(f1, f2, f3), xlogscale = TRUE, ylogscale = TRUE, xlim=c(.1, 1), ylim=c(.1, 1), ynoise=FALSE)
  +   par(mfrow = c(1,1))
  +   if (requireNamespace ("ggplot2", quietly = TRUE)) {
  +     ppcomp(list(f1, f2, f3), xlogscale = TRUE, ylogscale = TRUE, xlim=c(.1, 1), ylim=c(.1, 1), plotstyle = "ggplot")
  +     ppcomp(list(f1, f2, f3), xlogscale = TRUE, ylogscale = TRUE, xlim=c(.1, 1), ylim=c(.1, 1), plotstyle = "ggplot")
  +   }
  +   
  +   par(mfrow=1:2, mar=c(4,4,2,1))
  +   ppcompcens(list(f1c, f2c), xlogscale = TRUE, ylogscale = TRUE, xlim=c(.1, 1), ylim=c(.1, 1))
  +   ppcompcens(list(f1c, f2c), xlogscale = TRUE, ylogscale = TRUE, xlim=c(.1, 1), ylim=c(.1, 1), ynoise=FALSE)
  +   par(mfrow = c(1,1))
  +   if (requireNamespace ("ggplot2", quietly = TRUE)) {
  +     ppcompcens(list(f1c, f2c), xlogscale = TRUE, ylogscale = TRUE, xlim=c(.1, 1), ylim=c(.1, 1), plotstyle = "ggplot")
  +     ppcompcens(list(f1c, f2c), xlogscale = TRUE, ylogscale = TRUE, xlim=c(.1, 1), ylim=c(.1, 1), plotstyle = "ggplot")
  +   }
  +   
  + }
  > 
  > proc.time()
     user  system elapsed 
     1.18    0.20    1.37 
Running the tests in 'tests/t-weird-qqcomp-cens.R' failed.
Complete output:
  > library(fitdistrplus)
  Loading required package: MASS
  Loading required package: survival
  > set.seed(123)
  > n <- 20
  > visualize <- FALSE # TRUE for manual tests with visualization of results
  > if (visualize)
  + {
  +   # (1) test qqcomp/qqcompcens on a good example
  +   #
  +   
  +   
  +   x <- rlnorm(n, 0, 1)
  +   dx <- data.frame(left=x, right=x)
  +   dx$right[1:(n/2)*2] <- NA
  +   dx$left[2:(n/4)*4-1] <- NA
  +   
  +   f1 <- fitdist(x, "lnorm")
  +   f1c <- fitdistcens(dx, "lnorm")
  +   
  +   f3 <- fitdist(x, "lnorm", fix.arg=list(sdlog=1))
  +   f3c <- fitdistcens(dx, "lnorm", fix.arg=list(sdlog=1))
  +   
  +   
  +   
  +   par(mfrow=1:2)
  +   qqcomp(f1)
  +   qqcompcens(f1c)
  +   par(mfrow = c(1,1))
  +   if(requireNamespace ("ggplot2", quietly = TRUE)) {
  +     qqcomp(f1, plotstyle = "ggplot")
  +     qqcompcens(f1c, plotstyle = "ggplot")
  +   }
  +   
  +   #test log-scale
  +   par(mfrow=1:2, mar=c(4,4,2,1))
  +   qqcomp(f1, xlogscale = TRUE, ylogscale = TRUE)
  +   qqcompcens(f1c, xlogscale = TRUE, ylogscale = TRUE)
  +   par(mfrow = c(1,1))
  +   if (requireNamespace ("ggplot2", quietly = TRUE)) {
  +     qqcomp(f1, xlogscale = TRUE, ylogscale = TRUE, plotstyle = "ggplot")
  +     qqcompcens(f1c, xlogscale = TRUE, ylogscale = TRUE, plotstyle = "ggplot")
  +   }
  +   
  +   # (2) test qqcomp/qqcompcens on a weird example
  +   #
  +   
  +   
  +   f2 <- fitdist(x, "unif")
  +   f2c <- fitdistcens(dx, "unif")
  +   
  +   par(mfrow=1:2, mar=c(4,4,2,1))
  +   qqcomp(list(f1, f2, f3))
  +   qqcompcens(list(f1c, f2c, f3c))
  +   par(mfrow = c(1,1))
  +   if (requireNamespace ("ggplot2", quietly = TRUE)) {
  +     qqcomp(list(f1, f2, f3), plotstyle = "ggplot")
  +     qqcompcens(list(f1c, f2c, f3c), plotstyle = "ggplot")
  +   }
  +   
  +   #test log-scale
  +   par(mfrow=1:2, mar=c(4,4,2,1))
  +   qqcomp(list(f1, f2, f3), xlogscale = TRUE, ylogscale = TRUE)
  +   qqcompcens(list(f1c, f2c, f3c), xlogscale = TRUE, ylogscale = TRUE)
  +   par(mfrow = c(1,1))
  +   if (requireNamespace ("ggplot2", quietly = TRUE)) {
  +     qqcomp(list(f1, f2, f3), xlogscale = TRUE, ylogscale = TRUE, plotstyle = "ggplot")
  +     qqcompcens(list(f1c, f2c, f3c), xlogscale = TRUE, ylogscale = TRUE, plotstyle = "ggplot")
  +   }
  +   
  +   #test y noise
  +   par(mfrow=1:2, mar=c(4,4,2,1))
  +   qqcomp(list(f1, f2, f3))
  +   qqcomp(list(f1, f2, f3), ynoise=FALSE)
  +   par(mfrow = c(1,1))
  +   if (requireNamespace ("ggplot2", quietly = TRUE)) {
  +     qqcomp(list(f1, f2, f3), ynoise=FALSE, plotstyle = "ggplot")
  +   }
  +   
  +   par(mfrow=1:2, mar=c(4,4,2,1))
  +   qqcompcens(list(f1c, f2c, f3c))
  +   qqcompcens(list(f1c, f2c, f3c), ynoise=FALSE)
  +   par(mfrow = c(1,1))
  +   if (requireNamespace ("ggplot2", quietly = TRUE)) {
  +     qqcompcens(list(f1c, f2c, f3c), plotstyle = "ggplot")
  +   }
  +   
  +   #test log-scale y-noise
  +   par(mfrow=1:2, mar=c(4,4,2,1))
  +   qqcomp(list(f1, f2, f3), xlogscale = TRUE, ylogscale = TRUE)
  +   qqcomp(list(f1, f2, f3), xlogscale = TRUE, ylogscale = TRUE, ynoise=FALSE)
  +   par(mfrow = c(1,1))
  +   if (requireNamespace ("ggplot2", quietly = TRUE)) {
  +     qqcomp(list(f1, f2, f3), xlogscale = TRUE, ylogscale = TRUE, ynoise=FALSE, plotstyle = "ggplot")
  +   }
  +   
  +   par(mfrow=1:2, mar=c(4,4,2,1))
  +   qqcompcens(list(f1c, f2c), xlogscale = TRUE, ylogscale = TRUE)
  +   qqcompcens(list(f1c, f2c), xlogscale = TRUE, ylogscale = TRUE, ynoise=FALSE)
  +   par(mfrow = c(1,1))
  +   if (requireNamespace ("ggplot2", quietly = TRUE)) {
  +     qqcompcens(list(f1c, f2c), xlogscale = TRUE, ylogscale = TRUE, plotstyle = "ggplot")
  +   }
  +   
  + }
  > 
  > proc.time()
     user  system elapsed 
     1.17    0.20    1.31