The R Project SVN R

Rev

Rev 77335 | Rev 77514 | Go to most recent revision | Details | Compare with Previous | Last modification | View Log | RSS feed

Rev Author Line No. Line
56186 murdoch 1 
 
77431 ripley 2 
R Under development (unstable) (2019-11-17 r77430) -- "Unsuffered Consequences"
75937 ripley 3 
Copyright (C) 2019 The R Foundation for Statistical Computing
68809 ripley 4 
Platform: x86_64-pc-linux-gnu (64-bit)
56186 murdoch 5 
 
6 
R is free software and comes with ABSOLUTELY NO WARRANTY.
7 
You are welcome to redistribute it under certain conditions.
8 
Type 'license()' or 'licence()' for distribution details.
9 
 
10 
  Natural language support but running in an English locale
11 
 
12 
R is a collaborative project with many contributors.
13 
Type 'contributors()' for more information and
14 
'citation()' on how to cite R or R packages in publications.
15 
 
16 
Type 'demo()' for some demos, 'help()' for on-line help, or
17 
'help.start()' for an HTML browser interface to help.
18 
Type 'q()' to quit R.
19 
 
20 
> pkgname <- "stats4"
21 
> source(file.path(R.home("share"), "R", "examples-header.R"))
22 
> options(warn = 1)
23 
> library('stats4')
24 
>
61787 ripley 25 
> base::assign(".oldSearch", base::search(), pos = 'CheckExEnv')
73819 hornik 26 
> base::assign(".old_wd", base::getwd(), pos = 'CheckExEnv')
56186 murdoch 27 
> cleanEx()
28 
> nameEx("mle")
29 
> ### * mle
30 
>
31 
> flush(stderr()); flush(stdout())
32 
>
33 
> ### Name: mle
34 
> ### Title: Maximum Likelihood Estimation
35 
> ### Aliases: mle
36 
> ### Keywords: models
37 
>
38 
> ### ** Examples
39 
>
40 
> ## Avoid printing to unwarranted accuracy
41 
> od <- options(digits = 5)
42 
> x <- 0:10
43 
> y <- c(26, 17, 13, 12, 20, 5, 9, 8, 5, 4, 8)
44 
>
45 
> ## Easy one-dimensional MLE:
61164 ripley 46 
> nLL <- function(lambda) -sum(stats::dpois(y, lambda, log = TRUE))
56186 murdoch 47 
> fit0 <- mle(nLL, start = list(lambda = 5), nobs = NROW(y))
48 
> # For 1D, this is preferable:
49 
> fit1 <- mle(nLL, start = list(lambda = 5), nobs = NROW(y),
50 
+             method = "Brent", lower = 1, upper = 20)
51 
> stopifnot(nobs(fit0) == length(y))
52 
>
53 
> ## This needs a constrained parameter space: most methods will accept NA
56281 ripley 54 
> ll <- function(ymax = 15, xhalf = 6) {
55 
+     if(ymax > 0 && xhalf > 0)
56 
+       -sum(stats::dpois(y, lambda = ymax/(1+x/xhalf), log = TRUE))
57 
+     else NA
58 
+ }
56186 murdoch 59 
> (fit <- mle(ll, nobs = length(y)))
60 
 
61 
Call:
62 
mle(minuslogl = ll, nobs = length(y))
63 
 
64 
Coefficients:
65 
   ymax   xhalf
66 
24.9931  3.0571
67 
> mle(ll, fixed = list(xhalf = 6))
68 
 
69 
Call:
70 
mle(minuslogl = ll, fixed = list(xhalf = 6))
71 
 
72 
Coefficients:
73 
  ymax  xhalf
74 
19.288  6.000
75 
> ## alternative using bounds on optimization
76 
> ll2 <- function(ymax = 15, xhalf = 6)
61164 ripley 77 
+     -sum(stats::dpois(y, lambda = ymax/(1+x/xhalf), log = TRUE))
56186 murdoch 78 
> mle(ll2, method = "L-BFGS-B", lower = rep(0, 2))
79 
 
80 
Call:
81 
mle(minuslogl = ll2, method = "L-BFGS-B", lower = rep(0, 2))
82 
 
83 
Coefficients:
84 
   ymax   xhalf
85 
24.9994  3.0558
86 
>
87 
> AIC(fit)
88 
[1] 61.208
89 
> BIC(fit)
90 
[1] 62.004
91 
>
92 
> summary(fit)
93 
Maximum likelihood estimation
94 
 
