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
> nsample <- 1000
> nsample <- 10
> 
> 
> #(1) beta distribution
> #
> 
> x <- rbeta(nsample, 3, 3/4)
> 
> llsurface(data = x, distr = "beta", plot.arg=c("shape1", "shape2"),
+           min.arg=c(0.1, 0.1), max.arg=c(7, 3))
> llsurface(data = x, distr = "beta", plot.arg=c("shape1", "shape2"),
+           min.arg=c(0.1, 0.1), max.arg=c(7, 3), back.col = FALSE )
> points(3, 3/4, pch="+", col="red")
> 
> 
> llcurve(data = x, distr = "beta", plot.arg = "shape1", min.arg = 0.1, max.arg = 7, 
+         fix.arg = list(shape2 = 3/4), lseq=100, col = "blue")
> llcurve(data = x, distr = "beta", plot.arg = "shape2", min.arg = 0.1, max.arg = 7, 
+         fix.arg = list(shape1 = 3), lseq=100, col = "red")
> 
> 
> 
> 
> #test
> 
> 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)
 [1]  0.6754801  0.8966304  0.5610926  0.8940287  0.3066702  0.6169070
 [7]  0.5754737  0.7896227  0.9212068  0.0179939 -0.8101866 -2.3404130
[13] -0.5159069 -2.2941043 -0.1289104 -0.6445948 -0.5468519 -1.2926766
[19] -2.9453115  0.1161824
> 
> 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)))
[1]   6.255106 -11.402774
> grlnL(c(3, 4), x)
[1] -6.255106 11.402774
> grlnL(c(3, 3/4), x)
[1]  0.6473635 -5.1145425
> 
> 
> ctr <- list(trace=0, REPORT=1, maxit=1000)
> 
> bfgs_gr <- mledist(x, dist="beta", optim.method="BFGS", gr=grlnL, control=ctr)
> bfgs <- mledist(x, dist="beta", optim.method="BFGS", control=ctr)
> 
> 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)
> 
> getval <- function(x)
+   c(x$estimate, loglik=x$loglik, x$counts)
> 
> cbind(NM=getval(nm), NMgrad=getval(nm_gr), CG=getval(cg), 
+       CGgrad=getval(cg_gr), BFGS=getval(bfgs), BFGSgrad=getval(bfgs_gr))
                NM    NMgrad         CG     CGgrad      BFGS  BFGSgrad
shape1    3.576531  3.576531   3.576207   3.576211  3.576208  3.576206
shape2    1.187960  1.187960   1.187941   1.187942  1.187941  1.187941
loglik    4.587472  4.587472   4.587472   4.587472  4.587472  4.587472
function 47.000000 47.000000 226.000000 215.000000 14.000000 13.000000
gradient        NA        NA  83.000000  87.000000  5.000000  5.000000
> 
> 
> 
> llsurface(data = x, distr = "beta", plot.arg = c("shape1", "shape2"),
+           min.arg = c(0.1, 0.1), max.arg = c(7, 3), pal.col = heat.colors(50))
> points(bfgs$estimate[1], bfgs$estimate[2], pch="+", col="red")
> points(3, 3/4, pch="x", col="green")
> 
> 
> 
> proc.time()
   user  system elapsed 
   1.98    0.21    2.17