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. > 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 3.26 0.45 3.64