The R Project SVN R

Rev

Rev 88768 | Show entire file | Ignore whitespace | Details | Blame | Last modification | View Log | RSS feed

Rev 88768 Rev 88779
Line 1... Line 1...
1 
 
 1 
 
2 
R Under development (unstable) (2025-09-02 r88767) -- "Unsuffered Consequences"
 2 
R Under development (unstable) (2025-09-03 r88777) -- "Unsuffered Consequences"
3 
Copyright (C) 2025 The R Foundation for Statistical Computing
 3 
Copyright (C) 2025 The R Foundation for Statistical Computing
4 
Platform: x86_64-apple-darwin21.6.0
 4 
Platform: x86_64-apple-darwin21.6.0
5 
 
 5 
 
6 
R is free software and comes with ABSOLUTELY NO WARRANTY.
 6 
R is free software and comes with ABSOLUTELY NO WARRANTY.
7 
You are welcome to redistribute it under certain conditions.
 7 
You are welcome to redistribute it under certain conditions.
Line 55... Line 55...
55 
> # For 1D, this is preferable:
 55 
> # For 1D, this is preferable:
56 
> fit1 <- mle(nLL, start = list(lambda = 5), nobs = NROW(ec50),
 56 
> fit1 <- mle(nLL, start = list(lambda = 5), nobs = NROW(ec50),
57 
+             method = "Brent", lower = 1, upper = 20)
 57 
+             method = "Brent", lower = 1, upper = 20)
58 
>
58 
>
59 
> ## This needs a constrained parameter space: most methods will accept NA
 59 
> ## This needs a constrained parameter space: most methods will accept NA
60 
> ll <- function(ymax = 15, xhalf = 6) {
 60 
> mll1 <- function(ymax = 15, xhalf = 6) {
61 
+     if(ymax > 0 && xhalf > 0)
61 
+     if(ymax > 0 && xhalf > 0)
62 
+       with(ec50, -sum(stats::dpois(y, lambda = ymax/(1+x/xhalf), log = TRUE)))
 62 
+       with(ec50, -sum(stats::dpois(y, lambda = ymax/(1+x/xhalf), log = TRUE)))
63 
+     else NA
 63 
+     else NA
64 
+ }
 64 
+ }
65 
> (fit <- mle(ll, nobs = NROW(ec50)))
 65 
> (fit <- mle(mll1, nobs = NROW(ec50)))
66 
 
 66 
 
67 
Call:
 67 
Call:
68 
mle(minuslogl = ll, nobs = NROW(ec50))
 68 
mle(minuslogl = mll1, nobs = NROW(ec50))
69 
 
 69 
 
70 
Coefficients:
 70 
Coefficients:
71 
   ymax   xhalf
71 
   ymax   xhalf
72 
24.9931  3.0571
72 
24.9931  3.0571
73 
> mle(ll, fixed = list(xhalf = 6))
 73 
> mle(mll1, fixed = list(xhalf = 6))
74 
 
 74 
 
75 
Call:
 75 
Call:
76 
mle(minuslogl = ll, fixed = list(xhalf = 6))
 76 
mle(minuslogl = mll1, fixed = list(xhalf = 6))
77 
 
 77 
 
78 
Coefficients:
 78 
Coefficients:
79 
  ymax  xhalf
79 
  ymax  xhalf
80 
19.288  6.000
80 
19.288  6.000
81 
>
81 
>
82 
> ## Alternative using bounds on optimization
 82 
> ## Alternative using bounds on optimization
83 
> ll2 <- function(ymax = 15, xhalf = 6)
 83 
> mll2 <- function(ymax = 15, xhalf = 6)
84 
+     with(ec50, -sum(stats::dpois(y, lambda = ymax/(1+x/xhalf), log = TRUE)))
 84 
+     with(ec50, -sum(stats::dpois(y, lambda = ymax/(1+x/xhalf), log = TRUE)))
85 
> mle(ll2, lower = rep(0, 2))
 85 
> mle(mll2, lower = rep(0, 2))
86 
 
 86 
 
87 
Call:
 87 
Call:
88 
mle(minuslogl = ll2, lower = rep(0, 2))
 88 
mle(minuslogl = mll2, lower = rep(0, 2))
89 
 
 89 
 
90 
Coefficients:
 90 
Coefficients:
91 
   ymax   xhalf
91 
   ymax   xhalf
92 
24.9994  3.0558
92 
24.9994  3.0558
93 
>
93 
>
Line 98... Line 98...
98 
>
98 
>
99 
> summary(fit)
 99 
> summary(fit)
100 
Maximum likelihood estimation
 100 
Maximum likelihood estimation
101 
 
 101 
 
102 
Call:
 102 
Call:
103 
mle(minuslogl = ll, nobs = NROW(ec50))
 103 
mle(minuslogl = mll1, nobs = NROW(ec50))
104 
 
 104 
 
105 
Coefficients:
 105 
Coefficients:
106 
      Estimate Std. Error
 106 
      Estimate Std. Error
107 
ymax   24.9931     4.2244
 107 
