| |
| R version 3.6.2 Patched (2020-02-12 r77795) -- "Dark and Stormy Night" |
| Copyright (C) 2020 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) |
| > 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)) |
| > # For 1D, this is preferable: |
| > fit1 <- mle(nLL, start = list(lambda = 5), nobs = NROW(y), |
| + method = "Brent", lower = 1, upper = 20) |
| > stopifnot(nobs(fit0) == length(y)) |
| > |
| > ## 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, method = "L-BFGS-B", lower = rep(0, 2)) |
| |
| Call: |
| mle(minuslogl = ll2, method = "L-BFGS-B", 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(ll, method = "L-BFGS-B", lower = c(0, 0))) |
| |
| Call: |
| mle(minuslogl = ll, method = "L-BFGS-B", 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 |
| > |
| > 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.567 0.006 0.577 0 0 |
| > grDevices::dev.off() |
| null device |
| 1 |
| > ### |
| > ### Local variables: *** |
| > ### mode: outline-minor *** |
| > ### outline-regexp: "\\(> \\)?### [*]+" *** |
| > ### End: *** |
| > quit('no') |
