The R Project SVN
R
(root)/branches/R-4-0-branch/tests/Examples/stats4-Ex.Rout.save – Rev 79826
Rev 79729 | Blame | Compare with Previous | Last modification | View Log | Download | RSS feed
R version 4.0.3 Patched (2021-01-08 r79820) -- "Bunny-Wunnies Freak Out"Copyright (C) 2021 The R Foundation for Statistical ComputingPlatform: 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 localeR 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 xhalf24.9931 3.0571> mle(ll, fixed = list(xhalf = 6))Call:mle(minuslogl = ll, fixed = list(xhalf = 6))Coefficients:ymax xhalf19.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 xhalf24.9994 3.0558>> AIC(fit)[1] 61.208> BIC(fit)[1] 62.004>> summary(fit)Maximum likelihood estimationCall:mle(minuslogl = ll, nobs = length(y))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(ll2, lower = c(0, 0)))Call:mle(minuslogl = ll2, lower = c(0, 0))Coefficients:ymax xhalf24.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 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(ll, 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(ll, lower = c(0, 4), fixed=list(ymax=23)) # has max on boundaryCall:mle(minuslogl = ll, fixed = list(ymax = 23), lower = c(0, 4))Coefficients:ymax xhalf23 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 MLECall:lm(formula = y ~ x)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, lower=c(-Inf,-Inf, 0.01)))Maximum likelihood estimationCall:mle(minuslogl = mll, 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, lower=list(sigma = 0.01))) # alternative specificationMaximum likelihood estimationCall:mle(minuslogl = mll, 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, 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, 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: 2.238 0.022 2.331 0 0> grDevices::dev.off()null device1> ###> ### Local variables: ***> ### mode: outline-minor ***> ### outline-regexp: "\\(> \\)?### [*]+" ***> ### End: ***> quit('no')