ymax   24.9931     4.2244
108 
xhalf   3.0571     1.0348
 108 
xhalf   3.0571     1.0348
Line 122... Line 122...
122 
xhalf  1.6616  6.4792
 122 
xhalf  1.6616  6.4792
123 
>
123 
>
124 
> ## Use bounded optimization
 124 
> ## Use bounded optimization
125 
> ## The lower bounds are really > 0,
 125 
> ## The lower bounds are really > 0,
126 
> ## but we use >=0 to stress-test profiling
 126 
> ## but we use >=0 to stress-test profiling
127 
> (fit2 <- mle(ll2, lower = c(0, 0)))
 127 
> (fit2 <- mle(mll2, lower = c(0, 0)))
128 
 
 128 
 
129 
Call:
 129 
Call:
130 
mle(minuslogl = ll2, lower = c(0, 0))
 130 
mle(minuslogl = mll2, lower = c(0, 0))
131 
 
 131 
 
132 
Coefficients:
 132 
Coefficients:
133 
   ymax   xhalf
133 
   ymax   xhalf
134 
24.9994  3.0558
134 
24.9994  3.0558
135 
> plot(profile(fit2), absVal = FALSE)
 135 
> plot(profile(fit2), absVal = FALSE)
136 
>
136 
>
137 
> ## A better parametrization:
 137 
> ## A better parametrization:
138 
> ll3 <- function(lymax = log(15), lxhalf = log(6))
 138 
> mll3 <- function(lymax = log(15), lxhalf = log(6))
139 
+     with(ec50, -sum(stats::dpois(y, lambda = exp(lymax)/(1+x/exp(lxhalf)), log = TRUE)))
 139 
+     with(ec50, -sum(stats::dpois(y, lambda = exp(lymax)/(1+x/exp(lxhalf)), log = TRUE)))
140 
> (fit3 <- mle(ll3))
 140 
> (fit3 <- mle(mll3))
141 
 
 141 
 
142 
Call:
 142 
Call:
143 
mle(minuslogl = ll3)
 143 
mle(minuslogl = mll3)
144 
 
 144 
 
145 
Coefficients:
 145 
Coefficients:
146 
 lymax lxhalf
146 
 lymax lxhalf
147 
3.2189 1.1170
147 
3.2189 1.1170
148 
> plot(profile(fit3), absVal = FALSE)
 148 
> plot(profile(fit3), absVal = FALSE)
Line 151... Line 151...
151 
         2.5 %  97.5 %
 151 
         2.5 %  97.5 %
152 
lymax  17.8815 34.6186
 152 
lymax  17.8815 34.6186
153 
lxhalf  1.6615  6.4794
 153 
lxhalf  1.6615  6.4794
154 
>
154 
>
155 
> # Regression tests for bounded cases (this was broken in R 3.x)
 155 
> # Regression tests for bounded cases (this was broken in R 3.x)
156 
> fit4 <- mle(ll, lower = c(0, 4)) # has max on boundary
 156 
> fit4 <- mle(mll1, lower = c(0, 4)) # has max on boundary
157 
> confint(fit4)
 157 
> confint(fit4)
158 
Profiling...
 158 
Profiling...
159 
       2.5 %  97.5 %
 159 
       2.5 %  97.5 %
160 
ymax  17.446 26.5081
 160 
ymax  17.446 26.5081
161 
xhalf     NA  6.9109
 161 
xhalf     NA  6.9109
162 
>
162 
>
163 
> ## direct check that fixed= and constraints work together
 163 
> ## direct check that fixed= and constraints work together
164 
> mle(ll, lower = c(0, 4), fixed=list(ymax=23)) # has max on boundary
 164 
> mle(mll1, lower = c(0, 4), fixed=list(ymax=23)) # has max on boundary
165 
 
 165 
 
166 
Call:
 166 
Call:
167 
mle(minuslogl = ll, fixed = list(ymax = 23), lower = c(0, 4))
 167 
mle(minuslogl = mll1, fixed = list(ymax = 23), lower = c(0, 4))
168 
 
 168 
 
169 
Coefficients:
 169 
Coefficients:
170 
 ymax xhalf
170 
 ymax xhalf
171 
   23     4
171 
   23     4
172 
>
172 
>
Line 184... Line 184...
184 
+      names(b0) <- colnames(X)
 184 
+      names(b0) <- colnames(X)
185 
+      function(b=b0, sigma=1)
 185 
+      function(b=b0, sigma=1)
186 
+          -sum(dnorm(y, X %*% b, sigma, log=TRUE))
 186 
+          -sum(dnorm(y, X %*% b, sigma, log=TRUE))
187 
+ }
 187 
+ }
188 
>
188 
>
189 
> mll <- LM_mll(y ~ x, lr)
 189 
> mll_lm <- LM_mll(y ~ x, lr)
190 
>
190 
>
191 
> summary(lm(y~x, data=lr)) # for comparison -- notice variance bias in MLE
 191 
