R Under development (unstable) (2025-06-30 r88369 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.87613420 -0.35048197 0.94684135 0.52082585 0.92596452 0.07772637 [7] 0.94935142 -0.76656641 0.91099077 0.35341876 -2.02554426 0.29866219 [13] -5.14252122 -0.43615005 -3.12355125 0.07538168 -6.72437889 0.41860340 [19] -2.64673146 -0.18337799 > > 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] 4.444205 -19.489608 > grlnL(c(3, 4), x) [1] -4.444205 19.489608 > grlnL(c(3, 3/4), x) [1] 2.458265 2.972292 > > > 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 0.9993314 0.9993314 0.9993461 0.9993472 0.9993548 0.9993470 shape2 0.3897140 0.3897140 0.3896660 0.3896655 0.3896650 0.3896652 loglik 6.2342447 6.2342447 6.2342448 6.2342448 6.2342448 6.2342448 function 41.0000000 41.0000000 113.0000000 106.0000000 9.0000000 20.0000000 gradient NA NA 47.0000000 33.0000000 6.0000000 5.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.70 0.23 1.92