>>>>> "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