> summary(lm(y~x, data=lr)) # for comparison -- notice variance bias in MLE
192 
 
 192 
 
193 
Call:
 193 
Call:
194 
lm(formula = y ~ x, data = lr)
 194 
lm(formula = y ~ x, data = lr)
Line 206... Line 206...
206 
 
 206 
 
207 
Residual standard error: 0.787 on 8 degrees of freedom
 207 
Residual standard error: 0.787 on 8 degrees of freedom
208 
Multiple R-squared:  0.937,	Adjusted R-squared:  0.929
208 
Multiple R-squared:  0.937,	Adjusted R-squared:  0.929
209 
F-statistic:  119 on 1 and 8 DF,  p-value: 4.39e-06
 209 
F-statistic:  119 on 1 and 8 DF,  p-value: 4.39e-06
210 
 
 210 
 
211 
> summary(mle(mll, lower=c(-Inf,-Inf, 0.01)))
 211 
> summary(mle(mll_lm, lower=c(-Inf,-Inf, 0.01)))
212 
Maximum likelihood estimation
 212 
Maximum likelihood estimation
213 
 
 213 
 
214 
Call:
 214 
Call:
215 
mle(minuslogl = mll, lower = c(-Inf, -Inf, 0.01))
 215 
mle(minuslogl = mll_lm, lower = c(-Inf, -Inf, 0.01))
216 
 
 216 
 
217 
Coefficients:
 217 
Coefficients:
218 
              Estimate Std. Error
 218 
              Estimate Std. Error
219 
b.(Intercept) 0.092667   0.480869
 219 
b.(Intercept) 0.092667   0.480869
220 
b.x           0.946061   0.077499
 220 
b.x           0.946061   0.077499
221 
sigma         0.703919   0.157400
 221 
sigma         0.703919   0.157400
222 
 
 222 
 
223 
-2 log L: 21.357
223 
-2 log L: 21.357
224 
> summary(mle(mll, lower=list(sigma = 0.01))) # alternative specification
 224 
> summary(mle(mll_lm, lower=list(sigma = 0.01))) # alternative specification
225 
Maximum likelihood estimation
 225 
Maximum likelihood estimation
226 
 
 226 
 
227 
Call:
 227 
Call:
228 
mle(minuslogl = mll, lower = list(sigma = 0.01))
 228 
mle(minuslogl = mll_lm, lower = list(sigma = 0.01))
229 
 
 229 
 
230 
Coefficients:
 230 
Coefficients:
231 
              Estimate Std. Error
 231 
              Estimate Std. Error
232 
b.(Intercept) 0.092667   0.480869
 232 
b.(Intercept) 0.092667   0.480869
233 
b.x           0.946061   0.077499
 233 
b.x           0.946061   0.077499
234 
sigma         0.703919   0.157400
 234 
sigma         0.703919   0.157400
235 
 
 235 
 
236 
-2 log L: 21.357
236 
-2 log L: 21.357
237 
>
237 
>
238 
> confint(mle(mll, lower=list(sigma = 0.01)))
 238 
> confint(mle(mll_lm, lower=list(sigma = 0.01)))
239 
Profiling...
 239 
Profiling...
240 
                 2.5 % 97.5 %
 240 
                 2.5 % 97.5 %
241 
b.(Intercept) -0.94831 1.1336
 241 
b.(Intercept) -0.94831 1.1336
242 
b.x            0.77829 1.1138
 242 
b.x            0.77829 1.1138
243 
sigma          0.48017 1.1755
 243 
sigma          0.48017 1.1755
244 
> plot(profile(mle(mll, lower=list(sigma = 0.01))))
 244 
> plot(profile(mle(mll_lm, lower=list(sigma = 0.01))))
245 
>
245 
>
246 
> Binom_mll <- function(x, n)
 246 
> Binom_mll <- function(x, n)
247 
+ {
 247 
+ {
248 
+     force(x); force(n) ## beware lazy evaluation
 248 
+     force(x); force(n) ## beware lazy evaluation
249 
+     function(p=.5) -dbinom(x, n, p, log=TRUE)
 249 
+     function(p=.5) -dbinom(x, n, p, log=TRUE)
Line 327... Line 327...
327 
> ### * <FOOTER>
 327 
> ### * <FOOTER>
328 
> ###
 328 
> ###
329 
> cleanEx()
 329 
> cleanEx()
330 
> options(digits = 7L)
 330 
> options(digits = 7L)
331 
> base::cat("Time elapsed: ", proc.time() - base::get("ptime", pos = 'CheckExEnv'),"\n")
 331 
> base::cat("Time elapsed: ", proc.time() - base::get("ptime", pos = 'CheckExEnv'),"\n")
332 
Time elapsed:  0.62 0.024 0.661 0 0
332 
Time elapsed:  0.574 0.015 0.593 0 0
333 
> grDevices::dev.off()
 333 
> grDevices::dev.off()
334 
null device
334 
null device
335 
          1
335 
          1
336 
> ###
 336 
> ###
337 
> ### Local variables: ***
 337 
> ### Local variables: ***