[R] optmizing with monotone stepfunctions?

Martin Maechler maechler at stat.math.ethz.ch
Wed Jan 10 14:52:18 CET 2001


>>>>> "Jens" == Jens Oehlschlaegel <jens.oehlschlaegel at bbdo-interone.de> writes:

    Jens> Before re-inventing the wheel I would like to ask: does anyone
    Jens> know about an optimizer in R which can reliably identify which
    Jens> value of X (Xopt) leads to Y (Yopt) closest to Ytarget in

    Jens> Y <- MonotoneStepFun(X)

    Jens> optionally with the restriction that Yopt <= Ytarget (at least if
    Jens> any Y <= Ytarget, otherwise any Yopt > Ytarget would be the
    Jens> preferred answer)

    Jens> If none is known, I will write one.

Do you mean monotone (also called "isotone" or "isotonic") regression?
It is known that the solution of the monotone regression problem
*is* a step function.

Now, MASS in the VR bundle does isotonic regression "inside" isoMDS().

I have recently written some code even using a class with a print and plot
method, relying on MASS' code.
I've actually wanted to propose to Bill and Brian to incorporate a form of
this into MASS 
-- or to allow to put their C code (and my R code) into "Standard R"... 

Here is my code : Function definitions and example calls 

-------------- next part --------------
### Simplification of MASS' Shepard()   [plot:]
##
## `` Given a dissimilarity d and configuration x, compute Shepard plot ''
###
### which is in /usr/local/app/R/R_local/src/VR/MASS/R/isoMDS.R
### C is in /usr/local/app/R/R_local/src/VR/MASS/src/MASS.c --> "VR_mds_"

## Note that `Shepard()' in MASS is not really documented ..
require(MASS)

isoreg <- function(x, y)
{
  ##-- strong simplification of MASS' Shepard()
  ord <- order(x) #- this is the only thing we need from `x'
  pn <- length(ord)
  y <- y[ord]
  Z <- .C("VR_mds_fn",
	  as.double(y),
	  yf = double(pn),
	  pn = pn,
	  ssq = double(1),
	  pd = as.integer(order(ord)-1),
          ## the next 4 args are just dummies here:
	  x = double(1),## as.double(x),
	  integer(1),	## "pr" as.integer(n),
	  integer(1),	## "pncol" as.integer(k),
          double (1),	## double(n*k),
          do.deriv = FALSE
          , DUP = FALSE, PACKAGE = "MASS"
	  )
  structure(list(x = x[ord], y = y, yf = Z$yf, ord = ord, call = match.call()),
            class = "isoreg")
}

print.isoreg <- function(x, digits = getOption("digits"), ...) {
  cat("Isotonic regression: ")#, x$call)
  str(x)
  invisible(x)
}

plot.isoreg <-
function(x, main = paste("Isotonic regression",deparse(x$call)),
         col.wise = FALSE,
         mar = .1 + c(3.5,2.5,1,1), 
         mgp = c(1.6, 0.7, 0),
         grid = length(x$x) < 12,
         col = "red",
         cex.main = par("cex.main"),
         col.main = par("col.main"),
         font.main = par("font.main"), ...) {
    if(!is.null(main)) main.wid <- 2
    op <- par(mfcol = if(col.wise) 1:2 else 2:1,
              oma = c(0,0, main.wid, 0), mar = mar, mgp = mgp) 
    on.exit(par(op))

    x0 <- x$x ; x0 <- c(x0[1] - min(diff(x0)), x0)
    cy <- cumsum(c(0, x$y))
    cf <- cumsum(c(0, x$yf))
    i <- cy == cf

    ## Plot of "Data" + Fit
    plot(x0, c(NA,x$y), ..., ylab = "x$y")
    lines (x$x, x$yf, col = col, lwd = 2, type = "S")
    points(x$x[i[-1]], x$yf[i[-1]], col = col, pch = 18, cex = 1.5)
    if(grid) {
        u <- par("usr") 
        abline(v = floor(u[1]):ceiling(u[2]),
               h = floor(u[3]):ceiling(u[4]), col = "gray70", lty = 3)
    }
    ## Cumulative Plot
    plot (x0, cy, type = "n", ylab = "cumsum(x$y)", ylim = range(cy,cf), ...)
    lines(x0, cf, col = col, lwd = 1.5)
    points(x0[i], cy[i], col = col, cex = 1.5, pch = 19)
    if(grid) abline(v = x0[i], col = "gray70", lty = 3, xpd = !col.wise)
    points(x0[-1], cy[-1])# over draw

    if(!is.null(main))
        mtext(main, side = 3, outer = TRUE,
              cex = cex.main, col = col.main, font = font.main)
    invisible()
}

s <- isoreg(1:9, c(1,0,4,3,3,5,4,2,0))
plot(s)
## `same' plot,  "proving" that only ranks of `x' are important
plot(isoreg(2^(1:9), c(1,0,4,3,3,5,4,2,0)), log = "x")

plot(s3 <- isoreg(   1:9,  y3 <- c(1,0,4,3,3,5,4,2, 3)))# last "3", not "0"

## experiment a bit (C-c C-j)
plot(isoreg(sample(9),  y3))

plot(ir <- isoreg(sample(10), sample(10, replace = TRUE)), col.wise = FALSE)

##> nice examples (use with x = 1:10):
plot(s4 <- isoreg(1:10, y4 <- c(5, 9, 1:2, 5:8, 3, 8)))
s4


More information about the R-help mailing list