tester <- function (y, k, func, grad) { qmaj <- function (x) func (y) + grad (y) * (x - y) + .5 * k * (x - y) ^ 2 ybnd <- y - 2 * grad (y) / k up <- max (y, ybnd) lw <- min (y, ybnd) x <- seq (lw, up, length = 100) s <- paste ("K = ", formatC(k, digits = 1, format= "f")) plot (x, func (x), col = "RED", lwd = 2, type = "l", ylab = "f,g", main = s) lines (x, qmaj (x), col = "BLUE", lwd = 2) abline (v = up) abline (v = lw) } f <- function (x) log (1 + exp (x)) g <- function (x) exp (x) / (1 + exp (x)) a <- function (x) (x ^ 3 - x) / 6 b <- function (x) (3 * x ^ 2 - 1) / 6 cubicSublevel <- function (y, a, sharp = FALSE) { f <- polynomial (a) g <- deriv (f) h <- deriv (g) i <- deriv (h) dfy <- predict (g, y) dgy <- predict (h, y) dhy <- predict (i, y) disk <- dgy ^ 2 - 4 * (2 / 3) * dhy * dfy if (disk <= 0) k <- max (0, dgy) else { r <- sort (solve (polynomial (c((2 / 3) * dhy * dfy, -dgy, 1)))) if ((r[1] >= 0) && (r[2] >= 0)) k <- dgy if ((r[1] <= 0) && (r[2] >= 0)) k <- r[2] if ((r[1] <= 0) && (r[2] <= 0)) stop ("unbounded") } return (y - dfy / k) }