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 <- 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.8980844  0.6929612  0.9156693  0.9112897  0.7032783  0.8625993
 [7]  0.9262892  0.8655583  0.4999978  0.9358063 -2.3673151 -0.8676297
[13] -2.7721655 -2.6542755 -0.9036533 -1.8639669 -3.1369915 -1.8969510
[19] -0.3984138 -3.6453771
> 
> 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.211534 -20.506739
> grlnL(c(3, 4), x)
[1] -8.211534 20.506739
> grlnL(c(3, 3/4), x)
[1] -1.309064  3.989423
> 
> 
> 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    8.021163  8.021163    8.020062    8.020310  8.003343  8.020332
shape2    1.035606  1.035606    1.035568    1.035587  1.034187  1.035591
loglik   11.783666 11.783666   11.783666   11.783666 11.783656 11.783666
function 53.000000 53.000000 1819.000000 1287.000000 50.000000 19.000000
gradient        NA        NA  911.000000  515.000000 48.000000 10.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.34    0.17    1.50