# [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> 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 - 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):ceiling(u),
h = floor(u):ceiling(u), 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
```