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 > > > > ### *