R version 2.13.0 Under development (unstable) (2011-02-26 r54603) Copyright (C) 2011 The R Foundation for Statistical Computing ISBN 3-900051-07-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. 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) > ## 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 > (fit1 <- 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(fit1), absVal=FALSE) > > ## a better parametrization: > ll2 <- function(lymax = log(15), lxhalf = log(6)) + -sum(stats::dpois(y, lambda=exp(lymax)/(1+x/exp(lxhalf)), log=TRUE)) > (fit2 <- mle(ll2)) Call: mle(minuslogl = ll2) Coefficients: lymax lxhalf 3.2189 1.1170 > plot(profile(fit2), absVal = FALSE) > exp(confint(fit2)) 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 > > > > ### *