The R Project SVN R

R Under development (unstable) (2025-09-03 r88777) -- "Unsuffered Consequences"

Copyright (C) 2025 The R Foundation for Statistical Computing

Platform: x86_64-apple-darwin21.6.0

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.

> pkgname <- "stats4"

> source(file.path(R.home("share"), "R", "examples-header.R"))

> options(warn = 1)

> library('stats4')

> 

> base::assign(".oldSearch", base::search(), pos = 'CheckExEnv')

> base::assign(".old_wd", base::getwd(), pos = 'CheckExEnv')

> cleanEx()

> nameEx("mle")

> ### * mle

> 

> flush(stderr()); flush(stdout())

> 

> ### Name: mle

> ### Title: Maximum Likelihood Estimation

> ### Aliases: mle

> ### Keywords: models

> 

> ### ** Examples

> 

> ## Avoid printing to unwarranted accuracy

> od <- options(digits = 5)

> 

> ## Simulated EC50 experiment with count data

> ec50 <- data.frame(

+ x = 0:10,

+ y = c(26, 17, 13, 12, 20, 5, 9, 8, 5, 4, 8))

> 

> ## Easy one-dimensional MLE:

> nLL <- function(lambda) with(ec50,-sum(stats::dpois(y, lambda, log = TRUE)))

> fit0 <- mle(nLL, start = list(lambda = 5), nobs = NROW(ec50))

> 

> ## sanity check --- notice that "nobs" must be input

> ## (not guaranteed to be meaningful for any likelihood)

> stopifnot(nobs(fit0) == length(ec50$y))

> 

> 

> # For 1D, this is preferable:

> fit1 <- mle(nLL, start = list(lambda = 5), nobs = NROW(ec50),

+             method = "Brent", lower = 1, upper = 20)

> 

> ## This needs a constrained parameter space: most methods will accept NA

> mll1 <- function(ymax = 15, xhalf = 6) {

+     if(ymax > 0 && xhalf > 0) 

+       with(ec50, -sum(stats::dpois(y, lambda = ymax/(1+x/xhalf), log = TRUE)))

+     else NA

+ }

> (fit <- mle(mll1, nobs = NROW(ec50)))

Call:

mle(minuslogl = mll1, nobs = NROW(ec50))

Coefficients:

   ymax   xhalf 

24.9931  3.0571 

> mle(mll1, fixed = list(xhalf = 6))

Call:

mle(minuslogl = mll1, fixed = list(xhalf = 6))

Coefficients:

  ymax  xhalf 

19.288  6.000 

> 

> ## Alternative using bounds on optimization

> mll2 <- function(ymax = 15, xhalf = 6)

+     with(ec50, -sum(stats::dpois(y, lambda = ymax/(1+x/xhalf), log = TRUE)))

> mle(mll2, lower = rep(0, 2))

Call:

mle(minuslogl = mll2, lower = rep(0, 2))

Coefficients:

   ymax   xhalf 

24.9994  3.0558 

> 

> AIC(fit)

[1] 61.208

> BIC(fit)

[1] 62.004

> 

> summary(fit)

Maximum likelihood estimation

Call:

mle(minuslogl = mll1, nobs = NROW(ec50))

Coefficients:

      Estimate Std. Error

ymax   24.9931     4.2244

xhalf   3.0571     1.0348

-2 log L: 57.208 

> logLik(fit)

'log Lik.' -28.604 (df=2)

> vcov(fit)

         ymax   xhalf

ymax  17.8459 -3.7206

xhalf -3.7206  1.0708

> plot(profile(fit), absVal = FALSE)

> confint(fit)

Profiling...

        2.5 %  97.5 %

ymax  17.8845 34.6194

xhalf  1.6616  6.4792

> 

> ## Use bounded optimization

> ## The lower bounds are really > 0,

> ## but we use >=0 to stress-test profiling

> (fit2 <- mle(mll2, lower = c(0, 0)))

Call:

mle(minuslogl = mll2, lower = c(0, 0))

Coefficients:

   ymax   xhalf 

24.9994  3.0558 

> plot(profile(fit2), absVal = FALSE)

> 

> ## A better parametrization:

> mll3 <- function(lymax = log(15), lxhalf = log(6))

+     with(ec50, -sum(stats::dpois(y, lambda = exp(lymax)/(1+x/exp(lxhalf)), log = TRUE)))

> (fit3 <- mle(mll3))

Call:

mle(minuslogl = mll3)

Coefficients:

 lymax lxhalf 

3.2189 1.1170 

> plot(profile(fit3), absVal = FALSE)

> exp(confint(fit3))

Profiling...

         2.5 %  97.5 %

lymax  17.8815 34.6186

lxhalf  1.6615  6.4794

> 

> # Regression tests for bounded cases (this was broken in R 3.x)

> fit4 <- mle(mll1, lower = c(0, 4)) # has max on boundary

> confint(fit4)

Profiling...

       2.5 %  97.5 %

ymax  17.446 26.5081

xhalf     NA  6.9109

> 

> ## direct check that fixed= and constraints work together

> mle(mll1, lower = c(0, 4), fixed=list(ymax=23)) # has max on boundary

Call:

mle(minuslogl = mll1, fixed = list(ymax = 23), lower = c(0, 4))

Coefficients:

 ymax xhalf 

   23     4 

> 

> ## Linear regression using MLE

> lr <- data.frame(

+ x = 1:10, 

+ y = c(0.48, 2.24, 2.22, 5.15, 4.64, 5.53, 7, 8.8, 7.67, 9.23))

