The R Project SVN R

Rev

Rev 59039 | Rev 67599 | Go to most recent revision | Show entire file | Ignore whitespace | Details | Blame | Last modification | View Log | RSS feed

Rev 59039 Rev 61162
Line 44... Line 44...
44 
od <- options(digits = 5)
 44 
od <- options(digits = 5)
45 
x <- 0:10
 45 
x <- 0:10
46 
y <- c(26, 17, 13, 12, 20, 5, 9, 8, 5, 4, 8)
 46 
y <- c(26, 17, 13, 12, 20, 5, 9, 8, 5, 4, 8)
47 
 
 47 
 
48 
## Easy one-dimensional MLE:
 48 
## Easy one-dimensional MLE:
49 
nLL <- function(lambda) -sum(stats::dpois(y, lambda, log=TRUE))
 49 
nLL <- function(lambda) -sum(stats::dpois(y, lambda, log = TRUE))
50 
fit0 <- mle(nLL, start = list(lambda = 5), nobs = NROW(y))
 50 
fit0 <- mle(nLL, start = list(lambda = 5), nobs = NROW(y))
51 
# For 1D, this is preferable:
 51 
# For 1D, this is preferable:
52 
fit1 <- mle(nLL, start = list(lambda = 5), nobs = NROW(y),
 52 
fit1 <- mle(nLL, start = list(lambda = 5), nobs = NROW(y),
53 
            method = "Brent", lower = 1, upper = 20)
 53 
            method = "Brent", lower = 1, upper = 20)
54 
stopifnot(nobs(fit0) == length(y))
 54 
stopifnot(nobs(fit0) == length(y))
Line 61... Line 61...
61 
}
 61 
}
62 
(fit <- mle(ll, nobs = length(y)))
 62 
(fit <- mle(ll, nobs = length(y)))
63 
mle(ll, fixed = list(xhalf = 6))
 63 
mle(ll, fixed = list(xhalf = 6))
64 
## alternative using bounds on optimization
 64 
## alternative using bounds on optimization
65 
ll2 <- function(ymax = 15, xhalf = 6)
 65 
ll2 <- function(ymax = 15, xhalf = 6)
66 
    -sum(stats::dpois(y, lambda=ymax/(1+x/xhalf), log = TRUE))
 66 
    -sum(stats::dpois(y, lambda = ymax/(1+x/xhalf), log = TRUE))
67 
mle(ll2, method = "L-BFGS-B", lower = rep(0, 2))
 67 
mle(ll2, method = "L-BFGS-B", lower = rep(0, 2))
68 
 
 68 
 
69 
AIC(fit)
 69 
AIC(fit)
70 
BIC(fit)
 70 
BIC(fit)
71 
 
 71 
 
Line 77... Line 77...
77 
 
 77 
 
78 
## Use bounded optimization
 78 
## Use bounded optimization
79 
## The lower bounds are really > 0,
 79 
## The lower bounds are really > 0,
80 
## but we use >=0 to stress-test profiling
 80 
## but we use >=0 to stress-test profiling
81 
(fit2 <- mle(ll, method = "L-BFGS-B", lower = c(0, 0)))
 81 
(fit2 <- mle(ll, method = "L-BFGS-B", lower = c(0, 0)))
82 
plot(profile(fit2), absVal=FALSE)
 82 
plot(profile(fit2), absVal = FALSE)
83 
 
 83 
 
84 
## a better parametrization:
 84 
## a better parametrization:
85 
ll3 <- function(lymax = log(15), lxhalf = log(6))
 85 
ll3 <- function(lymax = log(15), lxhalf = log(6))
86 
    -sum(stats::dpois(y, lambda=exp(lymax)/(1+x/exp(lxhalf)), log=TRUE))
 86 
    -sum(stats::dpois(y, lambda = exp(lymax)/(1+x/exp(lxhalf)), log = TRUE))
87 
(fit3 <- mle(ll3))
 87 
(fit3 <- mle(ll3))
88 
plot(profile(fit3), absVal = FALSE)
 88 
plot(profile(fit3), absVal = FALSE)
89 
exp(confint(fit3))
 89 
exp(confint(fit3))
90 
 
 90 
 
91 
options(od)
 91 
options(od)