tests/Examples/stats4-Ex.Rout.save

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')