95 
Call:
96 
mle(minuslogl = ll, nobs = length(y))
97 
 
98 
Coefficients:
99 
      Estimate Std. Error
100 
ymax   24.9931     4.2244
101 
xhalf   3.0571     1.0348
102 
 
103 
-2 log L: 57.208
104 
> logLik(fit)
105 
'log Lik.' -28.604 (df=2)
106 
> vcov(fit)
107 
         ymax   xhalf
108 
ymax  17.8459 -3.7206
109 
xhalf -3.7206  1.0708
110 
> plot(profile(fit), absVal = FALSE)
111 
> confint(fit)
112 
Profiling...
113 
        2.5 %  97.5 %
114 
ymax  17.8845 34.6194
115 
xhalf  1.6616  6.4792
116 
>
56281 ripley 117 
> ## Use bounded optimization
118 
> ## The lower bounds are really > 0,
119 
> ## but we use >=0 to stress-test profiling
56186 murdoch 120 
> (fit2 <- mle(ll, method = "L-BFGS-B", lower = c(0, 0)))
121 
 
122 
Call:
123 
mle(minuslogl = ll, method = "L-BFGS-B", lower = c(0, 0))
124 
 
125 
Coefficients:
126 
   ymax   xhalf
127 
24.9994  3.0558
61164 ripley 128 
> plot(profile(fit2), absVal = FALSE)
56186 murdoch 129 
>
130 
> ## a better parametrization:
131 
> ll3 <- function(lymax = log(15), lxhalf = log(6))
61164 ripley 132 
+     -sum(stats::dpois(y, lambda = exp(lymax)/(1+x/exp(lxhalf)), log = TRUE))
56186 murdoch 133 
> (fit3 <- mle(ll3))
134 
 
135 
Call:
136 
mle(minuslogl = ll3)
137 
 
138 
Coefficients:
139 
 lymax lxhalf
140 
3.2189 1.1170
141 
> plot(profile(fit3), absVal = FALSE)
142 
> exp(confint(fit3))
143 
Profiling...
144 
         2.5 %  97.5 %
145 
lymax  17.8815 34.6186
146 
lxhalf  1.6615  6.4794
147 
>
148 
> options(od)
149 
>
150 
>
151 
>
152 
> cleanEx()
153 
> nameEx("update-methods")
154 
> ### * update-methods
155 
>
156 
> flush(stderr()); flush(stdout())
157 
>
158 
> ### Name: update-methods
159 
> ### Title: Methods for Function 'update' in Package 'stats4'
160 
> ### Aliases: update-methods update,ANY-method update,mle-method
161 
> ### Keywords: methods
162 
>
163 
> ### ** Examples
164 
>
165 
> x <- 0:10
166 
> y <- c(26, 17, 13, 12, 20, 5, 9, 8, 5, 4, 8)
61164 ripley 167 
> ll <- function(ymax = 15, xhalf = 6)
168 
+     -sum(stats::dpois(y, lambda = ymax/(1+x/xhalf), log = TRUE))
56186 murdoch 169 
> fit <- mle(ll)
59447 ripley 170 
Warning in stats::dpois(y, lambda = ymax/(1 + x/xhalf), log = TRUE) :
171 
  NaNs produced
56186 murdoch 172 
> ## note the recorded call contains ..1, a problem with S4 dispatch
61164 ripley 173 
> update(fit, fixed = list(xhalf = 3))
56186 murdoch 174 
 
175 
Call:
176 
mle(minuslogl = ll, fixed = ..1)
177 
 
178 
Coefficients:
179 
    ymax    xhalf
180 
25.19609  3.00000
181 
>
182 
>
183 
>
184 
> ### * <FOOTER>
185 
> ###
73819 hornik 186 
> cleanEx()
62439 ripley 187 
> options(digits = 7L)
61787 ripley 188 
> base::cat("Time elapsed: ", proc.time() - base::get("ptime", pos = 'CheckExEnv'),"\n")
77431 ripley 189 
Time elapsed:  0.566 0.02 0.589 0 0
56186 murdoch 190 
> grDevices::dev.off()
191 
null device
192 
          1
193 
> ###
194 
> ### Local variables: ***
195 
> ### mode: outline-minor ***
196 
> ### outline-regexp: "\\(> \\)?### [*]+" ***
197 
> ### End: ***
198 
> quit('no')