R Under development (unstable) (2025-01-06 r87534 ucrt) -- "Unsuffered Consequences"
Copyright (C) 2025 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
>
> nsample <- 500
> nsample <- 10
> visualize <- FALSE # TRUE for manual tests with visualization of results
> set.seed(1234)
>
>
> # (1) check input data
> #
> data(smokedfish)
> fitsf <- fitdistcens(smokedfish, "lnorm")
> summary(fitsf)
Fitting of the distribution ' lnorm ' By maximum likelihood on censored data
Parameters
estimate Std. Error
meanlog -3.627606 0.4637122
sdlog 3.544570 0.4876610
Loglikelihood: -90.65154 AIC: 185.3031 BIC: 190.5725
Correlation matrix:
meanlog sdlog
meanlog 1.0000000 -0.4325873
sdlog -0.4325873 1.0000000
>
> #fitted coef is -3.63 3.54, fitted loglik is -90.7
>
> smokedfish2 <- tail(smokedfish, 20)
>
> try(fitdistcens(rbind(smokedfish2, c(NA, NA)) , "lnorm"))
Error in checkCensoredDataFrameNAInfNan(censdata) :
censdata contain two NA values on the same line.
>
> try(fitdistcens(rbind(smokedfish2, c(3, Inf)) , "lnorm"))
Error in checkCensoredDataFrameNAInfNan(censdata) :
censdata contain Inf (infinite) values.
>
> try(fitdistcens(rbind(smokedfish2, c(NaN, 3)) , "lnorm"))
Error in checkCensoredDataFrameNAInfNan(censdata) :
censdata contain NaN (not a numeric) values.
>
> try(fitdistcens(rbind(smokedfish2, c(5, 3)) , "lnorm"))
Error in checkCensoredDataFrameNAInfNan(censdata) :
each censdata$left value must be less or equal to the corresponding censdata$right value
>
> fitsf <- fitdistcens(smokedfish, "lnorm", calcvcov = FALSE)
>
>
>
> # (6) custom optimisation function - example with the genetic algorithm
> #
> data(fluazinam)
> log10EC50 <- log10(fluazinam)
>
> summary(fitdistcens(log10EC50, "logis"))
Fitting of the distribution ' logis ' By maximum likelihood on censored data
Parameters
estimate Std. Error
location 2.1518291 0.3222830
scale 0.6910423 0.1745231
Loglikelihood: -20.55391 AIC: 45.10781 BIC: 46.38593
Correlation matrix:
location scale
location 1.00000000 0.05097494
scale 0.05097494 1.00000000
>
> #fitted coef is 2.152 0.691 , fitted loglik is -20.6
>
> #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=1, 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.11, Build Date: 2024-10-03)
## 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.
##
Tue Jan 7 14:37:36 2025
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 1.471610e+00
1 1.468136e+00
'wait.generations' limit reached.
No significant improvement in 10 generations.
Solution Fitness Value: 1.468136e+00
Parameters at the Solution (parameter, gradient):
X[ 1] : 2.151909e+00 G[ 1] : -1.092935e-06
X[ 2] : 6.909667e-01 G[ 2] : 1.301687e-06
Solution Found Generation 1
Number of Generations Run 12
Tue Jan 7 14:37:38 2025
Total run time : 0 hours 0 minutes and 2 seconds
>
> summary(fit.with.genoud)
Fitting of the distribution ' logis ' By maximum likelihood on censored data
Parameters
estimate Std. Error
location 2.1519089 0.3222567
scale 0.6909667 0.1744837
Loglikelihood: -20.55391 AIC: 45.10781 BIC: 46.38593
Correlation matrix:
location scale
location 1.0000000 0.0510634
scale 0.0510634 1.0000000
>
>
> # (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)
+
+ #fitted coef is 5, fitted loglik is 609
+ }
>
>
> # (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.840868
Fixed parameters:
value
shape 1.5
> f1$fix.arg
$shape
[1] 1.5
>
> #fitted coef is 1.5, fitted loglik is -14.7
>
> 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.840868
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
>
>
> #fitted coef is 3.385 0.496, fitted loglik is -139
>
>
> # (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.412509
rate 19.140930
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.31 0.17 2.46