Wrapper function to compute maximum likelihood estimators (MLE)
of any distribution implemented in R
.
Usage
maxlogL(
x,
dist = "dnorm",
fixed = NULL,
link = NULL,
start = NULL,
lower = NULL,
upper = NULL,
optimizer = "nlminb",
control = NULL,
StdE_method = c("optim", "numDeriv"),
silent = FALSE,
...
)
Arguments
- x
a vector with data to be fitted. This argument must be a matrix with hierarchical distributions.
- dist
a length-one character vector with the name of density/mass function of interest. The default value is
'dnorm'
, to compute maximum likelihood estimators of normal distribution.- fixed
a list with fixed/known parameters of distribution of interest. Fixed parameters must be passed with its name.
- link
a list with names of parameters to be linked, and names of the link function object. For names of parameters, please visit documentation of density/mass function. There are three link functions available:
log_link
,logit_link
andNegInv_link
.- start
a numeric vector with initial values for the parameters to be estimated.
- lower
a numeric vector with lower bounds, with the same length of argument
start
(for box-constrained optimization).- upper
a numeric vector with upper bounds, with the same length of argument
start
(for box-constrained optimization).- optimizer
a length-one character vector with the name of optimization routine.
nlminb
,optim
,DEoptim
andga
are available; custom optimization routines can also be implemented.nlminb
is the default routine.- control
control parameters of the optimization routine. See, e.g.,
optim
'scontrol
list,nlminb
'scontrol
list, andDEoptim.control
for DEoptim.- StdE_method
a length-one character vector with the routine for Hessian matrix computation. The This is needed for standard error estimation. The options available are
"optim"
and"numDeriv"
. For further information, visitoptim
orhessian
.- silent
logical. If TRUE, warnings of
maxlogL
are suppressed.- ...
further arguments to be supplied to the optimizer.
Value
A list with class "maxlogL"
containing the following lists:
- fit
A list with output information about estimation.
- inputs
A list with all input arguments.
- outputs
A list with some output additional information:
Number of parameters.
Sample size
Standard error computation method.
Details
maxlogL
computes the likelihood function corresponding to
the distribution specified in argument dist
and maximizes it through
optim
, nlminb
or
DEoptim
. maxlogL
generates an S3 object of class maxlogL
.
Noncentrality parameters must be named as ncp
in the distribution.
Note
The following generic functions can be used with a maxlogL
object:
summary, print, AIC, BIC, logLik
.
References
Nelder JA, Mead R (1965). “A Simplex Method for Function Minimization.” The Computer Journal, 7(4), 308–313. ISSN 0010-4620, doi:10.1093/comjnl/7.4.308 , https://academic.oup.com/comjnl/article-lookup/doi/10.1093/comjnl/7.4.308.
Fox PA, Hall AP, Schryer NL (1978). “The PORT Mathematical Subroutine Library.” ACM Transactions on Mathematical Software, 4(2), 104–126. ISSN 00983500, doi:10.1145/355780.355783 , https://dl.acm.org/doi/10.1145/355780.355783.
Nash JC (1979). Compact Numerical Methods for Computers. Linear Algebra and Function Minimisation, 2nd Edition edition. Adam Hilger, Bristol.
Dennis JE, Gay DM, Walsh RE (1981). “An Adaptive Nonlinear Least-Squares Algorithm.” ACM Transactions on Mathematical Software, 7(3), 348–368. ISSN 00983500, doi:10.1145/355958.355965 , https://dl.acm.org/doi/10.1145/355958.355965.
See also
summary.maxlogL
, optim
, nlminb
,
DEoptim
, DEoptim.control
,
maxlogLreg
, bootstrap_maxlogL
Other maxlogL:
cum_hazard.maxlogL()
,
expected_value.maxlogL()
,
maxlogLreg()
Author
Jaime Mosquera Gutiérrez, jmosquerag@unal.edu.co
Examples
library(EstimationTools)
#----------------------------------------------------------------------------
# Example 1: estimation with one fixed parameter
x <- rnorm(n = 10000, mean = 160, sd = 6)
theta_1 <- maxlogL(x = x, dist = 'dnorm', control = list(trace = 1),
link = list(over = "sd", fun = "log_link"),
fixed = list(mean = 160))
#> 0: 44191.390: 1.00000
#> 1: 32573.039: 2.00000
#> 2: 32371.180: 1.92526
#> 3: 32244.045: 1.77840
#> 4: 32237.083: 1.80803
#> 5: 32236.976: 1.80485
#> 6: 32236.976: 1.80476
#> 7: 32236.976: 1.80476
summary(theta_1)
#> _______________________________________________________________
#> Optimization routine: nlminb
#> Standard Error calculation: Hessian from optim
#> _______________________________________________________________
#> AIC BIC
#> 64473.95 64473.95
#> _______________________________________________________________
#> Estimate Std. Error Z value Pr(>|z|)
#> sd 6.07851 0.04298 141.4 <2e-16 ***
#> ---
#> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#> _______________________________________________________________
#> Note: p-values valid under asymptotic normality of estimators
#> ---
#----------------------------------------------------------------------------
# Example 2: both parameters of normal distribution mapped with logarithmic
# function
theta_2 <- maxlogL(x = x, dist = "dnorm",
link = list(over = c("mean","sd"),
fun = c("log_link","log_link")))
summary(theta_2)
#> _______________________________________________________________
#> Optimization routine: nlminb
#> Standard Error calculation: Hessian from optim
#> _______________________________________________________________
#> AIC BIC
#> 64477.59 64492.01
#> _______________________________________________________________
#> Estimate Std. Error Z value Pr(>|z|)
#> mean 160.03662 0.06078 2632.9 <2e-16 ***
#> sd 6.07840 0.04298 141.4 <2e-16 ***
#> ---
#> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#> _______________________________________________________________
#> Note: p-values valid under asymptotic normality of estimators
#> ---
#--------------------------------------------------------------------------------
# Example 3: parameter estimation in ZIP distribution
if (!require('gamlss.dist')) install.packages('gamlss.dist')
#> Loading required package: gamlss.dist
library(gamlss.dist)
z <- rZIP(n=1000, mu=6, sigma=0.08)
theta_3 <- maxlogL(x = z, dist = 'dZIP', start = c(0, 0),
lower = c(-Inf, -Inf), upper = c(Inf, Inf),
optimizer = 'optim',
link = list(over=c("mu", "sigma"),
fun = c("log_link", "logit_link")))
summary(theta_3)
#> _______________________________________________________________
#> Optimization routine: optim
#> Standard Error calculation: Hessian from optim
#> _______________________________________________________________
#> AIC BIC
#> 4832.682 4842.498
#> _______________________________________________________________
#> Estimate Std. Error Z value Pr(>|z|)
#> mu 6.01217 0.08139 73.872 <2e-16 ***
#> sigma 0.07869 0.00865 9.096 <2e-16 ***
#> ---
#> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#> _______________________________________________________________
#> Note: p-values valid under asymptotic normality of estimators
#> ---
#--------------------------------------------------------------------------------
# Example 4: parameter estimation with fixed noncentrality parameter.
y_2 <- rbeta(n = 1000, shape1 = 2, shape2 = 3)
theta_41 <- maxlogL(x = y_2, dist = "dbeta",
link = list(over = c("shape1", "shape2"),
fun = c("log_link","log_link")))
summary(theta_41)
#> _______________________________________________________________
#> Optimization routine: nlminb
#> Standard Error calculation: Hessian from optim
#> _______________________________________________________________
#> AIC BIC
#> -455.0972 -450.1894
#> _______________________________________________________________
#> Estimate Std. Error Z value Pr(>|z|)
#> shape1 1.90876 0.08001 23.86 <2e-16 ***
#> shape2 2.93345 0.12854 22.82 <2e-16 ***
#> ---
#> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#> _______________________________________________________________
#> Note: p-values valid under asymptotic normality of estimators
#> ---
# It is also possible define 'ncp' as fixed parameter
theta_42 <- maxlogL(x = y_2, dist = "dbeta", fixed = list(ncp = 0),
link = list(over = c("shape1", "shape2"),
fun = c("log_link","log_link")) )
summary(theta_42)
#> _______________________________________________________________
#> Optimization routine: nlminb
#> Standard Error calculation: Hessian from optim
#> _______________________________________________________________
#> AIC BIC
#> -455.0972 -450.1894
#> _______________________________________________________________
#> Estimate Std. Error Z value Pr(>|z|)
#> shape1 1.90876 0.08001 23.86 <2e-16 ***
#> shape2 2.93345 0.12854 22.82 <2e-16 ***
#> ---
#> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#> _______________________________________________________________
#> Note: p-values valid under asymptotic normality of estimators
#> ---
#--------------------------------------------------------------------------------