R Under development (unstable) (2025-07-02 r88374 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 <- 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.8783655  0.6145503  0.6383075  0.9141985  0.7882158  0.8266081
 [7]  0.9210092  0.9084386  0.9046738  0.7616455 -2.0551152 -0.6386573
[13] -0.7008729 -2.7309452 -1.2846269 -1.5367851 -2.9385715 -2.5846253
[19] -2.4997799 -1.1454623
> 
> 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]   8.156013 -18.115442
> grlnL(c(3, 4), x)
[1] -8.156013 18.115442
> grlnL(c(3, 3/4), x)
[1] -1.253543  1.598125
> 
> 
> 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   11.293895 11.293895   11.288921   11.29306  11.266650 11.293047
shape2    1.553652  1.553652    1.553195    1.55361   1.551079  1.553615
loglik   11.930286 11.930286   11.930286   11.93029  11.930274 11.930286
function 47.000000 47.000000 2001.000000 1719.00000 104.000000 22.000000
gradient        NA        NA 1001.000000  689.00000 102.000000 15.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 
   2.15    0.40    2.50