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