The R Project SVN
R
(root)/trunk/tests/Examples/stats4-Ex.Rout.save @ 90215 – Rev 88779
Rev 88768 | Blame | Compare with Previous | Last modification | View Log | Download | RSS feed
R Under development (unstable) (2025-09-03 r88777) -- "Unsuffered Consequences"Copyright (C) 2025 The R Foundation for Statistical ComputingPlatform: x86_64-apple-darwin21.6.0R 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 xhalf24.9931 3.0571> mle(mll1, fixed = list(xhalf = 6))Call:mle(minuslogl = mll1, fixed = list(xhalf = 6))Coefficients:ymax xhalf19.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 xhalf24.9994 3.0558>> AIC(fit)[1] 61.208> BIC(fit)[1] 62.004>> summary(fit)Maximum likelihood estimationCall:mle(minuslogl = mll1, nobs = NROW(ec50))Coefficients:Estimate Std. Errorymax 24.9931 4.2244xhalf 3.0571 1.0348-2 log L: 57.208> logLik(fit)'log Lik.' -28.604 (df=2)> vcov(fit)ymax xhalfymax 17.8459 -3.7206xhalf -3.7206 1.0708> plot(profile(fit), absVal = FALSE)> confint(fit)Profiling...2.5 % 97.5 %ymax 17.8845 34.6194xhalf 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 xhalf24.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 lxhalf3.2189 1.1170> plot(profile(fit3), absVal = FALSE)> exp(confint(fit3))Profiling...2.5 % 97.5 %lymax 17.8815 34.6186lxhalf 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.5081xhalf NA 6.9109>> ## direct check that fixed= and constraints work together> mle(mll1, lower = c(0, 4), fixed=list(ymax=23)) # has max on boundaryCall:mle(minuslogl = mll1, fixed = list(ymax = 23), lower = c(0, 4))Coefficients:ymax xhalf23 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 MLECall:lm(formula = y ~ x, data = lr)Residuals:Min 1Q Median 3Q Max-0.937 -0.500 -0.211 0.278 1.273Coefficients:Estimate Std. Error t value Pr(>|t|)(Intercept) 0.0927 0.5376 0.17 0.87x 0.9461 0.0866 10.92 4.4e-06 ***---Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1Residual standard error: 0.787 on 8 degrees of freedomMultiple R-squared: 0.937, Adjusted R-squared: 0.929F-statistic: 119 on 1 and 8 DF, p-value: 4.39e-06> summary(mle(mll_lm, lower=c(-Inf,-Inf, 0.01)))Maximum likelihood estimationCall:mle(minuslogl = mll_lm, lower = c(-Inf, -Inf, 0.01))Coefficients:Estimate Std. Errorb.(Intercept) 0.092667 0.480869b.x 0.946061 0.077499sigma 0.703919 0.157400-2 log L: 21.357> summary(mle(mll_lm, lower=list(sigma = 0.01))) # alternative specificationMaximum likelihood estimationCall:mle(minuslogl = mll_lm, lower = list(sigma = 0.01))Coefficients:Estimate Std. Errorb.(Intercept) 0.092667 0.480869b.x 0.946061 0.077499sigma 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.1336b.x 0.77829 1.1138sigma 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:p0.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:p0.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 xhalf25.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 device1> ###> ### Local variables: ***> ### mode: outline-minor ***> ### outline-regexp: "\\(> \\)?### [*]+" ***> ### End: ***> quit('no')