From 4fa8de2761a02a469610fa95ac3244fc9d6a3118 Mon Sep 17 00:00:00 2001 From: Christophe Dutang Date: Fri, 25 Oct 2024 16:27:14 +0200 Subject: [PATCH] scale the logLik --- NEWS.md | 1 + R/mledist.R | 39 ++--- tests/t-mledist-nocens.R | 316 ++++++++++++++++++++++----------------- 3 files changed, 196 insertions(+), 160 deletions(-) diff --git a/NEWS.md b/NEWS.md index 39af747..628f82e 100644 --- a/NEWS.md +++ b/NEWS.md @@ -9,6 +9,7 @@ BUG FIX - mgedist() may suffer a numerical issue for Anderson-Darling GoF metrics. All GoF metrics now take care of numerical issue, such as log(0) or 1/0, and are properly scaled by the sample sized to avoid large sample size issues. Thanks for Ethan Chapman for reporting the bug. - the default starting value for the gamma distribution was wrongly computed for the rate parameter. Thanks for Wendy Martin for reporting the bug. +- mledist() may suffer a scaling issue and the objective function is properly scaled by the sample sized to avoid large sample size issues. # fitdistrplus 1.2-1 diff --git a/R/mledist.R b/R/mledist.R index bbd2c8d..ff2cc38 100644 --- a/R/mledist.R +++ b/R/mledist.R @@ -64,6 +64,7 @@ mledist <- function (data, distr, start=NULL, fix.arg=NULL, optim.method="defaul if (!(is.numeric(data) & length(data)>1)) stop(txt1uncens) checkUncensoredNAInfNan(data) + data.size <- length(data) } else { cens <- TRUE @@ -75,7 +76,7 @@ mledist <- function (data, distr, start=NULL, fix.arg=NULL, optim.method="defaul pdistname<-paste("p", distname, sep="") if (!exists(pdistname, mode="function")) stop(paste("The ", pdistname, " function must be defined to apply maximum likelihood to censored data")) - + data.size <- NROW(censdata) } if (cens) { @@ -153,7 +154,7 @@ mledist <- function (data, distr, start=NULL, fix.arg=NULL, optim.method="defaul ############# MLE fit using optim or custom.optim ########## - # definition of the function to minimize : - log likelihood + # definition of the function to minimize : minus average log likelihood # for non censored data if (!cens && is.null(weights)) { # the argument names are: @@ -163,12 +164,12 @@ mledist <- function (data, distr, start=NULL, fix.arg=NULL, optim.method="defaul # - ddistnam for distribution name if ("log" %in% argddistname){ fnobj <- function(par, fix.arg, obs, ddistnam){ - -sum(do.call(ddistnam, c(list(obs), as.list(par), as.list(fix.arg), log=TRUE) ) ) + -mean(do.call(ddistnam, c(list(obs), as.list(par), as.list(fix.arg), log=TRUE) ) ) } } else{ fnobj <- function(par, fix.arg, obs, ddistnam) { - -sum(log(do.call(ddistnam, c(list(obs), as.list(par), as.list(fix.arg)) ) ) ) + -mean(log(do.call(ddistnam, c(list(obs), as.list(par), as.list(fix.arg)) ) ) ) } } } @@ -177,22 +178,22 @@ mledist <- function (data, distr, start=NULL, fix.arg=NULL, optim.method="defaul argpdistname<-names(formals(pdistname)) if (("log" %in% argddistname) & ("log.p" %in% argpdistname)) fnobjcens <- function(par, fix.arg, rcens, lcens, icens, ncens, ddistnam, pdistnam) - -sum(do.call(ddistnam, c(list(ncens), as.list(par), as.list(fix.arg), list(log=TRUE)))) - - sum(do.call(pdistnam, c(list(lcens), as.list(par), as.list(fix.arg), list(log=TRUE)))) - - sum(do.call(pdistnam, c(list(rcens), as.list(par), as.list(fix.arg), list(lower.tail=FALSE), list(log=TRUE)))) - - sum(log(do.call(pdistnam, c(list(icens$right), as.list(par), as.list(fix.arg))) - # without log=TRUE here + -mean(do.call(ddistnam, c(list(ncens), as.list(par), as.list(fix.arg), list(log=TRUE)))) - + mean(do.call(pdistnam, c(list(lcens), as.list(par), as.list(fix.arg), list(log=TRUE)))) - + mean(do.call(pdistnam, c(list(rcens), as.list(par), as.list(fix.arg), list(lower.tail=FALSE), list(log=TRUE)))) - + mean(log(do.call(pdistnam, c(list(icens$right), as.list(par), as.list(fix.arg))) - # without log=TRUE here do.call(pdistnam, c(list(icens$left), as.list(par), as.list(fix.arg))) )) # without log=TRUE here else fnobjcens <- function(par, fix.arg, rcens, lcens, icens, ncens, ddistnam, pdistnam) - -sum(log(do.call(ddistnam, c(list(ncens), as.list(par), as.list(fix.arg))))) - - sum(log(do.call(pdistnam, c(list(lcens), as.list(par), as.list(fix.arg))))) - - sum(log(1-do.call(pdistnam, c(list(rcens), as.list(par), as.list(fix.arg))))) - - sum(log(do.call(pdistnam, c(list(icens$right), as.list(par), as.list(fix.arg))) - + -mean(log(do.call(ddistnam, c(list(ncens), as.list(par), as.list(fix.arg))))) - + mean(log(do.call(pdistnam, c(list(lcens), as.list(par), as.list(fix.arg))))) - + mean(log(1-do.call(pdistnam, c(list(rcens), as.list(par), as.list(fix.arg))))) - + mean(log(do.call(pdistnam, c(list(icens$right), as.list(par), as.list(fix.arg))) - do.call(pdistnam, c(list(icens$left), as.list(par), as.list(fix.arg))) )) }else if(!cens && !is.null(weights)) { fnobj <- function(par, fix.arg, obs, ddistnam) { - -sum(weights * log(do.call(ddistnam, c(list(obs), as.list(par), as.list(fix.arg)) ) ) ) + -mean(weights * log(do.call(ddistnam, c(list(obs), as.list(par), as.list(fix.arg)) ) ) ) } }else if(cens && !is.null(weights)) { @@ -203,10 +204,10 @@ mledist <- function (data, distr, start=NULL, fix.arg=NULL, optim.method="defaul p3 <- log(1-do.call(pdistnam, c(list(rcens), as.list(par), as.list(fix.arg)))) p4 <- log(do.call(pdistnam, c(list(icens$right), as.list(par), as.list(fix.arg))) - do.call(pdistnam, c(list(icens$left), as.list(par), as.list(fix.arg))) ) - - sum(weights[irow.ncens] * p1) - - sum(weights[irow.lcens] * p2) - - sum(weights[irow.rcens] * p3) - - sum(weights[irow.icens] * p4) + - mean(weights[irow.ncens] * p1) - + mean(weights[irow.lcens] * p2) - + mean(weights[irow.rcens] * p3) - + mean(weights[irow.icens] * p4) } } @@ -385,7 +386,7 @@ mledist <- function (data, distr, start=NULL, fix.arg=NULL, optim.method="defaul res <- list(estimate = opt$par, convergence = opt$convergence, value=opt$value, hessian = opt$hessian, optim.function=opt.fun, optim.method=meth, fix.arg = fix.arg, fix.arg.fun = fix.arg.fun, weights = weights, - counts=opt$counts, optim.message=opt$message, loglik = -opt$value, + counts=opt$counts, optim.message=opt$message, loglik = -opt$value*data.size, vcov=varcovar) } else # Try to minimize the minus (log-)likelihood using a user-supplied optim function @@ -428,7 +429,7 @@ mledist <- function (data, distr, start=NULL, fix.arg=NULL, optim.method="defaul res <- list(estimate = opt$par, convergence = opt$convergence, value=opt$value, hessian = opt$hessian, optim.function = custom.optim, optim.method = method.cust, fix.arg = fix.arg, fix.arg.fun = fix.arg.fun, weights = weights, - counts=opt$counts, optim.message=opt$message, loglik = -opt$value, + counts=opt$counts, optim.message=opt$message, loglik = -opt$value*data.size, vcov=varcovar) } diff --git a/tests/t-mledist-nocens.R b/tests/t-mledist-nocens.R index 6acb418..bd42bd0 100644 --- a/tests/t-mledist-nocens.R +++ b/tests/t-mledist-nocens.R @@ -10,6 +10,8 @@ set.seed(1234) x1 <- rnorm(n=100) mledist(x1,"norm") +#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 @@ -17,6 +19,9 @@ 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) mledist(x1,"gumbel",start=list(a=10,b=5), silent=FALSE) +#fitted coef around -0.6267919 0.8564855, fitted loglik -139.3812 + + # (3) fit a discrete distribution (Poisson) # @@ -24,6 +29,8 @@ set.seed(1234) x2 <- rpois(n=30,lambda = 2) mledist(x2,"pois") +#fitted coef around 1.7, fitted loglik -46.18434 + # (4) fit a finite-support distribution (beta) # @@ -31,6 +38,7 @@ set.seed(1234) x3 <- rbeta(n=100,shape1=5, shape2=10) mledist(x3,"beta") +#fitted coef around 4.859798 10.918841, fitted loglik 78.33052 # (5) fit frequency distributions on USArrests dataset. # @@ -40,6 +48,8 @@ mledist(x4, "pois", silent=TRUE) mledist(x4, "pois", silent=FALSE) mledist(x4, "nbinom") +#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), @@ -49,35 +59,35 @@ mledist(x4, "nbinom") 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") - + cat(i, try(mledist(x*10^i, "cauchy")$estimate, silent=TRUE), "\n") + # (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) + -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") + cat(i, try(mledist(x*10^i, "cauchy")$estimate, silent=TRUE), "\n") # (8) normal mixture @@ -86,49 +96,50 @@ cat(i, try(mledist(x*10^i, "cauchy")$estimate, silent=TRUE), "\n") #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) + 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) + 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)) - - - + lower=c(0, 0, 0, 0, 0)) +#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) - - -#simulate a sample -x4 <- rllogis(1000, 6, 2) - -#fit -mledist(x4, "llogis", start=list(shape=10, rate=1), lower=1, upper=Inf) - + 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 } @@ -137,54 +148,54 @@ mledist(x4, "llogis", start=list(shape=10, rate=1), lower=1, upper=Inf) # 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) - + + + 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) + } @@ -192,54 +203,54 @@ fgenoud=fgenoud$estimate) # 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) - - + + + 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) + + } @@ -279,6 +290,8 @@ f2 <- mledist(xval, "pois", weights=xtab, start=list(lambda=mean(x))) f1$estimate f2$estimate #should be identical +#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)))) #test non integer weights @@ -352,6 +365,9 @@ mledist(x3, "zmgeom", fix.arg=initp1, start=list(p2=1/2))[c("estimate", "fix.arg mledist(x4, "zmgeom", fix.arg=initp1, start=list(p2=1/2))[c("estimate", "fix.arg")] +#fitted coef around 0.09044335, fitted loglik -298.5162 + + # (17) test the component optim.message x <- rnorm(100) #change parameter to obtain unsuccessful convergence @@ -373,3 +389,21 @@ mledist(x, "beta", lower=0, optim.method="Nelder-Mead") #constrOptim, Nelder-Mea #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 + + +# (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 +}