R Under development (unstable) (2024-10-16 r87241 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(mbbefd)
Loading required package: fitdistrplus
Loading required package: MASS
Loading required package: survival
Loading required package: alabama
Loading required package: numDeriv
Loading required package: Rcpp
Package:  mbbefd
Version:  0.8.12
Date:     
BugReport: https://github.com/spedygiorgio/mbbefd/issues


> library(fitdistrplus)
> 
> #____________________________________________________________
> #gbeta
> n <- 1e3
> nboot <- 100
> nboot <- 10
> set.seed(12345)
> x <- rgbeta(n, 2, 2, 5/2)
> 
> initpar <- list(shape0=2, shape1=2, shape2=5/2)
> 
> if(FALSE)
+ {
+ #____________________________________________________________
+ #test all computation methods
+ ctr <- list(trace=0, REPORT=1, maxit=1000)
+ reslist <- NULL
+ for(meth in c("BFGS", "Nelder", "CG")) #CG with FR update
+ {
+   nograd$time <- system.time(nograd <- mledist(x, dist="gbeta", optim.method=meth, 
+                 control=ctr, start=initpar))[3]
+   reslist <- c(reslist, list(nograd))
+ }
+ for(type in 2:3) #CG with PR or BS updates
+ {
+   nograd$time <- system.time(nograd <- mledist(x, dist="gbeta", optim.method="CG", 
+                 control=c(ctr, type=type), start=initpar))[3] 
+   reslist <- c(reslist, list(nograd)) 
+ }
+ fullname <- c("BFGS", "NM", paste("CG", c("FR", "PR", "BS")))
+ names(reslist) <- fullname
+ 
+ dgbeta2 <- function(x, shape0, shape1, shape2, log=FALSE)
+   dgbeta(x, exp(shape0), exp(shape1), exp(shape2), log=log)
+ 
+ initpar2 <- lapply(initpar, log)
+ 
+ for(meth in c("BFGS", "Nelder", "CG")) #CG with FR update
+ {
+   nograd$time <- system.time(nograd <- mledist(x, dist="gbeta2", optim.method="BFGS", 
+                           control=ctr, start=initpar2))[3]
+   nograd$estimate <- exp(nograd$estimate)
+   reslist <- c(reslist, list(nograd))
+ }
+ for(type in 2:3) #CG with PR or BS updates
+ {
+   nograd$time <- system.time(nograd <- mledist(x, dist="gbeta2", optim.method="CG", 
+                         control=c(ctr, type=type), start=initpar2))[3] 
+   nograd$estimate <- exp(nograd$estimate)
+   reslist <- c(reslist, list(nograd)) 
+ }
+ names(reslist)[(length(fullname)+1):length(reslist)] <- paste("exp.", fullname)
+ 
+ getval <- function(x)
+   c(x$estimate, loglik=x$loglik, x$counts, x$time)
+ 
+ resNM <- sapply(reslist[grep("NM", names(reslist))], getval)
+ resCG <- sapply(reslist[grep("CG", names(reslist))], getval)
+ resBFGS <- sapply(reslist[grep("BFGS", names(reslist))], getval)  
+ rownames(resNM) <- rownames(resCG) <- rownames(resBFGS) <- c(paste("fitted", c("shape0", "shape1", "shape2")), "fitted loglik", "func. eval. nb.", "grad. eval. nb.", "time (sec)")
+ 
+ 
+ #____________________________________________________________
+ #empirical check of the log-likelihood computation
+ head(cbind(do.call("dgbeta", c(list(x), initpar, NULL, log=TRUE)), 
+ log(do.call("dgbeta", c(list(x), initpar, NULL) ) )
+ ))
+ colSums(cbind(do.call("dgbeta", c(list(x), initpar, NULL, log=TRUE)), 
+               log(do.call("dgbeta", c(list(x), initpar, NULL) ) )
+ ))
+ 
+ #____________________________________________________________
+ #test with starting values equal theoretical values
+ 
+ f1 <- fitDR(x, "oigbeta", "mle")
+ 
+ f1 <- fitdist(x, "gbeta", method="mle", start=initpar) # , control=list(trace=3, REPORT=1)) 
+ summary(f1)
+ cdfcomp(f1, do.points=FALSE, ylogscale = TRUE)
+ lines(0:100/100, pgbeta(0:100/100, 2, 2, 5/2), col="green")
+ 
+ denscomp(f1)
+ lines(0:100/100, dgbeta(0:100/100, 2, 2, 5/2), col="green")
+ 
+ #____________________________________________________________
+ #look at the log-likelihood function around the estimated value
+ par(mfrow=c(1,3))
+ llsurface(plot.min=c(0.1, 0.1), plot.max=c(5, 4), nlevels=20,
+                          plot.arg=c("shape1", "shape2"), fix.arg=as.list(f1$estimate[1]),
+                          plot.np=50, obs=x, distr="gbeta", plot.type="contour")
+ points(f1$estimate["shape1"], f1$estimate["shape2"], pch="+", col="red")
+ points(2, 5/2, pch="x", col="green")
+ llsurface(plot.min=c(0.1, 0.1), plot.max=c(6, 4), nlevels=20,
+        plot.arg=c("shape0", "shape2"), fix.arg=as.list(f1$estimate[2]),
+        plot.np=50, obs=x, distr="gbeta", plot.type="contour")
+ points(f1$estimate["shape0"], f1$estimate["shape2"], pch="+", col="red")
+ points(2, 5/2, pch="x", col="green")
+ llsurface(plot.min=c(0.1, 0.1), plot.max=c(5, 6), nlevels=20,
+        plot.arg=c("shape1", "shape0"), fix.arg=as.list(f1$estimate[3]),
+        plot.np=50, obs=x, distr="gbeta", plot.type="contour")
+ points(f1$estimate["shape1"], f1$estimate["shape0"], pch="+", col="red")
+ points(2, 2, pch="x", col="green")
+ 
+ 
+ par(mfrow=c(1,3))
+ llcurve(plot.min=0.1, plot.max=5, plot.arg="shape0", fix.arg=as.list(f1$estimate[-1]), plot.np=50, 
+         obs=x, distr="gbeta", enhance=FALSE)
+ abline(v=c(2, f1$estimate["shape0"]), col=c("green", "red"))
+ llcurve(plot.min=0.1, plot.max=4, plot.arg="shape1", fix.arg=as.list(f1$estimate[-2]), plot.np=50, 
+                       obs=x, distr="gbeta", enhance=FALSE)
+ abline(v=c(2, f1$estimate["shape1"]), col=c("green", "red"))
+ llcurve(plot.min=0.1, plot.max=4, plot.arg="shape2", fix.arg=as.list(f1$estimate[-3]), plot.np=50, 
+                       obs=x, distr="gbeta", enhance=FALSE)
+ abline(v=c(5/2, f1$estimate["shape2"]), col=c("green", "red"))
+ 
+ 
+ #bootstrap
+ b1 <- bootdist(f1, niter=nboot, silent=TRUE)
+ summary(b1)
+ 
+ plot(b1, enhance=TRUE, trueval=c(2, 2, 5/2))
+ }
> 
> #____________________________________________________________
> #init value
> 
> 
> 
> 
> s00 <- optimize(function(z) 
+   (mbbefd:::Theil.emp(x) - mbbefd:::Theil.theo.shape0(z, obs=x))^2, lower=0.01, upper=20)$minimum
> initpar1 <- c(list(shape0=1), as.list(fitdist(x, "beta", method="mme")$estimate))
Warning messages:
1: In cov2cor(varcovar) :
  diag(V) had non-positive or NA entries; the non-finite result may be dubious
2: In sqrt(diag(varcovar)) : NaNs produced
> initpar2 <- c(list(shape0=s00), as.list(fitdist(x^s00, "beta", method="mme")$estimate))
> 
> 
> 
> 
> #____________________________________________________________
> 
> fitdist(x, "gbeta", method="mle", start=initpar1)
Fitting of the distribution ' gbeta ' by maximum likelihood 
Parameters:
        estimate Std. Error
shape0 4.2327147  1.4141502
shape1 0.8631648  0.3481856
shape2 2.9658930  0.2832509
> fitdist(x, "gbeta", method="mle", start=initpar2)
Fitting of the distribution ' gbeta ' by maximum likelihood 
Parameters:
        estimate Std. Error
shape0 4.2345065  1.4087431
shape1 0.8626725  0.3465092
shape2 2.9659882  0.2824930
> 
> f2 <- fitdist(x, "gbeta", method="mle", start=initpar1)
> summary(f2)
Fitting of the distribution ' gbeta ' by maximum likelihood 
Parameters : 
        estimate Std. Error
shape0 4.2327147  1.4141502
shape1 0.8631648  0.3481856
shape2 2.9658930  0.2832509
Loglikelihood:  407.0968   AIC:  -808.1936   BIC:  -793.4703 
Correlation matrix:
           shape0    shape1     shape2
shape0  1.0000000 -0.995292  0.8699008
shape1 -0.9952920  1.000000 -0.8320690
shape2  0.8699008 -0.832069  1.0000000

> cdfcomp(f2, do.points=FALSE, ylogscale = TRUE)
> lines(0:100/100, pgbeta(0:100/100, 2, 2, 5/2), col="green")
> 
> proc.time()
   user  system elapsed 
   2.42    0.25    2.65