R Under development (unstable) (2026-01-21 r89314 ucrt) -- "Unsuffered Consequences"
Copyright (C) 2026 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 <- 1000
> nsample <- 10
> 
> set.seed(123) # here just to make random sampling reproducible
> 
> 
> #(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.6075044  0.9437128  0.7551198  0.7761835  0.9459930  0.8611853
 [7]  0.2021912  0.8493394  0.7697487  0.3459894 -0.6211948 -4.4557160
[13] -1.1145620 -1.2187381 -4.9050493 -1.8486153 -0.0246526 -1.7292425
[19] -1.1855091 -0.1743470
> 
> 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]   7.056968 -17.277627
> grlnL(c(3, 4), x)
[1] -7.056968 17.277627
> grlnL(c(3, 3/4), x)
[1] -0.1544983  0.7603107
> 
> 
> 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.119724  3.119724   3.1193176   3.1192761  3.1194693  3.119574
shape2    0.732445  0.732445   0.7324581   0.7324524  0.7324777  0.732491
loglik    6.894663  6.894663   6.8946631   6.8946631  6.8946631  6.894663
function 49.000000 49.000000 322.0000000 824.0000000 30.0000000 11.000000
gradient        NA        NA 161.0000000 277.0000000 17.0000000  6.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.29    0.07    1.37