- ## compute the 95 % profile likelihood CIs
- ## (not very clean... but seems to work -- EP 2003-03-29)
- CI <- matrix(NA, 2, 2)
- foo <- function(p) dev(p, para[2]) - 3.84 - Dev
- inc <- 1e-2
- lo <- para[1] - inc
- up <- para[1] + inc
- while (foo(lo) < 0) lo <- lo - inc
- while (foo(up) < 0) up <- up + inc
- CI[1, 1] <- uniroot(foo, lower = lo, upper = para[1])$root
- if (CI[1, 1] < 0) CI[1, 1] <- 0
- CI[1, 2] <- uniroot(foo, lower = para[1], upper = up)$root
- foo <- function(p) dev(para[1], p) - 3.84 - Dev
- lo <- para[2] - inc
- up <- para[2] + inc
- while (foo(lo) < 0) lo <- lo - inc
- while (foo(up) < 0) up <- up + inc
- CI[2, 1] <- uniroot(foo, lower = lo, upper = para[2])$root
- CI[2, 2] <- uniroot(foo, lower = para[2], upper = up)$root
+
+ ## 95% profile likelihood CIs
+
+ ## which: index of the parameter (1 or 2)
+ ## s: sign of the increment (-1 or +1)
+ foo <- function(which, s) {
+ i <- 0.1
+
+ if (which == 1) {
+ p <- para[1] + s * i
+ bar <- function() dev(p, para[2])
+ } else { # which == 2
+ p <- para[2] + s * i
+ bar <- function() dev(para[1], p)
+ }
+
+ while (i > 1e-9) {
+ while (bar() < Dev + 3.84) p <- p + s * i
+ p <- p - s * i
+ i <- i / 10
+ }
+ p
+ }
+
+ CI <- mapply(foo, c(1, 2, 1, 2), c(-1, -1, 1, 1))
+ dim(CI) <- c(2, 2)
+