| |
| R Under development (unstable) (2022-03-19 r81942) -- "Unsuffered Consequences" |
| Copyright (C) 2022 The R Foundation for Statistical Computing |
| Platform: x86_64-pc-linux-gnu (64-bit) |
| |
| 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. |
| |
| Natural language support but running in an English locale |
| |
| 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 |
| > x <- 0:10 |
| > y <- c(26, 17, 13, 12, 20, 5, 9, 8, 5, 4, 8) |
| > |
| > ## Easy one-dimensional MLE: |
| > nLL <- function(lambda) -sum(stats::dpois(y, lambda, log = TRUE)) |
| > fit0 <- mle(nLL, start = list(lambda = 5), nobs = NROW(y)) |
| > |
| > ## sanity check --- notice that "nobs" must be input |
| > ## (not guaranteed to be meaningful for any likelihood) |
| > stopifnot(nobs(fit0) == length(y)) |
| > |
| > |
| > # For 1D, this is preferable: |
| > fit1 <- mle(nLL, start = list(lambda = 5), nobs = NROW(y), |
| + method = "Brent", lower = 1, upper = 20) |
| > |
| > ## This needs a constrained parameter space: most methods will accept NA |
| > ll <- function(ymax = 15, xhalf = 6) { |
| + if(ymax > 0 && xhalf > 0) |
| + -sum(stats::dpois(y, lambda = ymax/(1+x/xhalf), log = TRUE)) |
| + else NA |
| + } |
| > (fit <- mle(ll, nobs = length(y))) |
| |
| Call: |
| mle(minuslogl = ll, nobs = length(y)) |
| |
| Coefficients: |
| ymax xhalf |
| 24.9931 3.0571 |
| > mle(ll, fixed = list(xhalf = 6)) |
| |
| Call: |
| mle(minuslogl = ll, fixed = list(xhalf = 6)) |
| |
| Coefficients: |
| ymax xhalf |
| 19.288 6.000 |
| > |
| > ## Alternative using bounds on optimization |
| > ll2 <- function(ymax = 15, xhalf = 6) |
| + -sum(stats::dpois(y, lambda = ymax/(1+x/xhalf), log = TRUE)) |
| > mle(ll2, lower = rep(0, 2)) |
| |
| Call: |
| mle(minuslogl = ll2, 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 = ll, nobs = length(y)) |
| |
| 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(ll2, lower = c(0, 0))) |
| |
| Call: |
| mle(minuslogl = ll2, lower = c(0, 0)) |
| |
| Coefficients: |
| ymax xhalf |
| 24.9994 3.0558 |
| > plot(profile(fit2), absVal = FALSE) |
| > |
| > ## A better parametrization: |
| > ll3 <- function(lymax = log(15), lxhalf = log(6)) |
| + -sum(stats::dpois(y, lambda = exp(lymax)/(1+x/exp(lxhalf)), log = TRUE)) |
| > (fit3 <- mle(ll3)) |
| |
| Call: |
| mle(minuslogl = ll3) |
| |
| 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(ll, 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(ll, lower = c(0, 4), fixed=list(ymax=23)) # has max on boundary |
| |
| Call: |
| mle(minuslogl = ll, fixed = list(ymax = 23), lower = c(0, 4)) |
| |
| Coefficients: |
| ymax xhalf |
| 23 4 |
| > |
| > ## Linear regression using MLE |
| > 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_mll(y ~ x) |
| > |
| > summary(lm(y~x)) # for comparison -- notice variance bias in MLE |
| |
| Call: |
| lm(formula = y ~ x) |
| |
| 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, lower=c(-Inf,-Inf, 0.01))) |
| Maximum likelihood estimation |
| |
| Call: |
| mle(minuslogl = mll, 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, lower=list(sigma = 0.01))) # alternative specification |
| Maximum likelihood estimation |
| |
| Call: |
| mle(minuslogl = mll, 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, 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, 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: 2.16 0.024 2.2 0 0 |
| > grDevices::dev.off() |
| null device |
| 1 |
| > ### |
| > ### Local variables: *** |
| > ### mode: outline-minor *** |
| > ### outline-regexp: "\\(> \\)?### [*]+" *** |
| > ### End: *** |
| > quit('no') |