Skip to contents

Maximum Likelihood Estimation for parametric distributions

Source: R/maxlogL.R
maxlogL.Rd

[Maturing]

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 and NegInv_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 and gaare available; custom optimization routines can also be implemented. nlminb is the default routine.

control

control parameters of the optimization routine. Please, visit documentation of selected optimizer for further information.

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, visit optim or hessian.

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:     44192.989:  1.00000
#>   1:     32573.255:  2.00000
#>   2:     32371.463:  1.92527
#>   3:     32244.360:  1.77845
#>   4:     32237.402:  1.80806
#>   5:     32237.296:  1.80488
#>   6:     32237.296:  1.80479
#>   7:     32237.296:  1.80479
summary(theta_1)
#> _______________________________________________________________
#> Optimization routine: nlminb 
#> Standard Error calculation: Hessian from optim 
#> _______________________________________________________________
#>        AIC      BIC
#>   64474.59 64474.59
#> _______________________________________________________________
#>    Estimate  Std. Error Z value Pr(>|z|)    
#> sd   6.07870    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
#>   64478.3 64492.72
#> _______________________________________________________________
#>      Estimate  Std. Error Z value Pr(>|z|)    
#> mean 160.03260    0.06079  2632.7   <2e-16 ***
#> sd     6.07861    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
#>   4831.278 4841.093
#> _______________________________________________________________
#>       Estimate  Std. Error Z value Pr(>|z|)    
#> mu     6.040718   0.081435   74.18   <2e-16 ***
#> sigma  0.075754   0.008502    8.91   <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
#>   -462.4313 -457.5235
#> _______________________________________________________________
#>        Estimate  Std. Error Z value Pr(>|z|)    
#> shape1   1.90433    0.07979   23.87   <2e-16 ***
#> shape2   2.96006    0.12979   22.81   <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
#>   -462.4313 -457.5235
#> _______________________________________________________________
#>        Estimate  Std. Error Z value Pr(>|z|)    
#> shape1   1.90433    0.07979   23.87   <2e-16 ***
#> shape2   2.96006    0.12979   22.81   <2e-16 ***
#> ---
#> Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#> _______________________________________________________________
#> Note: p-values valid under asymptotic normality of estimators 
#> ---


#--------------------------------------------------------------------------------