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 84... Line 84...
84 
# For 1D, this is preferable:
 84 
# For 1D, this is preferable:
85 
fit1 <- mle(nLL, start = list(lambda = 5), nobs = NROW(ec50),
 85 
fit1 <- mle(nLL, start = list(lambda = 5), nobs = NROW(ec50),
86 
            method = "Brent", lower = 1, upper = 20)
 86 
            method = "Brent", lower = 1, upper = 20)
87 
 
 87 
 
88 
## This needs a constrained parameter space: most methods will accept NA
 88 
## This needs a constrained parameter space: most methods will accept NA
89 
ll <- function(ymax = 15, xhalf = 6) {
 89 
mll1 <- function(ymax = 15, xhalf = 6) {
90 
    if(ymax > 0 && xhalf > 0)
90 
    if(ymax > 0 && xhalf > 0)
91 
      with(ec50, -sum(stats::dpois(y, lambda = ymax/(1+x/xhalf), log = TRUE)))
 91 
      with(ec50, -sum(stats::dpois(y, lambda = ymax/(1+x/xhalf), log = TRUE)))
92 
    else NA
 92 
    else NA
93 
}
 93 
}
94 
(fit <- mle(ll, nobs = NROW(ec50)))
 94 
(fit <- mle(mll1, nobs = NROW(ec50)))
95 
mle(ll, fixed = list(xhalf = 6))
 95 
mle(mll1, fixed = list(xhalf = 6))
96 
 
 96 
 
97 
## Alternative using bounds on optimization
 97 
## Alternative using bounds on optimization
98 
ll2 <- function(ymax = 15, xhalf = 6)
 98 
mll2 <- function(ymax = 15, xhalf = 6)
99 
    with(ec50, -sum(stats::dpois(y, lambda = ymax/(1+x/xhalf), log = TRUE)))
 99 
    with(ec50, -sum(stats::dpois(y, lambda = ymax/(1+x/xhalf), log = TRUE)))
100 
mle(ll2, lower = rep(0, 2))
 100 
mle(mll2, lower = rep(0, 2))
101 
 
 101 
 
102 
AIC(fit)
 102 
AIC(fit)
103 
BIC(fit)
 103 
BIC(fit)
104 
 
 104 
 
105 
summary(fit)
 105 
summary(fit)
Line 109... Line 109...
109 
confint(fit)
 109 
confint(fit)
110 
 
 110 
 
111 
## Use bounded optimization
 111 
## Use bounded optimization
112 
## The lower bounds are really > 0,
 112 
## The lower bounds are really > 0,
113 
## but we use >=0 to stress-test profiling
 113 
## but we use >=0 to stress-test profiling
114 
(fit2 <- mle(ll2, lower = c(0, 0)))
 114 
(fit2 <- mle(mll2, lower = c(0, 0)))
115 
plot(profile(fit2), absVal = FALSE)
 115 
plot(profile(fit2), absVal = FALSE)
116 
 
 116 
 
117 
## A better parametrization:
 117 
## A better parametrization:
118 
ll3 <- function(lymax = log(15), lxhalf = log(6))
 118 
mll3 <- function(lymax = log(15), lxhalf = log(6))
119 
    with(ec50, -sum(stats::dpois(y, lambda = exp(lymax)/(1+x/exp(lxhalf)), log = TRUE)))
 119 
    with(ec50, -sum(stats::dpois(y, lambda = exp(lymax)/(1+x/exp(lxhalf)), log = TRUE)))
120 
(fit3 <- mle(ll3))
 120 
(fit3 <- mle(mll3))
121 
plot(profile(fit3), absVal = FALSE)
 121 
plot(profile(fit3), absVal = FALSE)
122 
exp(confint(fit3))
 122 
exp(confint(fit3))
123 
 
 123 
 
124 
# Regression tests for bounded cases (this was broken in R 3.x)
 124 
# Regression tests for bounded cases (this was broken in R 3.x)
125 
fit4 <- mle(ll, lower = c(0, 4)) # has max on boundary
 125 
fit4 <- mle(mll1, lower = c(0, 4)) # has max on boundary
126 
confint(fit4)
 126 
confint(fit4)
127 
 
 127 
 
128 
## direct check that fixed= and constraints work together
 128 
## direct check that fixed= and constraints work together
129 
mle(ll, lower = c(0, 4), fixed=list(ymax=23)) # has max on boundary
 129 
mle(mll1, lower = c(0, 4), fixed=list(ymax=23)) # has max on boundary
130 
 
 130 
 
131 
## Linear regression using MLE
 131 
## Linear regression using MLE
132 
lr <- data.frame(
 132 
lr <- data.frame(
133 
x = 1:10,
133 
x = 1:10,
134 
y = c(0.48, 2.24, 2.22, 5.15, 4.64, 5.53, 7, 8.8, 7.67, 9.23))
 134 
y = c(0.48, 2.24, 2.22, 5.15, 4.64, 5.53, 7, 8.8, 7.67, 9.23))
Line 142... Line 142...
142 
     names(b0) <- colnames(X)
 142 
     names(b0) <- colnames(X)
143 
     function(b=b0, sigma=1)
 143 
     function(b=b0, sigma=1)
144 
         -sum(dnorm(y, X \%*\% b, sigma, log=TRUE))
 144 
         -sum(dnorm(y, X \%*\% b, sigma, log=TRUE))
145 
}
 145 
}
146 
 
 146 
 
147 
mll <- LM_mll(y ~ x, lr)
 147 
mll_lm <- LM_mll(y ~ x, lr)
148 
 
 148 
 
149 
summary(lm(y~x, data=lr)) # for comparison -- notice variance bias in MLE
 149 
summary(lm(y~x, data=lr)) # for comparison -- notice variance bias in MLE
150 
summary(mle(mll, lower=c(-Inf,-Inf, 0.01)))
 150 
summary(mle(mll_lm, lower=c(-Inf,-Inf, 0.01)))
151 
summary(mle(mll, lower=list(sigma = 0.01))) # alternative specification
 151 
summary(mle(mll_lm, lower=list(sigma = 0.01))) # alternative specification
152 
 
 152 
 
153 
confint(mle(mll, lower=list(sigma = 0.01)))
 153 
confint(mle(mll_lm, lower=list(sigma = 0.01)))
154 
plot(profile(mle(mll, lower=list(sigma = 0.01))))
 154 
plot(profile(mle(mll_lm, lower=list(sigma = 0.01))))
155 
 
 155 
 
156 
Binom_mll <- function(x, n)
 156 
Binom_mll <- function(x, n)
157 
{
 157 
{
158 
    force(x); force(n) ## beware lazy evaluation
 158 
    force(x); force(n) ## beware lazy evaluation
159 
    function(p=.5) -dbinom(x, n, p, log=TRUE)
 159 
    function(p=.5) -dbinom(x, n, p, log=TRUE)