R Under development (unstable) (2025-01-06 r87534 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 > > 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] 1.418309 $hessian mean sd mean 1.00126 0.000000 sd 0.00000 2.002538 $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 > > #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 3 9.75 2 9.5 4 8 4.5 6.625 1.5 6.875 -0.25 5.6875 2.25 2.8125 2.375 -0.53125 -2.5 1.875 -6 -0.5 $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.348589 Scaled convergence tolerance is 1.09503e-07 Stepsize computed as 1.000000 11 5 10 6 BUILD 3 8.869361 5.605410 9 6 8 6.5 EXTENSION 5 7.348589 4.163219 8 7.5 7 8.75 EXTENSION 7 5.605410 3.629945 5 9.25 2.5 10.875 REFLECTION 9 4.163219 3.422121 4 11.5 5 10.25 LO-REDUCTION 11 3.629945 3.422121 3 10.75 4 10.25 LO-REDUCTION 13 3.484496 3.422121 3 9.75 2 9.5 EXTENSION 15 3.427361 3.285277 4 8 4.5 6.625 EXTENSION 17 3.422121 3.228640 1.5 6.875 REFLECTION 19 3.285277 2.972002 4 4 3.5 5.375 LO-REDUCTION 21 3.228640 2.972002 0.5 5.625 REFLECTION 23 3.007904 2.752000 2.25 5.8125 HI-REDUCTION 25 2.972002 2.752000 1.25 4.5625 1.125 3.40625 EXTENSION 27 2.879777 2.366250 1.53125 5.164062 HI-REDUCTION 29 2.752000 2.366250 2.15625 2.945312 1.742188 3.615234 LO-REDUCTION 31 2.726103 2.366250 1.335938 1.857422 1.238281 0.2041016 $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) > # > > set.seed(1234) > x2 <- rpois(n=30,lambda = 2) > mledist(x2,"pois") $estimate lambda 1.7 $convergence [1] 0 $value [1] 1.539478 $hessian lambda lambda 0.5882357 $optim.function [1] "optim" $optim.method [1] "BFGS" $fix.arg NULL $fix.arg.fun NULL $weights NULL $counts function gradient 4 1 $optim.message NULL $loglik [1] -46.18434 $vcov NULL > > #fitted coef around 1.7, fitted loglik -46.18434 > > # (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] -0.7833052 $hessian shape1 shape2 shape1 0.16295311 -0.06542752 shape2 -0.06542752 0.03047900 $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 > > #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.005856175 $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.005856175 $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.987921e-07 mu -3.987921e-07 1.282596e-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. > > 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 -0.0002897257 0.001193096 > > > # (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)) > fit1 $estimate w m1 s1 m2 s2 9.394775e+05 7.863282e+00 2.661790e+00 1.240012e+02 2.928074e+00 $convergence [1] 7 $value [1] -11.32679 $hessian w m1 s1 m2 s2 w 2.00030720 -0.162971396 0.03754193 -0.32871928 0.293953944 m1 -0.16297140 0.049555201 0.02988310 -0.01307104 0.003907968 s1 0.03754193 0.029883098 -0.05285673 0.03239808 -0.032171477 m2 -0.32871928 -0.013071041 0.03239808 0.11732246 0.174422901 s2 0.29395394 0.003907968 -0.03217148 0.17442290 0.516194874 $optim.function [1] "constrOptim" $optim.method [1] "Nelder-Mead" $fix.arg NULL $fix.arg.fun NULL $weights NULL $counts function gradient 20914 NA $optim.message [1] "Barrier algorithm ran out of iterations and did not converge" $loglik [1] 22653.58 $vcov NULL > > #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.108684 2.001231 $convergence [1] 0 $value [1] -0.503337 $hessian shape rate shape 0.015023006 -0.006843682 rate -0.006843682 6.108444812 $optim.function [1] "constrOptim" $optim.method [1] "Nelder-Mead" $fix.arg NULL $fix.arg.fun NULL $weights NULL $counts function gradient 112 NA $optim.message NULL $loglik [1] 503.337 $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 9.87 > f2$estimate #should be identical lambda 9.87 > > #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 20.19 $convergence [1] 1 $value [1] 14.9111 $hessian lambda lambda 0.04952947 $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] 1.418309 $hessian mean sd mean 1.00126 0.000000 sd 0.00000 2.002538 $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] 1.418309 $hessian mean sd mean 1.0013830523 0.0002319865 sd 0.0002319865 2.0031537814 $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] 1.418309 $hessian sd sd 2.002538 $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 5 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] 1.597803 $hessian sd sd 1.398537 $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 13 8 $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] 1.597803 $hessian sd sd 1.398536 $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 13 8 $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.011586 $convergence [1] 0 $value [1] 1.430463 $hessian sd sd 1.954498 $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 6 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.9987434 $convergence [1] 0 $value [1] 1.418309 $hessian sd sd 2.008833 $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 2 1 $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.1232226 $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.0949462 $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.09494427 $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] 654433429 $hessian mean sd mean 8.940697e-02 -9470.463 sd -9.470463e+03 513939307.988 $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] -65443342909 $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] 0.9985353 $hessian rate rate 0.9970769 $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] -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] 0.9985353 $hessian [,1] [1,] 0.9970769 $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] 0.9985353 $hessian rate rate 0.9970769 $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] 0.9985353 $hessian rate rate 0.9970769 $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] 0.9985353 $hessian rate rate 0.9970769 $optim.function [1] "optim" $optim.method [1] "L-BFGS-B" $fix.arg NULL $fix.arg.fun NULL $weights NULL $counts function gradient 11 11 $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] -0.1869748 $hessian shape1 shape2 shape1 0.5475563 -0.2661467 shape2 -0.2661467 0.2085442 $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] -0.1869748 $hessian shape1 shape2 shape1 0.5536355 -0.2663521 shape2 -0.2663521 0.2068112 $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.666900 2.568436 $convergence [1] 0 $value [1] -0.1869748 $hessian shape1 shape2 shape1 0.5475752 -0.2661527 shape2 -0.2661527 0.2085483 $optim.function [1] "optim" $optim.method [1] "L-BFGS-B" $fix.arg NULL $fix.arg.fun NULL $weights NULL $counts function gradient 6 6 $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. > > > # (19) large sample size issue > > if(FALSE) + { + set.seed(123) + 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 > > set.seed(123) > 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.709077 iter 2 value 2.598940 iter 3 value 2.594783 iter 4 value 2.431621 iter 5 value 2.420711 iter 6 value 2.415262 iter 7 value 2.415034 iter 8 value 2.415031 iter 9 value 2.415031 final value 2.415031 converged $estimate rate shift 0.2429195 1.0033039 $convergence [1] 0 $value [1] 2.415031 $hessian rate shift rate 16.946891 1.058807 shift 1.058807 -353.634854 $optim.function [1] "optim" $optim.method [1] "L-BFGS-B" $fix.arg NULL $fix.arg.fun NULL $weights NULL $counts function gradient 20 20 $optim.message [1] "CONVERGENCE: REL_REDUCTION_OF_F <= FACTR*EPSMCH" $loglik [1] -2415.031 $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 = 5.119917 Scaled convergence tolerance is 7.62927e-08 Stepsize computed as 0.100000 BUILD 3 5.536598 5.019917 EXTENSION 5 5.119917 4.199078 EXTENSION 7 5.019917 3.678126 EXTENSION 9 4.199078 2.519428 LO-REDUCTION 11 3.678126 2.519428 HI-REDUCTION 13 2.860618 2.519428 HI-REDUCTION 15 2.586843 2.519428 HI-REDUCTION 17 2.543317 2.508438 HI-REDUCTION 19 2.519428 2.488134 EXTENSION 21 2.508438 2.456584 HI-REDUCTION 23 2.488134 2.456584 EXTENSION 25 2.479061 2.416758 HI-REDUCTION 27 2.456584 2.416758 LO-REDUCTION 29 2.455083 2.416758 HI-REDUCTION 31 2.439583 2.416758 HI-REDUCTION 33 2.433202 2.416758 LO-REDUCTION 35 2.432005 2.416758 HI-REDUCTION 37 2.424344 2.416758 HI-REDUCTION 39 2.420647 2.416758 HI-REDUCTION 41 2.418847 2.416758 HI-REDUCTION 43 2.417969 2.416758 HI-REDUCTION 45 2.417923 2.416758 LO-REDUCTION 47 2.417540 2.416758 LO-REDUCTION 49 2.417069 2.416758 EXTENSION 51 2.416765 2.416081 REFLECTION 53 2.416758 2.416072 HI-REDUCTION 55 2.416345 2.416072 REFLECTION 57 2.416081 2.415832 HI-REDUCTION 59 2.416072 2.415832 HI-REDUCTION 61 2.415939 2.415832 HI-REDUCTION 63 2.415916 2.415832 HI-REDUCTION 65 2.415879 2.415832 HI-REDUCTION 67 2.415859 2.415832 HI-REDUCTION 69 2.415845 2.415832 HI-REDUCTION 71 2.415833 2.415827 HI-REDUCTION 73 2.415832 2.415821 HI-REDUCTION 75 2.415827 2.415817 LO-REDUCTION 77 2.415821 2.415816 HI-REDUCTION 79 2.415818 2.415816 HI-REDUCTION 81 2.415817 2.415816 LO-REDUCTION 83 2.415816 2.415816 EXTENSION 85 2.415816 2.415815 HI-REDUCTION 87 2.415816 2.415815 EXTENSION 89 2.415816 2.415815 LO-REDUCTION 91 2.415815 2.415815 EXTENSION 93 2.415815 2.415814 EXTENSION 95 2.415815 2.415813 REFLECTION 97 2.415814 2.415813 EXTENSION 99 2.415813 2.415813 LO-REDUCTION 101 2.415813 2.415813 LO-REDUCTION 103 2.415813 2.415813 HI-REDUCTION 105 2.415813 2.415813 LO-REDUCTION 107 2.415813 2.415813 Exiting from Nelder Mead minimizer 109 function evaluations used Nelder-Mead direct search function minimizer function value for initial parameters = 2.415030 Scaled convergence tolerance is 3.59868e-08 Stepsize computed as 0.100289 BUILD 3 99999999999999996864064668260046802.000000 2.415030 LO-REDUCTION 5 2.482311 2.415030 LO-REDUCTION 7 2.470904 2.415030 HI-REDUCTION 9 2.436529 2.415030 LO-REDUCTION 11 2.424735 2.415030 HI-REDUCTION 13 2.420256 2.415030 HI-REDUCTION 15 2.419485 2.415030 LO-REDUCTION 17 2.418397 2.415030 HI-REDUCTION 19 2.416955 2.415030 HI-REDUCTION 21 2.416348 2.415030 LO-REDUCTION 23 2.416283 2.415030 HI-REDUCTION 25 2.415683 2.415030 HI-REDUCTION 27 2.415396 2.415030 HI-REDUCTION 29 2.415264 2.415030 LO-REDUCTION 31 2.415255 2.415030 HI-REDUCTION 33 2.415160 2.415030 REFLECTION 35 2.415132 2.415015 HI-REDUCTION 37 2.415076 2.415015 REFLECTION 39 2.415030 2.414963 HI-REDUCTION 41 2.415015 2.414963 LO-REDUCTION 43 2.415007 2.414963 HI-REDUCTION 45 2.414985 2.414963 REFLECTION 47 2.414965 2.414940 LO-REDUCTION 49 2.414963 2.414940 HI-REDUCTION 51 2.414951 2.414940 HI-REDUCTION 53 2.414946 2.414940 REFLECTION 55 2.414942 2.414934 LO-REDUCTION 57 2.414940 2.414934 HI-REDUCTION 59 2.414938 2.414934 REFLECTION 61 2.414937 2.414932 HI-REDUCTION 63 2.414935 2.414932 HI-REDUCTION 65 2.414934 2.414932 LO-REDUCTION 67 2.414934 2.414932 REFLECTION 69 2.414933 2.414931 HI-REDUCTION 71 2.414932 2.414931 LO-REDUCTION 73 2.414932 2.414931 HI-REDUCTION 75 2.414932 2.414931 HI-REDUCTION 77 2.414931 2.414931 REFLECTION 79 2.414931 2.414931 HI-REDUCTION 81 2.414931 2.414931 HI-REDUCTION 83 2.414931 2.414931 HI-REDUCTION 85 2.414931 2.414931 REFLECTION 87 2.414931 2.414931 Exiting from Nelder Mead minimizer 89 function evaluations used Nelder-Mead direct search function minimizer function value for initial parameters = 2.414931 Scaled convergence tolerance is 3.59853e-08 Stepsize computed as 0.100330 BUILD 3 99999999999999996864064668260046802.000000 2.414931 LO-REDUCTION 5 2.482228 2.414931 LO-REDUCTION 7 2.470819 2.414931 HI-REDUCTION 9 2.436432 2.414931 LO-REDUCTION 11 2.424638 2.414931 HI-REDUCTION 13 2.420159 2.414931 HI-REDUCTION 15 2.419392 2.414931 LO-REDUCTION 17 2.418300 2.414931 HI-REDUCTION 19 2.416856 2.414931 HI-REDUCTION 21 2.416247 2.414931 LO-REDUCTION 23 2.416183 2.414931 HI-REDUCTION 25 2.415583 2.414931 HI-REDUCTION 27 2.415296 2.414931 HI-REDUCTION 29 2.415165 2.414931 LO-REDUCTION 31 2.415156 2.414931 HI-REDUCTION 33 2.415060 2.414931 HI-REDUCTION 35 2.415032 2.414931 HI-REDUCTION 37 2.415018 2.414931 HI-REDUCTION 39 2.415002 2.414931 HI-REDUCTION 41 2.414992 2.414931 HI-REDUCTION 43 2.414981 2.414931 HI-REDUCTION 45 2.414974 2.414931 HI-REDUCTION 47 2.414967 2.414931 HI-REDUCTION 49 2.414961 2.414931 HI-REDUCTION 51 2.414956 2.414931 HI-REDUCTION 53 2.414952 2.414931 HI-REDUCTION 55 2.414949 2.414931 HI-REDUCTION 57 2.414946 2.414931 HI-REDUCTION 59 2.414944 2.414931 HI-REDUCTION 61 2.414942 2.414931 HI-REDUCTION 63 2.414940 2.414931 HI-REDUCTION 65 2.414938 2.414931 HI-REDUCTION 67 2.414937 2.414931 HI-REDUCTION 69 2.414936 2.414931 HI-REDUCTION 71 2.414935 2.414931 HI-REDUCTION 73 2.414934 2.414931 HI-REDUCTION 75 2.414934 2.414931 HI-REDUCTION 77 2.414933 2.414931 HI-REDUCTION 79 2.414933 2.414931 HI-REDUCTION 81 2.414933 2.414931 HI-REDUCTION 83 2.414932 2.414931 HI-REDUCTION 85 2.414932 2.414931 HI-REDUCTION 87 2.414932 2.414931 HI-REDUCTION 89 2.414932 2.414931 HI-REDUCTION 91 2.414931 2.414931 HI-REDUCTION 93 2.414931 2.414931 HI-REDUCTION 95 2.414931 2.414931 HI-REDUCTION 97 2.414931 2.414931 HI-REDUCTION 99 2.414931 2.414931 HI-REDUCTION 101 2.414931 2.414931 HI-REDUCTION 103 2.414931 2.414931 REFLECTION 105 2.414931 2.414930 HI-REDUCTION 107 2.414931 2.414930 HI-REDUCTION 109 2.414931 2.414930 HI-REDUCTION 111 2.414931 2.414930 HI-REDUCTION 113 2.414931 2.414930 REFLECTION 115 2.414931 2.414930 Exiting from Nelder Mead minimizer 117 function evaluations used $estimate rate shift 0.2429235 1.0033039 $convergence [1] 0 $value [1] 2.415031 $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 198 NA $optim.message NULL $loglik [1] -2415.031 $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 10.73 0.35 11.09