R Under development (unstable) (2024-07-10 r86888 ucrt) -- "Unsuffered Consequences" Copyright (C) 2024 The R Foundation for Statistical Computing Platform: x86_64-w64-mingw32/x64 R is free software and comes with ABSOLUTELY NO WARRANTY. You are welcome to redistribute it under certain conditions. Type 'license()' or 'licence()' for distribution details. R is a collaborative project with many contributors. Type 'contributors()' for more information and 'citation()' on how to cite R or R packages in publications. Type 'demo()' for some demos, 'help()' for on-line help, or 'help.start()' for an HTML browser interface to help. Type 'q()' to quit R. > library(fitdistrplus) Loading required package: MASS Loading required package: survival > nbboot <- 101 > nbboot <- 11 > nsample <- 10 > visualize <- FALSE # TRUE for manual tests with visualization of results > > # (1) Fit of a normal distribution to fluazinam data in log10 > # followed by nonparametric bootstrap > # > data(fluazinam) > (d1 <-log10(fluazinam)) left right 1 0.5797836 0.5797836 2 1.5263393 1.5263393 3 1.9395193 1.9395193 4 3.2304489 NA 5 2.8061800 2.8061800 6 3.0625820 NA 7 2.0530784 2.0530784 8 2.1105897 2.1105897 9 2.7678976 2.7678976 10 3.2685780 NA 11 0.2041200 0.2041200 12 0.6812412 0.6812412 13 1.9138139 1.9138139 14 2.1903317 2.1903317 > f1 <- fitdistcens(d1, "norm") > b1 <- bootdistcens(f1, niter = nbboot, silent=TRUE) > b1 <- bootdistcens(f1, niter = nbboot, silent=FALSE) > > > b1 Parameter values obtained with nonparametric bootstrap mean sd 1 2.530599 1.4083597 2 1.883131 1.0327915 3 1.760712 1.0674224 4 2.100774 1.2287994 5 1.835305 1.2607802 6 1.744712 1.0550995 7 2.218063 1.3671278 8 2.019242 1.3127166 9 2.434384 0.8926819 10 2.447645 1.4436840 11 2.415643 1.4894924 > summary(b1) Nonparametric bootstrap medians and 95% percentile CI Median 2.5% 97.5% mean 2.100774 1.7487121 2.50986 sd 1.260780 0.9277093 1.47804 > plot(b1) > > # (3) Estimation of the standard deviation of a normal distribution > # by maximum likelihood with the mean fixed at 0.1 using the argument fix.arg > # followed by nonparametric bootstrap > # > f1b <- fitdistcens(d1, "norm", start=list(sd=1.5), fix.arg=list(mean=0.1)) > b1b <- bootdistcens(f1b, niter=nbboot) > > summary(b1b) Nonparametric bootstrap medians and 95% percentile CI Median 2.5% 97.5% 2.561314 1.977308 3.781893 > plot(b1b) > > # (4) Comparison of fitdist and fitdistcens and bootdist and bootdistcens > # for non censored data > 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) > fx1<-fitdist(x1,"norm",method="mle") > cx1<-bootdist(fx1,bootmethod="nonparam", niter=nbboot) > xx1<-data.frame(left=x1,right=x1) > fxx1<-fitdistcens(xx1,"norm") > summary(fx1) Fitting of the distribution ' norm ' by maximum likelihood Parameters : estimate Std. Error mean 10.411111 1.118886 sd 4.747033 0.791172 Loglikelihood: -53.57625 AIC: 111.1525 BIC: 112.9332 Correlation matrix: mean sd mean 1 0 sd 0 1 > summary(fxx1) Fitting of the distribution ' norm ' By maximum likelihood on censored data Parameters estimate Std. Error mean 10.411111 1.118886 sd 4.747033 0.791172 Loglikelihood: -53.57625 AIC: 111.1525 BIC: 112.9332 Correlation matrix: mean sd mean 1 0 sd 0 1 > cdfcomp(fx1) > cdfcompcens(fxx1) > cxx1<-bootdistcens(fxx1, niter=nbboot) > summary(cx1) Nonparametric bootstrap medians and 95% percentile CI Median 2.5% 97.5% mean 10.461011 7.784975 11.162961 sd 4.546917 3.719441 5.579375 > summary(cxx1) Nonparametric bootstrap medians and 95% percentile CI Median 2.5% 97.5% mean 9.883519 8.640787 11.341877 sd 4.303428 3.123505 5.409085 > > > > # (5) fixing parameters > # > set.seed(1234) > x <- rexp(nsample, 5) > x <- data.frame(left=x, right=x+.1) > > f1 <- fitdistcens(x, "gamma", fix.arg=list(shape=1.5)) > b1 <- bootdistcens(f1, niter=nbboot) > plot(b1) > > f1 <- fitdistcens(x, "gamma", fix.arg=function(x) list(shape=1.5)) > b1 <- bootdistcens(f1, niter=nbboot) > plot(b1) > > # (6) efficiency of parallel operation > if (visualize) # too long to run on CRAN and forbidden due to parallel computing + { + niter <- 5001 + data(fluazinam) + d1 <-log10(fluazinam) + f1 <- fitdistcens(d1, "norm") + + for (cli in 1:4) + { + print(cli) + ptm <- proc.time() + print(summary(bootdistcens(f1, niter = niter, parallel = "snow", ncpus = cli))) + print(proc.time() - ptm) + } + # not available on Windows + for (cli in 1:4) + { + print(cli) + ptm <- proc.time() + print(summary(bootdistcens(f1, niter = niter, parallel = "multicore", ncpus = cli))) + print(proc.time() - ptm) + } + } > > # (5) with weights (not yet available, test of error message) > # > 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 > summary(bootdistcens(fa, niter = nbboot)) Nonparametric bootstrap medians and 95% percentile CI Median 2.5% 97.5% meanlog 3.3771835 3.3310508 3.5324526 sdlog 0.4821702 0.4609298 0.5453199 > > > proc.time() user system elapsed 2.25 0.37 2.57