R Under development (unstable) (2026-01-22 r89323 ucrt) -- "Unsuffered Consequences" Copyright (C) 2026 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 > set.seed(123) # here just to make random sampling reproducible > > data("groundbeef") > serving <- groundbeef$serving > > # (1) basic fit of a normal distribution with maximum likelihood estimation > # > > x1 <- rnorm(n=100) > mledist(x1,"norm") $estimate mean sd 0.09040591 0.90824033 $convergence [1] 0 $value [1] 1.322692 $hessian mean sd mean 1.212267e+00 5.551115e-11 sd 5.551115e-11 2.424561e+00 $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] -132.2692 $vcov NULL > > #fitted coef around -0.1567617 0.9993707, fitted loglik -141.8309 > > # (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) + { + cat(a, b, "\n") + res <- 1/b*exp((a-x)/b)*exp(-exp((a-x)/b)) + if(any(is.infinite(log(res)))) + stop(log(res)) + res + } > mledist(x1,"gumbel",start=list(a=10,b=5), silent=TRUE) 10 5 11 5 10 6 9 6 8 6.5 8 7.5 7 8.75 5 9.25 2.5 10.875 4 11.5 5 10.25 3 10.75 4 10.25 4 9.25 3.5 8.75 4.5 7.75 4.75 6.5 3.25 6 2.375 4.375 3.625 2.125 3.53125 7.09375 1.15625 4.96875 -0.640625 4.203125 -1.796875 1.484375 -4.460938 -1.320312 -4.8125 1.3125 0.578125 3.609375 -0.578125 0.890625 -0.546875 -0.765625 -2.953125 -1.234375 -0.3046875 2.398438 0.9140625 1.804688 0.2363281 1.724609 -0.03710938 0.2167969 $estimate [1] NA NA $convergence [1] 100 $loglik [1] NA $hessian [1] NA $optim.function [1] "optim" $fix.arg NULL $optim.method [1] "Nelder-Mead" $fix.arg.fun NULL $counts [1] NA NA $vcov NULL > mledist(x1,"gumbel",start=list(a=10,b=5), silent=FALSE, control=list(trace=1), lower=0) 10 5 10 5 Nelder-Mead direct search function minimizer 10 5 function value for initial parameters = 7.002787 Scaled convergence tolerance is 1.0435e-07 Stepsize computed as 1.000000 11 5 10 6 BUILD 3 8.436053 5.413769 9 6 8 6.5 EXTENSION 5 7.002787 4.062969 8 7.5 7 8.75 EXTENSION 7 5.413769 3.592467 5 9.25 2.5 10.875 REFLECTION 9 4.062969 3.400994 4 11.5 5 10.25 LO-REDUCTION 11 3.592467 3.400994 3 10.75 4 10.25 LO-REDUCTION 13 3.467818 3.400994 4 9.25 3.5 8.75 EXTENSION 15 3.414846 3.262722 4.5 7.75 4.75 6.5 EXTENSION 17 3.400994 3.221595 3.25 6 2.375 4.375 EXTENSION 19 3.262722 2.675052 3.625 2.125 3.53125 7.09375 HI-REDUCTION 21 3.221595 2.675052 1.15625 4.96875 REFLECTION 23 3.110574 2.648370 0 2.25 2.648438 5.882812 HI-REDUCTION 25 2.899356 2.648370 0.8828125 3.460938 0 2.25 REFLECTION 27 2.675052 2.313586 1.697266 4.294922 HI-REDUCTION 29 2.648370 2.313586 1.423828 2.787109 1.557617 1.696289 REFLECTION 31 2.569083 2.246965 0.609375 1.953125 0.06542969 0.7822266 EXTENSION 33 2.313586 1.651485 0.6064453 0.1083984 $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 > > #fitted coef around -0.6267919 0.8564855, fitted loglik -139.3812 > > > # (3) fit a discrete distribution (Poisson) > # > > x2 <- rpois(n=30,lambda = 2) > mledist(x2,"pois") $estimate lambda 2.066667 $convergence [1] 0 $value [1] 1.663729 $hessian lambda lambda 0.4838712 $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] -49.91188 $vcov NULL > > #fitted coef around 1.7, fitted loglik -46.18434 > > # (4) fit a finite-support distribution (beta) > # > > x3 <- rbeta(n=100,shape1=5, shape2=10) > mledist(x3,"beta") $estimate shape1 shape2 5.344291 11.259568 $convergence [1] 0 $value [1] -0.7917456 $hessian shape1 shape2 shape1 0.14362910 -0.06207699 shape2 -0.06207699 0.03079685 $optim.function [1] "optim" $optim.method [1] "Nelder-Mead" $fix.arg NULL $fix.arg.fun NULL $weights NULL $counts function gradient 45 NA $optim.message NULL $loglik [1] 79.17456 $vcov NULL > > #fitted coef around 4.859798 10.918841, fitted loglik 78.33052 > > # (5) fit frequency distributions on USArrests dataset. > # > > x4 <- USArrests$Assault > mledist(x4, "pois", silent=TRUE) $estimate lambda 170.76 $convergence [1] 0 $value [1] 24.2341 $hessian lambda lambda 0.005856174 $optim.function [1] "optim" $optim.method [1] "BFGS" $fix.arg NULL $fix.arg.fun NULL $weights NULL $counts function gradient 2 1 $optim.message NULL $loglik [1] -1211.705 $vcov NULL > mledist(x4, "pois", silent=FALSE) $estimate lambda 170.76 $convergence [1] 0 $value [1] 24.2341 $hessian lambda lambda 0.005856174 $optim.function [1] "optim" $optim.method [1] "BFGS" $fix.arg NULL $fix.arg.fun NULL $weights NULL $counts function gradient 2 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] 5.806593 $hessian size mu size 3.518616e-02 -3.985701e-07 mu -3.985701e-07 1.282603e-04 $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 > > #fitted coef around 3.822579 170.747853, fitted loglik -290.3297 > > > # (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. > > 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) > #distribution function > pnorm2 <- function(q, w, m1, s1, m2, s2) + w*pnorm(q, m1, s1) + (1-w)*pnorm(q, 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) + } > > > #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), upper=c(1, Inf, Inf, Inf, Inf)) > fit1 $estimate w m1 s1 m2 s2 0.4977092 5.0051861 0.9941993 9.9848285 0.9940721 $convergence [1] 0 $value [1] 2.086962 $hessian w m1 s1 m2 s2 w 1.97996629 -0.166692761 0.03086102 -0.32402098 0.283795512 m1 -0.16669276 0.047950966 0.02998200 -0.01320996 0.003391253 s1 0.03086102 0.029981996 -0.05059679 0.03149497 -0.031331318 m2 -0.32402098 -0.013209963 0.03149497 0.11942163 0.168809354 s2 0.28379551 0.003391253 -0.03133132 0.16880935 0.498087013 $optim.function [1] "constrOptim" $optim.method [1] "Nelder-Mead" $fix.arg NULL $fix.arg.fun NULL $weights NULL $counts function gradient 636 NA $optim.message NULL $loglik [1] -4173.925 $vcov NULL > > fit1 <- fitdist(x, "norm2", start=list(w=1/3, m1=4, s1=2, m2=8, s2=2), + lower=c(0, 0, 0, 0, 0), upper=c(1, Inf, Inf, Inf, Inf)) > denscomp(fit1) > > > #fitted coef around 0.5003298 4.9842719 0.9909527 10.0296973 1.0024444 , fitted loglik -4185.114 > > # (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) + + #fitted coef around 6.334596 2.110411, fitted loglik -58.73242 + + #simulate a sample + x4 <- rllogis(1000, 6, 2) + + #fit + mledist(x4, "llogis", start=list(shape=10, rate=2), lower=1, upper=Inf) + + #fitted coef around 6.150001 2.005405, fitted loglik 511.873 + } 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 shape rate 6.008014 2.005417 $convergence [1] 0 $value [1] -0.4871117 $hessian shape rate shape 0.014592765 0.002855195 rate 0.002855195 6.325814925 $optim.function [1] "constrOptim" $optim.method [1] "Nelder-Mead" $fix.arg NULL $fix.arg.fun NULL $weights NULL $counts function gradient 118 NA $optim.message NULL $loglik [1] 487.1117 $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)) + + + require("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 10.3 > f2$estimate #should be identical lambda 10.3 > > #fitted coef around 9.74, fitted loglik -246.8505 > > #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 19.74667 $convergence [1] 1 $value [1] 14.89852 $hessian lambda lambda 0.05064146 $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] -4469.557 $vcov NULL > > > # (15) basic fit of a normal distribution with new fix.arg/start.arg > # > > x1 <- rnorm(n=100) > > #correct usage > mledist(x1,"norm") $estimate mean sd 0.07560489 1.05471103 $convergence [1] 0 $value [1] 1.472205 $hessian mean sd mean 0.8989448 0.000000 sd 0.0000000 1.797904 $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] -147.2205 $vcov NULL > mledist(x1,"norm", start=function(x) list(mean=0, sd=1)) $estimate mean sd 0.07542446 1.05481014 $convergence [1] 0 $value [1] 1.472205 $hessian mean sd mean 0.8987758787 0.0003074916 sd 0.0003074916 1.7970597157 $optim.function [1] "optim" $optim.method [1] "Nelder-Mead" $fix.arg NULL $fix.arg.fun NULL $weights NULL $counts function gradient 39 NA $optim.message NULL $loglik [1] -147.2205 $vcov NULL > mledist(x1,"norm", fix.arg=function(x) list(mean=mean(x))) $estimate sd 1.054711 $convergence [1] 0 $value [1] 1.472205 $hessian sd sd 1.797904 $optim.function [1] "optim" $optim.method [1] "BFGS" $fix.arg $fix.arg$mean [1] 0.07560489 $fix.arg.fun function (x) list(mean = mean(x)) $weights NULL $counts function gradient 6 1 $optim.message NULL $loglik [1] -147.2205 $vcov NULL > mledist(x1,"norm", fix.arg=list(mean=1/2)) $estimate sd 1.136894 $convergence [1] 0 $value [1] 1.547238 $hessian sd sd 1.547365 $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 10 6 $optim.message NULL $loglik [1] -154.7238 $vcov NULL > mledist(x1,"norm", fix.arg=list(mean=1/2), start=list(sd=1)) $estimate sd 1.1369 $convergence [1] 0 $value [1] 1.547238 $hessian sd sd 1.547324 $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 8 7 $optim.message NULL $loglik [1] -154.7238 $vcov NULL > mledist(x1,"norm", fix.arg=function(x) list(mean=0), start=list(sd=1)) $estimate sd 1.057417 $convergence [1] 0 $value [1] 1.474768 $hessian sd sd 1.788715 $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 9 6 $optim.message NULL $loglik [1] -147.4768 $vcov NULL > mledist(x1,"norm", fix.arg=function(x) list(mean=mean(x)), start=list(sd=1)) $estimate sd 1.054711 $convergence [1] 0 $value [1] 1.472205 $hessian sd sd 1.797904 $optim.function [1] "optim" $optim.method [1] "BFGS" $fix.arg $fix.arg$mean [1] 0.07560489 $fix.arg.fun function (x) list(mean = mean(x)) $weights NULL $counts function gradient 9 6 $optim.message NULL $loglik [1] -147.2205 $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 15 18 20 21 33 41 46 6 6 3 3 3 2 7 4 3 1 4 1 5 1 1 1 1 1 1 > #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.1192064 $fix.arg $fix.arg$p1 [1] 0.46 > mledist(x3, "zmgeom", fix.arg=initp1, start=list(p2=1/2))[c("estimate", "fix.arg")] $estimate p2 0.09021387 $fix.arg $fix.arg$p1 [1] 0.41 > mledist(x4, "zmgeom", fix.arg=initp1, start=list(p2=1/2))[c("estimate", "fix.arg")] $estimate p2 0.1216453 $fix.arg $fix.arg$p1 [1] 0.23 > > > #fitted coef around 0.09044335, fitted loglik -298.5162 > > > # (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] 654431974 $hessian mean sd mean 0.1192093 -9470.478 sd -9470.4777002 513938165.426 $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] -65443197415 $vcov NULL > > > # (18) management of bounds in optim/constrOptim > x <- rexp(100) > mledist(x, "exp") #optim, BFGS $estimate rate 1.091208 $convergence [1] 0 $value [1] 0.9127151 $hessian rate rate 0.8398196 $optim.function [1] "optim" $optim.method [1] "BFGS" $fix.arg NULL $fix.arg.fun NULL $weights NULL $counts function gradient 5 1 $optim.message NULL $loglik [1] -91.27151 $vcov NULL > mledist(x, "exp", optim.method="Brent", lower=0, upper=100) #optim, Brent $estimate rate 1.091208 $convergence [1] 0 $value [1] 0.9127151 $hessian [,1] [1,] 0.8398196 $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] -91.27151 $vcov NULL > mledist(x, "exp", optim.method="Nelder-Mead") #optim, Nelder-Mead $estimate rate 1.091208 $convergence [1] 0 $value [1] 0.9127151 $hessian rate rate 0.8398196 $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] -91.27151 $vcov NULL > mledist(x, "exp", lower=0, optim.method="Nelder-Mead") #constrOptim, Nelder-Mead $estimate rate 1.091208 $convergence [1] 0 $value [1] 0.9127151 $hessian rate rate 0.8398196 $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] -91.27151 $vcov NULL > mledist(x, "exp", lower=0, optim.method="BFGS") #optim, L-BFGS-B $estimate rate 1.091208 $convergence [1] 0 $value [1] 0.9127151 $hessian rate rate 0.8398196 $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] -91.27151 $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.522306 2.412763 $convergence [1] 0 $value [1] -0.1721336 $hessian shape1 shape2 shape1 0.6275431 -0.2891163 shape2 -0.2891163 0.2227355 $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] 17.21336 $vcov NULL > mledist(x, "beta", lower=0, optim.method="Nelder-Mead") #constrOptim, Nelder-Mead $estimate shape1 shape2 1.522454 2.412955 $convergence [1] 0 $value [1] -0.1721336 $hessian shape1 shape2 shape1 0.6132140 -0.2861454 shape2 -0.2861454 0.2221141 $optim.function [1] "constrOptim" $optim.method [1] "Nelder-Mead" $fix.arg NULL $fix.arg.fun NULL $weights NULL $counts function gradient 96 NA $optim.message NULL $loglik [1] 17.21336 $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.522396 2.412903 $convergence [1] 0 $value [1] -0.1721336 $hessian shape1 shape2 shape1 0.6274904 -0.2890972 shape2 -0.2890972 0.2227186 $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] 17.21336 $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. > > > # (19) large sample size issue > > if(FALSE) + { + obs <- rlnorm(1e7, 3, 2) + for(i in 2:7) + cat(i, try(mledist(obs[1:10^i], "lnorm", control=list(trace=0, REPORT=1))$estimate, silent=TRUE), "\n") + + # 2 3.143734 1.819549 + # 3 3.023688 1.955416 + # 4 2.995646 1.993392 + # 5 3.002682 2.000394 + # 6 2.999036 1.999887 + # 7 3.000222 1.999981 + } > > # (20) threshold parameter leading to infinite log-density > > dsexp <- function(x, rate, shift) + dexp(x-shift, rate=rate) > psexp <- function(x, rate, shift) + pexp(x-shift, rate=rate) > rsexp <- function(n, rate, shift) + rexp(n, rate=rate)+shift > x <- rsexp(1000, 1/4, 1) > sum(log(dsexp(x, 1, 1.995))) [1] -Inf > mledist(x, "sexp", start=list(rate=1, shift=0), upper= c(Inf, min(x)), control=list(trace=1, REPORT=1), optim="L-BFGS-B") iter 1 value 2.645185 $estimate [1] NA NA $convergence [1] 100 $loglik [1] NA $hessian [1] NA $optim.function [1] "optim" $fix.arg NULL $optim.method [1] "L-BFGS-B" $fix.arg.fun NULL $counts [1] NA NA $vcov NULL > mledist(x, "sexp", start=list(rate=1, shift=0), upper= c(Inf, min(x)), control=list(trace=1, REPORT=1)) Nelder-Mead direct search function minimizer function value for initial parameters = 4.916798 Scaled convergence tolerance is 7.3266e-08 Stepsize computed as 0.100000 BUILD 3 5.313168 4.816799 EXTENSION 5 4.916798 4.036583 EXTENSION 7 4.816799 3.535943 EXTENSION 9 4.036583 2.468648 LO-REDUCTION 11 3.535943 2.468648 HI-REDUCTION 13 2.881002 2.468648 HI-REDUCTION 15 2.634853 2.468648 LO-REDUCTION 17 2.497541 2.451760 HI-REDUCTION 19 2.469930 2.451760 HI-REDUCTION 21 2.468648 2.451760 EXTENSION 23 2.462041 2.431919 LO-REDUCTION 25 2.451760 2.431919 EXTENSION 27 2.438661 2.406387 EXTENSION 29 2.431919 2.373719 HI-REDUCTION 31 2.409903 2.373719 HI-REDUCTION 33 2.406387 2.373719 REFLECTION 35 2.399498 2.370311 HI-REDUCTION 37 2.385452 2.370311 HI-REDUCTION 39 2.378518 2.370311 HI-REDUCTION 41 2.375079 2.370311 REFLECTION 43 2.373719 2.368374 EXTENSION 45 2.370311 2.365356 HI-REDUCTION 47 2.368374 2.365356 LO-REDUCTION 49 2.367957 2.365356 HI-REDUCTION 51 2.366505 2.365356 LO-REDUCTION 53 2.366502 2.365356 HI-REDUCTION 55 2.366004 2.365356 HI-REDUCTION 57 2.365871 2.365356 EXTENSION 59 2.365773 2.365301 HI-REDUCTION 61 2.365473 2.365301 HI-REDUCTION 63 2.365358 2.365301 REFLECTION 65 2.365356 2.365233 HI-REDUCTION 67 2.365301 2.365233 HI-REDUCTION 69 2.365265 2.365233 LO-REDUCTION 71 2.365256 2.365233 HI-REDUCTION 73 2.365242 2.365233 LO-REDUCTION 75 2.365242 2.365233 HI-REDUCTION 77 2.365236 2.365233 HI-REDUCTION 79 2.365235 2.365233 REFLECTION 81 2.365234 2.365233 HI-REDUCTION 83 2.365233 2.365233 HI-REDUCTION 85 2.365233 2.365232 HI-REDUCTION 87 2.365233 2.365232 Exiting from Nelder Mead minimizer 89 function evaluations used Nelder-Mead direct search function minimizer function value for initial parameters = 2.364447 Scaled convergence tolerance is 3.5233e-08 Stepsize computed as 0.100285 BUILD 3 99999999999999996864064668260046802.000000 2.364447 LO-REDUCTION 5 2.425892 2.364447 LO-REDUCTION 7 2.417420 2.364447 HI-REDUCTION 9 2.385306 2.364447 LO-REDUCTION 11 2.374596 2.364447 HI-REDUCTION 13 2.369944 2.364447 HI-REDUCTION 15 2.368844 2.364447 LO-REDUCTION 17 2.367970 2.364447 HI-REDUCTION 19 2.366470 2.364447 HI-REDUCTION 21 2.365808 2.364447 LO-REDUCTION 23 2.365764 2.364447 HI-REDUCTION 25 2.365133 2.364447 HI-REDUCTION 27 2.364829 2.364447 HI-REDUCTION 29 2.364680 2.364447 REFLECTION 31 2.364674 2.364388 HI-REDUCTION 33 2.364527 2.364388 HI-REDUCTION 35 2.364464 2.364388 REFLECTION 37 2.364447 2.364358 HI-REDUCTION 39 2.364404 2.364358 HI-REDUCTION 41 2.364388 2.364358 LO-REDUCTION 43 2.364384 2.364358 REFLECTION 45 2.364375 2.364358 HI-REDUCTION 47 2.364366 2.364358 REFLECTION 49 2.364358 2.364347 HI-REDUCTION 51 2.364358 2.364347 HI-REDUCTION 53 2.364354 2.364347 HI-REDUCTION 55 2.364353 2.364347 HI-REDUCTION 57 2.364352 2.364347 HI-REDUCTION 59 2.364351 2.364347 REFLECTION 61 2.364350 2.364346 HI-REDUCTION 63 2.364348 2.364346 HI-REDUCTION 65 2.364347 2.364346 HI-REDUCTION 67 2.364347 2.364346 HI-REDUCTION 69 2.364347 2.364346 HI-REDUCTION 71 2.364346 2.364346 HI-REDUCTION 73 2.364346 2.364346 HI-REDUCTION 75 2.364346 2.364346 HI-REDUCTION 77 2.364346 2.364346 HI-REDUCTION 79 2.364346 2.364346 HI-REDUCTION 81 2.364346 2.364346 HI-REDUCTION 83 2.364346 2.364346 HI-REDUCTION 85 2.364346 2.364346 Exiting from Nelder Mead minimizer 87 function evaluations used Nelder-Mead direct search function minimizer function value for initial parameters = 2.364345 Scaled convergence tolerance is 3.52315e-08 Stepsize computed as 0.100325 BUILD 3 99999999999999996864064668260046802.000000 2.364345 LO-REDUCTION 5 2.425810 2.364345 LO-REDUCTION 7 2.417336 2.364345 HI-REDUCTION 9 2.385209 2.364345 LO-REDUCTION 11 2.374498 2.364345 HI-REDUCTION 13 2.369845 2.364345 HI-REDUCTION 15 2.368748 2.364345 LO-REDUCTION 17 2.367871 2.364345 HI-REDUCTION 19 2.366370 2.364345 HI-REDUCTION 21 2.365706 2.364345 LO-REDUCTION 23 2.365662 2.364345 HI-REDUCTION 25 2.365032 2.364345 HI-REDUCTION 27 2.364728 2.364345 HI-REDUCTION 29 2.364579 2.364345 HI-REDUCTION 31 2.364573 2.364345 LO-REDUCTION 33 2.364505 2.364345 HI-REDUCTION 35 2.364429 2.364345 HI-REDUCTION 37 2.364396 2.364345 HI-REDUCTION 39 2.364392 2.364345 HI-REDUCTION 41 2.364382 2.364345 HI-REDUCTION 43 2.364377 2.364345 HI-REDUCTION 45 2.364371 2.364345 HI-REDUCTION 47 2.364367 2.364345 HI-REDUCTION 49 2.364364 2.364345 HI-REDUCTION 51 2.364361 2.364345 HI-REDUCTION 53 2.364358 2.364345 HI-REDUCTION 55 2.364356 2.364345 HI-REDUCTION 57 2.364355 2.364345 HI-REDUCTION 59 2.364353 2.364345 HI-REDUCTION 61 2.364352 2.364345 HI-REDUCTION 63 2.364351 2.364345 HI-REDUCTION 65 2.364350 2.364345 HI-REDUCTION 67 2.364349 2.364345 HI-REDUCTION 69 2.364349 2.364345 HI-REDUCTION 71 2.364348 2.364345 HI-REDUCTION 73 2.364348 2.364345 HI-REDUCTION 75 2.364347 2.364345 HI-REDUCTION 77 2.364347 2.364345 HI-REDUCTION 79 2.364347 2.364345 HI-REDUCTION 81 2.364346 2.364345 HI-REDUCTION 83 2.364346 2.364345 HI-REDUCTION 85 2.364346 2.364345 HI-REDUCTION 87 2.364346 2.364345 HI-REDUCTION 89 2.364346 2.364345 HI-REDUCTION 91 2.364346 2.364345 HI-REDUCTION 93 2.364346 2.364345 HI-REDUCTION 95 2.364346 2.364345 HI-REDUCTION 97 2.364346 2.364345 HI-REDUCTION 99 2.364346 2.364345 HI-REDUCTION 101 2.364345 2.364345 HI-REDUCTION 103 2.364345 2.364345 HI-REDUCTION 105 2.364345 2.364345 HI-REDUCTION 107 2.364345 2.364345 HI-REDUCTION 109 2.364345 2.364345 HI-REDUCTION 111 2.364345 2.364345 HI-REDUCTION 113 2.364345 2.364345 HI-REDUCTION 115 2.364345 2.364345 HI-REDUCTION 117 2.364345 2.364345 HI-REDUCTION 119 2.364345 2.364345 REFLECTION 121 2.364345 2.364345 Exiting from Nelder Mead minimizer 123 function evaluations used $estimate rate shift 0.2555494 1.0032457 $convergence [1] 0 $value [1] 2.364446 $hessian rate shift rate 1.000002 -1.000000e+00 shift -1.000000 2.220446e-10 $optim.function [1] "constrOptim" $optim.method [1] "Nelder-Mead" $fix.arg NULL $fix.arg.fun NULL $weights NULL $counts function gradient 176 NA $optim.message NULL $loglik [1] -2364.446 $vcov NULL > > #fitted coef is 0.243 1.003 , fitted loglik -2415, fitted hessian is > # rate shift > # rate 1000 -1.00e+03 > # shift -1000 2.14e-05 > > > > > > proc.time() user system elapsed 3.76 0.45 4.15