minQuadratic <- function (y, f, g, k, lw, up) { func <- function (x, f, g, k) { return (f + g * (x - y) + 0.5 * k * (x - y) ^ 2) } fup <- func (up, f, g, k) flw <- func (lw, f, g, k) if (k <= 0) { if (fup <= flw) return (list (x = up, f = fup)) if (fup >= flw) return (list (x = lw, f = flw)) } xmn <- y - g / k fmn <- func (xmn, f, g, k) if (xmn >= up) return (list (x = up, f = fup)) if (xmn <= lw) return (list (x = lw, f = flw)) return (list (x = xmn, f = fmn)) } minimaxQ <- function (y, f, g, k, lw, up) { n <- length (f) res <- Inf for (i in 1:n) { xlw <- lw xup <- up fail <- FALSE for (j in 1:n) { if (j == i) next df <- f[i] - f[j] dg <- g[i] - g[j] if (dg > 0) xlw <- max (xlw, y - df / dg) if (dg < 0) xup <- min (xup, y - df / dg) if ((dg == 0) && (f[i] < f[j])) fail <- TRUE } if (xlw >= xup) fail <- TRUE if (!fail) h <- minQuadratic (y, f[i], g[i], k, xlw, xup) if (h\$f < res) { res <- h\$f sol <- h\$x } } return (sol) } minimaxB <- function (y, f, g, k, lw, up) { n <- length (f) for (i in 1:(n - 1)) { for (j in (i + 1):n) { df <- f[i] - f[j] dg <- g[i] - g[j] dk <- k[i] - k[j] rt <- y + solve (polynomial (c (df, dg, dk / 2))) rt <- rt[which ((rt <= up) & (rt >= lw) & (!is.complex (rt)))] } } rt <- sort (unique (c(lw, up, rt))) m <- length (rt) fmin <- Inf for (j in 1:(m - 1)) { x <- (rt[j] + rt[j + 1]) / 2 smax <- -Inf for (i in 1:n) { si <- f[i] + g[i] * (x - y) + 0.5 * k[i] * (x - y) ^ 2 if (si > smax) { smax <- si l <- i } } h <- minQuadratic (y, f[l], g[l], k[l], rt[j], rt[j + 1]) if (h\$f < fmin) { fmin <- h\$f xmin <- h\$x } } return (xmin) }