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.5594888  0.6225404  0.8363072  0.8920901  0.4496241  0.9441497
 [7]  0.6278384  0.8286124  0.4549052  0.9218214 -0.5125414 -0.6589899
[13] -1.6138967 -2.2610162 -0.3155062 -4.5275256 -0.6727954 -1.5521796
[19] -0.3236943 -2.9665812
> 
> 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.137378 -15.404726
> grlnL(c(3, 4), x)
[1] -7.137378 15.404726
> grlnL(c(3, 3/4), x)
[1] -0.2349083 -1.1125896
> 
> 
> 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    4.0986062  4.0986062   4.0988926   4.0988894  4.0988833  4.0988808
shape2    0.9645588  0.9645588   0.9645577   0.9645571  0.9645544  0.9645537
loglik    6.7949900  6.7949900   6.7949900   6.7949900  6.7949900  6.7949900
function 57.0000000 57.0000000 386.0000000 360.0000000 18.0000000 18.0000000
gradient         NA         NA 141.0000000 143.0000000  9.0000000  9.0000000
> 
> 
> 
> 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.51    0.35    1.86