The R Project SVN R

R version 2.14.0 alpha (2011-10-03 r57134)
Copyright (C) 2011 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
Platform: x86_64-unknown-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')
> 
> assign(".oldSearch", search(), 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 dpois(x, lambda, log) : 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>
> ###
> cat("Time elapsed: ", proc.time() - get("ptime", pos = 'CheckExEnv'),"\n")
Time elapsed:  0.763 0.011 0.779 0 0 
> grDevices::dev.off()
null device 
          1 
> ###
> ### Local variables: ***
> ### mode: outline-minor ***
> ### outline-regexp: "\\(> \\)?### [*]+" ***
> ### End: ***
> quit('no')