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 > > > data("groundbeef") > serving <- groundbeef$serving > > # (1) basic fit of a normal distribution with maximum likelihood estimation > # > > set.seed(1234) > x1 <- rnorm(n=100) > mledist(x1,"norm") $estimate mean sd -0.1567617 0.9993707 $convergence [1] 0 $value [1] 141.8309 $hessian mean sd mean 100.126 0.0000 sd 0.000 200.2538 $optim.function [1] "optim" $optim.method [1] "Nelder-Mead" $fix.arg NULL $fix.arg.fun NULL $weights NULL $counts function gradient 43 NA $optim.message NULL $loglik [1] -141.8309 $vcov NULL > > # (2) defining your own distribution functions, here for the Gumbel distribution > # for other distributions, see the CRAN task view dedicated to probability distributions > > dgumbel <- function(x,a,b) 1/b*exp((a-x)/b)*exp(-exp((a-x)/b)) > mledist(x1,"gumbel",start=list(a=10,b=5), silent=TRUE) $estimate a b -0.6267919 0.8564855 $convergence [1] 0 $value [1] 139.3812 $hessian a b a 136.31209 -61.53959 b -61.53959 270.09656 $optim.function [1] "optim" $optim.method [1] "Nelder-Mead" $fix.arg NULL $fix.arg.fun NULL $weights NULL $counts function gradient 91 NA $optim.message NULL $loglik [1] -139.3812 $vcov NULL > mledist(x1,"gumbel",start=list(a=10,b=5), silent=FALSE) $estimate a b -0.6267919 0.8564855 $convergence [1] 0 $value [1] 139.3812 $hessian a b a 136.31209 -61.53959 b -61.53959 270.09656 $optim.function [1] "optim" $optim.method [1] "Nelder-Mead" $fix.arg NULL $fix.arg.fun NULL $weights NULL $counts function gradient 91 NA $optim.message NULL $loglik [1] -139.3812 $vcov NULL Warning messages: 1: In log(do.call(ddistnam, c(list(obs), as.list(par), as.list(fix.arg)))) : NaNs produced 2: In log(do.call(ddistnam, c(list(obs), as.list(par), as.list(fix.arg)))) : NaNs produced > > # (3) fit a discrete distribution (Poisson) > # > > set.seed(1234) > x2 <- rpois(n=30,lambda = 2) > mledist(x2,"pois") $estimate lambda 1.7 $convergence [1] 0 $value [1] 46.18434 $hessian lambda lambda 17.64707 $optim.function [1] "optim" $optim.method [1] "BFGS" $fix.arg NULL $fix.arg.fun NULL $weights NULL $counts function gradient 6 1 $optim.message NULL $loglik [1] -46.18434 $vcov NULL > > # (4) fit a finite-support distribution (beta) > # > > set.seed(1234) > x3 <- rbeta(n=100,shape1=5, shape2=10) > mledist(x3,"beta") $estimate shape1 shape2 4.859798 10.918841 $convergence [1] 0 $value [1] -78.33052 $hessian shape1 shape2 shape1 16.295311 -6.542752 shape2 -6.542752 3.047900 $optim.function [1] "optim" $optim.method [1] "Nelder-Mead" $fix.arg NULL $fix.arg.fun NULL $weights NULL $counts function gradient 47 NA $optim.message NULL $loglik [1] 78.33052 $vcov NULL > > > # (5) fit frequency distributions on USArrests dataset. > # > > x4 <- USArrests$Assault > mledist(x4, "pois", silent=TRUE) $estimate lambda 170.76 $convergence [1] 0 $value [1] 1211.705 $hessian lambda lambda 0.2928087 $optim.function [1] "optim" $optim.method [1] "BFGS" $fix.arg NULL $fix.arg.fun NULL $weights NULL $counts function gradient 3 1 $optim.message NULL $loglik [1] -1211.705 $vcov NULL > mledist(x4, "pois", silent=FALSE) $estimate lambda 170.76 $convergence [1] 0 $value [1] 1211.705 $hessian lambda lambda 0.2928087 $optim.function [1] "optim" $optim.method [1] "BFGS" $fix.arg NULL $fix.arg.fun NULL $weights NULL $counts function gradient 3 1 $optim.message NULL $loglik [1] -1211.705 $vcov NULL > mledist(x4, "nbinom") $estimate size mu 3.822579 170.747853 $convergence [1] 0 $value [1] 290.3297 $hessian size mu size 1.759308e+00 -1.993783e-05 mu -1.993783e-05 6.412989e-03 $optim.function [1] "optim" $optim.method [1] "Nelder-Mead" $fix.arg NULL $fix.arg.fun NULL $weights NULL $counts function gradient 47 NA $optim.message NULL $loglik [1] -290.3297 $vcov NULL > > > # (6) scaling problem > # the simulated dataset (below) has particularly small values, hence without scaling (10^0), > # the optimization raises an error. The for loop shows how scaling by 10^i > # for i=1,...,6 makes the fitting procedure work correctly. > > set.seed(1234) > x2 <- rnorm(100, 1e-4, 2e-4) > for(i in 6:0) + cat(i, try(mledist(x*10^i, "cauchy")$estimate, silent=TRUE), "\n") 6 Error : object 'x' not found 5 Error : object 'x' not found 4 Error : object 'x' not found 3 Error : object 'x' not found 2 Error : object 'x' not found 1 Error : object 'x' not found 0 Error : object 'x' not found > > > # (7) scaling problem > # > > x <- c(-0.00707717, -0.000947418, -0.00189753, + -0.000474947, -0.00190205, -0.000476077, 0.00237812, 0.000949668, + 0.000474496, 0.00284226, -0.000473149, -0.000473373, 0, 0, 0.00283688, + -0.0037843, -0.0047506, -0.00238379, -0.00286807, 0.000478583, + 0.000478354, -0.00143575, 0.00143575, 0.00238835, 0.0042847, + 0.00237248, -0.00142281, -0.00142484, 0, 0.00142484, 0.000948767, + 0.00378609, -0.000472478, 0.000472478, -0.0014181, 0, -0.000946522, + -0.00284495, 0, 0.00331832, 0.00283554, 0.00141476, -0.00141476, + -0.00188947, 0.00141743, -0.00236351, 0.00236351, 0.00235794, + 0.00235239, -0.000940292, -0.0014121, -0.00283019, 0.000472255, + 0.000472032, 0.000471809, -0.0014161, 0.0014161, -0.000943842, + 0.000472032, -0.000944287, -0.00094518, -0.00189304, -0.000473821, + -0.000474046, 0.00331361, -0.000472701, -0.000946074, 0.00141878, + -0.000945627, -0.00189394, -0.00189753, -0.0057143, -0.00143369, + -0.00383326, 0.00143919, 0.000479272, -0.00191847, -0.000480192, + 0.000960154, 0.000479731, 0, 0.000479501, 0.000958313, -0.00383878, + -0.00240674, 0.000963391, 0.000962464, -0.00192586, 0.000481812, + -0.00241138, -0.00144963) > > > #only i == 0, no scaling, should not converge. > for(i in 6:0) + cat(i, try(mledist(x*10^i, "cauchy")$estimate, silent=TRUE), "\n") 6 -290.2967 1193.769 5 -29.03782 119.3099 4 -2.897257 11.93096 3 -0.2897257 1.193096 2 -0.02898858 0.1193277 1 -0.002897257 0.01193096 0 NA NA > > > # (8) normal mixture > # > > #mixture of two normal distributions > #density > dnorm2 <- function(x, w, m1, s1, m2, s2) + w*dnorm(x, m1, s1) + (1-w)*dnorm(x, m2, s2) > #numerically-approximated quantile function > qnorm2 <- function(p, w, m1, s1, m2, s2) + { + L2 <- function(x, prob) + (prob - pnorm2(x, w, m1, s1, m2, s2))^2 + sapply(p, function(pr) optimize(L2, c(-20, 30), prob=pr)$minimum) + } > #distribution function > pnorm2 <- function(q, w, m1, s1, m2, s2) + w*pnorm(q, m1, s1) + (1-w)*pnorm(q, m2, s2) > > > #basic normal distribution > x <- c(rnorm(1000, 5), rnorm(1000, 10)) > #MLE fit > fit1 <- mledist(x, "norm2", start=list(w=1/3, m1=4, s1=2, m2=8, s2=2), + lower=c(0, 0, 0, 0, 0)) > > > > > > > # (9) fit a Pareto and log-logistic distributions > # > > if(any(installed.packages()[,"Package"] == "actuar")) + { + require(actuar) + #simulate a sample + x4 <- rpareto(1000, 6, 2) + + #fit + mledist(x4, "pareto", start=list(shape=10, scale=10), lower=1, upper=Inf) + + + #simulate a sample + x4 <- rllogis(1000, 6, 2) + + #fit + mledist(x4, "llogis", start=list(shape=10, rate=1), lower=1, upper=Inf) + + } Loading required package: actuar Attaching package: 'actuar' The following object is masked _by_ '.GlobalEnv': dgumbel The following objects are masked from 'package:stats': sd, var The following object is masked from 'package:grDevices': cm $estimate [1] NA NA $convergence [1] 100 $loglik [1] NA $hessian [1] NA $optim.function [1] "constrOptim" $fix.arg NULL $optim.method [1] "Nelder-Mead" $fix.arg.fun NULL $counts [1] NA NA $vcov NULL > > > > # (10) custom optim for exponential distribution > # > if(any(installed.packages()[,"Package"] == "rgenoud") && FALSE) + { + + + mysample <- rexp(1000, 5) + mystart <- list(rate=8) + + fNM <- mledist(mysample, "exp", optim.method="Nelder-Mead") + fBFGS <- mledist(mysample, "exp", optim.method="BFGS") + fLBFGSB <- mledist(mysample, "exp", optim.method="L-BFGS-B", lower=0) + fSANN <- mledist(mysample, "exp", optim.method="SANN") + fCG <- try(mledist(mysample, "exp", optim.method="CG") ) + if(inherits(fCG, "try-error")) + fCG <- list(estimate=NA) + + #the warning tell us to use optimise... + + #to meet the standard 'fn' argument and specific name arguments, we wrap optimize, + myoptimize <- function(fn, par, ...) + { + res <- optimize(f=fn, ..., maximum=FALSE) + c(res, convergence=0, value=res$objective, par=res$minimum, hessian=NA) + } + + foptimize <- mledist(mysample, "exp", start=mystart, custom.optim=myoptimize, interval=c(0, 100)) + + + library(rgenoud) + + #wrap genoud function rgenoud package + mygenoud <- function(fn, par, ...) + { + res <- genoud(fn, starting.values=par, ...) + c(res, convergence=0, counts=NULL) + } + + + fgenoud <- mledist(mysample, "exp", start=mystart, custom.optim= mygenoud, nvars=1, + Domains=cbind(0, 10), boundary.enforcement=1, + hessian=TRUE, print.level=0) + + + c(NM=fNM$estimate, + BFGS=fBFGS$estimate, + LBFGSB=fLBFGSB$estimate, + SANN=fSANN$estimate, + CG=fCG$estimate, + optimize=foptimize$estimate, + fgenoud=fgenoud$estimate) + + } > > > # (11) custom optim for gamma distribution > # > if(any(installed.packages()[,"Package"] == "rgenoud") && FALSE) + { + + + mysample <- rgamma(1000, 5, 3) + mystart <- c(shape=10, rate=10) + + fNM <- mledist(mysample, "gamma", optim.method="Nelder-Mead") + fBFGS <- mledist(mysample, "gamma", optim.method="BFGS") + fLBFGSB <- mledist(mysample, "gamma", optim.method="L-BFGS-B", lower=0) + fSANN <- mledist(mysample, "gamma", optim.method="SANN") + fCG <- try( mledist(mysample, "gamma", optim.method="CG", control=list(maxit=1000)) ) + if(inherits(fCG, "try-error")) + fCG <- list(estimate=NA) + + fgenoud <- mledist(mysample, "gamma", start=mystart, custom.optim= mygenoud, nvars=2, + Domains=cbind(c(0,0), c(100,100)), boundary.enforcement=1, + hessian=TRUE, print.level=0) + + cbind(NM=fNM$estimate, + BFGS=fBFGS$estimate, + LBFGSB=fLBFGSB$estimate, + SANN=fSANN$estimate, + CG=fCG$estimate, + fgenoud=fgenoud$estimate) + + + data(groundbeef) + + fNM <- mledist(groundbeef$serving, "gamma", optim.method="Nelder-Mead") + fBFGS <- mledist(groundbeef$serving, "gamma", optim.method="BFGS") + fLBFGSB <- mledist(groundbeef$serving, "gamma", optim.method="L-BFGS-B", lower=0) + fSANN <- mledist(groundbeef$serving, "gamma", optim.method="SANN") + fCG <- try( mledist(groundbeef$serving, "gamma", optim.method="CG", control=list(maxit=10000)) ) + if(inherits(fCG, "try-error")) + fCG <- list(estimate=NA) + + fgenoud <- mledist(groundbeef$serving, "gamma", start=list(shape=4, rate=1), + custom.optim= mygenoud, nvars=2, max.generations=10, + Domains=cbind(c(0,0), c(10,10)), boundary.enforcement=1, + hessian=TRUE, print.level=0, P9=10) + + cbind(NM=fNM$estimate, + BFGS=fBFGS$estimate, + LBFGSB=fLBFGSB$estimate, + SANN=fSANN$estimate, + CG=fCG$estimate, + fgenoud=fgenoud$estimate) + + + } > > > > # (12) test error messages > # > > dnorm2 <- function(x, a) + "NA" > x <- rexp(10) > > #should get a one-line error > res <- mledist(x, "norm2", start=list(a=1)) > #as in > attr(try(log("a"), silent=TRUE), "condition") > > > try(mledist(c(serving, "a"), "gamma")) Error in mledist(c(serving, "a"), "gamma") : data must be a numeric vector of length greater than 1. > try(mledist(c(serving, NA), "gamma")) Error in checkUncensoredNAInfNan(data) : data contain NA values. > try(mledist(c(serving, Inf), "gamma")) Error in checkUncensoredNAInfNan(data) : data contain Inf (infinite) values. > try(mledist(c(serving, -Inf), "gamma")) Error in checkUncensoredNAInfNan(data) : data contain Inf (infinite) values. > try(mledist(c(serving, NaN), "gamma")) Error in checkUncensoredNAInfNan(data) : data contain NaN (not a numeric) values. > > > > # (13) weighted MLE > # > n <- 1e6 > n <- 1e2 > x <- rpois(n, 10) > xtab <- table(x) > xval <- sort(unique(x)) > f1 <- mledist(x, "pois", start=list(lambda=mean(x)), optim.method="Brent", lower=0, + upper=100, control=list(trace=1)) > f2 <- mledist(xval, "pois", weights=xtab, start=list(lambda=mean(x))) Warning message: In mledist(xval, "pois", weights = xtab, start = list(lambda = mean(x))) : weights are not taken into account in the default initial values > > f1$estimate lambda 9.87 > f2$estimate #should be identical lambda 9.87 > > #test discrete distrib > f2 <- try(mledist(xval, "pois", weights=1:length(xval), start=list(lambda=mean(x)))) Warning message: In mledist(xval, "pois", weights = 1:length(xval), start = list(lambda = mean(x))) : weights are not taken into account in the default initial values > #test non integer weights > f2 <- try(mledist(xval, "pois", weights=rep(1/3, length(xval)), start=list(lambda=mean(x)))) Error in mledist(xval, "pois", weights = rep(1/3, length(xval)), start = list(lambda = mean(x))) : weights should be a vector of (strictly) positive integers > f2 <- try(mledist(1:10, "pois", weights=c(rep(1, 9), 1.001), start=list(lambda=mean(x)))) Error in mledist(1:10, "pois", weights = c(rep(1, 9), 1.001), start = list(lambda = mean(x))) : weights should be a vector of (strictly) positive integers > f2 <- try(mledist(1:10, "pois", weights=c(rep(1, 9), 1.0000001), start=list(lambda=mean(x)))) Error in mledist(1:10, "pois", weights = c(rep(1, 9), 1.0000001), start = list(lambda = mean(x))) : weights should be a vector of (strictly) positive integers > > # (14) no convergence > # > n <- 1e2 > x <- c(rep(0, n), rpois(n, 10), rpois(n, 50)) > > mledist(x, "pois", optim.method="Nelder-Mead", control=list(maxit=10)) $estimate lambda 20.19 $convergence [1] 1 $value [1] 4473.33 $hessian lambda lambda 14.85884 $optim.function [1] "optim" $optim.method [1] "Nelder-Mead" $fix.arg NULL $fix.arg.fun NULL $weights NULL $counts function gradient 12 NA $optim.message NULL $loglik [1] -4473.33 $vcov NULL > > > > > # (15) basic fit of a normal distribution with new fix.arg/start.arg > # > > set.seed(1234) > x1 <- rnorm(n=100) > > #correct usage > mledist(x1,"norm") $estimate mean sd -0.1567617 0.9993707 $convergence [1] 0 $value [1] 141.8309 $hessian mean sd mean 100.126 0.0000 sd 0.000 200.2538 $optim.function [1] "optim" $optim.method [1] "Nelder-Mead" $fix.arg NULL $fix.arg.fun NULL $weights NULL $counts function gradient 43 NA $optim.message NULL $loglik [1] -141.8309 $vcov NULL > mledist(x1,"norm", start=function(x) list(mean=0, sd=1)) $estimate mean sd -0.1568775 0.9993092 $convergence [1] 0 $value [1] 141.8309 $hessian mean sd mean 100.13830523 0.02319864 sd 0.02319864 200.31537814 $optim.function [1] "optim" $optim.method [1] "Nelder-Mead" $fix.arg NULL $fix.arg.fun NULL $weights NULL $counts function gradient 47 NA $optim.message NULL $loglik [1] -141.8309 $vcov NULL > mledist(x1,"norm", fix.arg=function(x) list(mean=mean(x))) $estimate sd 0.9993707 $convergence [1] 0 $value [1] 141.8309 $hessian sd sd 200.2538 $optim.function [1] "optim" $optim.method [1] "BFGS" $fix.arg $fix.arg$mean [1] -0.1567617 $fix.arg.fun function (x) list(mean = mean(x)) $weights NULL $counts function gradient 8 1 $optim.message NULL $loglik [1] -141.8309 $vcov NULL > mledist(x1,"norm", fix.arg=list(mean=1/2)) $estimate sd 1.195858 $convergence [1] 0 $value [1] 159.7803 $hessian sd sd 139.8536 $optim.function [1] "optim" $optim.method [1] "BFGS" $fix.arg $fix.arg$mean [1] 0.5 $fix.arg.fun NULL $weights NULL $counts function gradient 27 6 $optim.message NULL $loglik [1] -159.7803 $vcov NULL > mledist(x1,"norm", fix.arg=list(mean=1/2), start=list(sd=1)) $estimate sd 1.195858 $convergence [1] 0 $value [1] 159.7803 $hessian sd sd 139.8536 $optim.function [1] "optim" $optim.method [1] "BFGS" $fix.arg $fix.arg$mean [1] 0.5 $fix.arg.fun NULL $weights NULL $counts function gradient 25 6 $optim.message NULL $loglik [1] -159.7803 $vcov NULL > mledist(x1,"norm", fix.arg=function(x) list(mean=0), start=list(sd=1)) $estimate sd 1.011589 $convergence [1] 0 $value [1] 143.0463 $hessian sd sd 195.4463 $optim.function [1] "optim" $optim.method [1] "BFGS" $fix.arg $fix.arg$mean [1] 0 $fix.arg.fun function (x) list(mean = 0) $weights NULL $counts function gradient 12 4 $optim.message NULL $loglik [1] -143.0463 $vcov NULL > mledist(x1,"norm", fix.arg=function(x) list(mean=mean(x)), start=list(sd=1)) $estimate sd 0.9993705 $convergence [1] 0 $value [1] 141.8309 $hessian sd sd 200.2539 $optim.function [1] "optim" $optim.method [1] "BFGS" $fix.arg $fix.arg$mean [1] -0.1567617 $fix.arg.fun function (x) list(mean = mean(x)) $weights NULL $counts function gradient 11 3 $optim.message NULL $loglik [1] -141.8309 $vcov NULL > > #wrong usage (see text message in util-checkparam.R) > try( mledist(x1,"norm", start=list(a=1/2)) ) #t3 Error in checkparamlist(arg_startfix$start.arg, arg_startfix$fix.arg, : 'start' must specify names which are arguments to 'distr'. > try( mledist(x1,"norm", start=function(x) list(a=0, b=1)) ) #t3 Error in checkparamlist(arg_startfix$start.arg, arg_startfix$fix.arg, : 'start' must specify names which are arguments to 'distr'. > try( mledist(x1,"norm", fix.arg=list(a=1/2)) ) #t4 Error in checkparamlist(arg_startfix$start.arg, arg_startfix$fix.arg, : 'fix.arg' must specify names which are arguments to 'distr'. > try( mledist(x1,"norm", fix.arg=function(x) list(a=0), start=list(sd=1)) ) #t4 Error in checkparamlist(arg_startfix$start.arg, arg_startfix$fix.arg, : 'fix.arg' must specify names which are arguments to 'distr'. > try( mledist(x1,"norm", start=matrix(1/2)) ) #t1 Error in manageparam(start.arg = start, fix.arg = fix.arg, obs = data, : Wrong type of argument for start > try( mledist(x1,"norm", fix.arg=matrix(1/2)) ) #t0 Error in manageparam(start.arg = start, fix.arg = fix.arg, obs = data, : Wrong type of argument for fix.arg > try( mledist(x1,"norm", fix.arg=matrix(1/2), start=matrix(1/2)) ) #t2 Error in manageparam(start.arg = start, fix.arg = fix.arg, obs = data, : Wrong type of argument for start > try( mledist(x1,"norm", fix.arg=function(x) list(mean=mean(x), sd=2), start=list(sd=1)) ) #t5 Error in checkparamlist(arg_startfix$start.arg, arg_startfix$fix.arg, : A distribution parameter cannot be specified both in 'start' and 'fix.arg'. > dabcnorm <- function(x, mean, sd) 1 > try( mledist(x1,"abcnorm", fix.arg=function(x) list(mean=mean(x))) ) #t8 Error in computing default starting values. Error in manageparam(start.arg = start, fix.arg = fix.arg, obs = data, : Error in startargdefault(obs, distname) : Unknown starting values for distribution abcnorm. > > > > # (16) relevant example for zero modified geometric distribution > # > dzmgeom <- function(x, p1, p2) + { + p1 * (x == 0) + (1-p1)*dgeom(x-1, p2) + } > rzmgeom <- function(n, p1, p2) + { + u <- rbinom(n, 1, 1-p1) #prob to get zero is p1 + u[u != 0] <- rgeom(sum(u != 0), p2)+1 + u + } > # check > # dzmgeom(0:5, 1/2, 1/10) > > x2 <- rzmgeom(100, 1/2, 1/10) > x3 <- rzmgeom(100, 1/3, 1/10) > x4 <- rzmgeom(100, 1/4, 1/10) > table(x2) x2 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 18 19 23 29 48 7 6 6 5 3 1 2 2 3 1 1 1 2 4 1 1 1 1 2 2 > #this is the MLE which converges almost surely and in distribution to the true value. > initp1 <- function(x) list(p1=mean(x == 0)) > > mledist(x2, "zmgeom", fix.arg=initp1, start=list(p2=1/2))[c("estimate", "fix.arg")] $estimate p2 0.1232217 $fix.arg $fix.arg$p1 [1] 0.48 > mledist(x3, "zmgeom", fix.arg=initp1, start=list(p2=1/2))[c("estimate", "fix.arg")] $estimate p2 0.09494948 $fix.arg $fix.arg$p1 [1] 0.38 > mledist(x4, "zmgeom", fix.arg=initp1, start=list(p2=1/2))[c("estimate", "fix.arg")] $estimate p2 0.09494494 $fix.arg $fix.arg$p1 [1] 0.23 > > > # (17) test the component optim.message > x <- rnorm(100) > #change parameter to obtain unsuccessful convergence > mledist(x, "norm", control=list(maxit=2), start=list(mean=1e5, sd=1), optim.method="L-BFGS-B", lower=0) $estimate mean sd 1.000000e+05 2.764093e+00 $convergence [1] 1 $value [1] 65443342926 $hessian mean sd mean 13.35144 -947049.1 sd -947049.14093 51393930814.7 $optim.function [1] "optim" $optim.method [1] "L-BFGS-B" $fix.arg NULL $fix.arg.fun NULL $weights NULL $counts function gradient 4 4 $optim.message [1] "NEW_X" $loglik [1] -65443342926 $vcov NULL > > > # (18) management of bounds in optim/constrOptim > x <- rexp(100) > mledist(x, "exp") #optim, BFGS $estimate rate 1.001466 $convergence [1] 0 $value [1] 99.85353 $hessian rate rate 99.70769 $optim.function [1] "optim" $optim.method [1] "BFGS" $fix.arg NULL $fix.arg.fun NULL $weights NULL $counts function gradient 8 1 $optim.message NULL $loglik [1] -99.85353 $vcov NULL > mledist(x, "exp", optim.method="Brent", lower=0, upper=100) #optim, Brent $estimate rate 1.001466 $convergence [1] 0 $value [1] 99.85353 $hessian [,1] [1,] 99.70769 $optim.function [1] "optim" $optim.method [1] "Brent" $fix.arg NULL $fix.arg.fun NULL $weights NULL $counts function gradient NA NA $optim.message NULL $loglik [1] -99.85353 $vcov NULL > mledist(x, "exp", optim.method="Nelder-Mead") #optim, Nelder-Mead $estimate rate 1.001466 $convergence [1] 0 $value [1] 99.85353 $hessian rate rate 99.70769 $optim.function [1] "optim" $optim.method [1] "Nelder-Mead" $fix.arg NULL $fix.arg.fun NULL $weights NULL $counts function gradient 22 NA $optim.message NULL $loglik [1] -99.85353 $vcov NULL > mledist(x, "exp", lower=0, optim.method="Nelder-Mead") #constrOptim, Nelder-Mead $estimate rate 1.001466 $convergence [1] 0 $value [1] 99.85353 $hessian rate rate 99.70769 $optim.function [1] "constrOptim" $optim.method [1] "Nelder-Mead" $fix.arg NULL $fix.arg.fun NULL $weights NULL $counts [1] 0 NA $optim.message NULL $loglik [1] -99.85353 $vcov NULL > mledist(x, "exp", lower=0, optim.method="BFGS") #optim, L-BFGS-B $estimate rate 1.001466 $convergence [1] 0 $value [1] 99.85353 $hessian rate rate 99.70769 $optim.function [1] "optim" $optim.method [1] "L-BFGS-B" $fix.arg NULL $fix.arg.fun NULL $weights NULL $counts function gradient 10 10 $optim.message [1] "CONVERGENCE: REL_REDUCTION_OF_F <= FACTR*EPSMCH" $loglik [1] -99.85353 $vcov NULL Warning message: In mledist(x, "exp", lower = 0, optim.method = "BFGS") : The BFGS method cannot be used with bounds without provided the gradient. The method is changed to L-BFGS-B. > > > x <- rbeta(100, 3/2, 7/3) > mledist(x, "beta", optim.method="Nelder") #optim, Nelder-Mead $estimate shape1 shape2 1.666939 2.568482 $convergence [1] 0 $value [1] -18.69748 $hessian shape1 shape2 shape1 54.75563 -26.61467 shape2 -26.61467 20.85442 $optim.function [1] "optim" $optim.method [1] "Nelder-Mead" $fix.arg NULL $fix.arg.fun NULL $weights NULL $counts function gradient 49 NA $optim.message NULL $loglik [1] 18.69748 $vcov NULL > mledist(x, "beta", lower=0, optim.method="Nelder-Mead") #constrOptim, Nelder-Mead $estimate shape1 shape2 1.666939 2.568482 $convergence [1] 0 $value [1] -18.69748 $hessian shape1 shape2 shape1 55.36355 -26.63521 shape2 -26.63521 20.68112 $optim.function [1] "constrOptim" $optim.method [1] "Nelder-Mead" $fix.arg NULL $fix.arg.fun NULL $weights NULL $counts function gradient 49 NA $optim.message NULL $loglik [1] 18.69748 $vcov NULL > #as the result of optim(c(-1.2,1), fr, method = "Nelder-Mead", hessian=TRUE, gr=NULL, lower=-Inf, upper=Inf) from optim() example > > mledist(x, "beta", lower=0, optim.method="BFGS") #optim, L-BFGS-B $estimate shape1 shape2 1.666899 2.568436 $convergence [1] 0 $value [1] -18.69748 $hessian shape1 shape2 shape1 54.75756 -26.61528 shape2 -26.61528 20.85484 $optim.function [1] "optim" $optim.method [1] "L-BFGS-B" $fix.arg NULL $fix.arg.fun NULL $weights NULL $counts function gradient 8 8 $optim.message [1] "CONVERGENCE: REL_REDUCTION_OF_F <= FACTR*EPSMCH" $loglik [1] 18.69748 $vcov NULL Warning message: In mledist(x, "beta", lower = 0, optim.method = "BFGS") : The BFGS method cannot be used with bounds without provided the gradient. The method is changed to L-BFGS-B. > > proc.time() user system elapsed 2.09 0.32 2.42