> 

> 

> LM_mll <- function(formula, data = environment(formula))

+ {

+      y <- model.response(model.frame(formula, data))

+      X <- model.matrix(formula, data)

+      b0 <- numeric(NCOL(X))

+      names(b0) <- colnames(X)

+      function(b=b0, sigma=1)

+          -sum(dnorm(y, X %*% b, sigma, log=TRUE))

+ }

> 

> mll_lm <- LM_mll(y ~ x, lr)

> 

> summary(lm(y~x, data=lr)) # for comparison -- notice variance bias in MLE

Call:

lm(formula = y ~ x, data = lr)

Residuals:

   Min     1Q Median     3Q    Max 

-0.937 -0.500 -0.211  0.278  1.273 

Coefficients:

            Estimate Std. Error t value Pr(>|t|)    

(Intercept)   0.0927     0.5376    0.17     0.87    

x             0.9461     0.0866   10.92  4.4e-06 ***

---

Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.787 on 8 degrees of freedom

Multiple R-squared:  0.937,	Adjusted R-squared:  0.929 

F-statistic:  119 on 1 and 8 DF,  p-value: 4.39e-06

> summary(mle(mll_lm, lower=c(-Inf,-Inf, 0.01)))

Maximum likelihood estimation

Call:

mle(minuslogl = mll_lm, lower = c(-Inf, -Inf, 0.01))

Coefficients:

              Estimate Std. Error

b.(Intercept) 0.092667   0.480869

b.x           0.946061   0.077499

sigma         0.703919   0.157400

-2 log L: 21.357 

> summary(mle(mll_lm, lower=list(sigma = 0.01))) # alternative specification

Maximum likelihood estimation

Call:

mle(minuslogl = mll_lm, lower = list(sigma = 0.01))

Coefficients:

              Estimate Std. Error

b.(Intercept) 0.092667   0.480869

b.x           0.946061   0.077499

sigma         0.703919   0.157400

-2 log L: 21.357 

> 

> confint(mle(mll_lm, lower=list(sigma = 0.01)))

Profiling...

                 2.5 % 97.5 %

b.(Intercept) -0.94831 1.1336

b.x            0.77829 1.1138

sigma          0.48017 1.1755

> plot(profile(mle(mll_lm, lower=list(sigma = 0.01))))

> 

> Binom_mll <- function(x, n)

+ {

+     force(x); force(n) ## beware lazy evaluation

+     function(p=.5) -dbinom(x, n, p, log=TRUE)

+ }

> 

> ## Likelihood functions for different x.

> ## This code goes wrong, if force(x) is not used in Binom_mll:

> 

> curve(Binom_mll(0, 10)(p), xname="p", ylim=c(0, 10))

> mll_list <- list(10)

> for (x in 1:10)

+     mll_list[[x]] <- Binom_mll(x, 10)

> for (mll in mll_list)

+     curve(mll(p), xname="p", add=TRUE)

> 

> mll <- Binom_mll(4,10)

> mle(mll, lower = 1e-16, upper = 1-1e-16) # limits must be inside (0,1)

Call:

mle(minuslogl = mll, lower = 1e-16, upper = 1 - 1e-16)

Coefficients:

  p 

0.4 

> 

> ## Boundary case: This works, but fails if limits are set closer to 0 and 1  

> mll <- Binom_mll(0, 10)

> mle(mll, lower=.005, upper=.995)

Call:

mle(minuslogl = mll, lower = 0.005, upper = 0.995)

Coefficients:

    p 

0.005 

> 

> ## Not run: 

> ##D ## We can use limits closer to the boundaries if we use the

> ##D ## drop-in replacement optimr() from the optimx package.

> ##D 

> ##D mle(mll, lower = 1e-16, upper = 1-1e-16, optim=optimx::optimr)

> ## End(Not run)

> 

> 

> options(od)

> 

> 

> 

> cleanEx()

> nameEx("update-methods")

> ### * update-methods

> 

> flush(stderr()); flush(stdout())

> 

> ### Name: update-methods

> ### Title: Methods for Function 'update' in Package 'stats4'

> ### Aliases: update-methods update,ANY-method update,mle-method

> ### Keywords: methods

> 

> ### ** Examples

> 

> x <- 0:10

> y <- c(26, 17, 13, 12, 20, 5, 9, 8, 5, 4, 8)

> ll <- function(ymax = 15, xhalf = 6)

+     -sum(stats::dpois(y, lambda = ymax/(1+x/xhalf), log = TRUE))

> fit <- mle(ll)

Warning in stats::dpois(y, lambda = ymax/(1 + x/xhalf), log = TRUE) :

  NaNs produced

> ## note the recorded call contains ..1, a problem with S4 dispatch

> update(fit, fixed = list(xhalf = 3))

Call:

mle(minuslogl = ll, fixed = ..1)

Coefficients:

    ymax    xhalf 

25.19609  3.00000 

> 

> 

> 

> ### * <FOOTER>

> ###

> cleanEx()

> options(digits = 7L)

> base::cat("Time elapsed: ", proc.time() - base::get("ptime", pos = 'CheckExEnv'),"\n")

Time elapsed:  0.574 0.015 0.593 0 0 

> grDevices::dev.off()

null device 

          1 

> ###

> ### Local variables: ***

> ### mode: outline-minor ***

> ### outline-regexp: "\\(> \\)?### [*]+" ***

> ### End: ***

> quit('no')