How can a make the sum of two step functions (R-stepfun) of class "stepfun"?

Question

From the example here I tried to make the sum as class "stepfun". I thought, as.stepfun is the right choice, but my ideas don't work. What is wrong?

y1 <- c(0, 1, 2, 0)
x1 <- c(1, 2, 3)
f1 <- stepfun(x = x1, y = y1)
print(class(f1))
# [1] "stepfun"  "function" #  OK!!!
plot(f1)

y2 <- c(0, 1, 0)
x2 <- c(1.5, 2.5)
f2 <- stepfun(x = x2, y = y2)
plot(f2)

fs <- function(x, f1, f2) {
  # y <- f1(x) + f2(x) #  OK
  # y <- as.stepfun(x = x, y = y, ties = "ordered", right = FALSE) # does not work
  # return(y) #           does not work
  return(f1(x) + f2(x))
}
print(class(fs)) # [1] "function"
# attributes(fs) # no new information...

fm <- function(x, f1, f2) {
  return(f1(x) * f2(x))
}
print(class(fm)) # [1] "function"

---

Example as. for data.frame which works as expected:

z <- c(1, 2)
class(z) #    [1] "numeric"
class(as.data.frame(z)) #    [1] "data.frame"

---

About internals of stepfun

function (x, y, f = as.numeric(right), ties = "ordered", right = FALSE) 
{
  if (is.unsorted(x)) 
    stop("stepfun: 'x' must be ordered increasingly")
  n <- length(x)
  if (n < 1) 
    stop("'x' must have length >= 1")
  n1 <- n + 1L
  if (length(y) != n1) 
    stop("'y' must be one longer than 'x'")
  rval <- approxfun(x, y[-if (right) 
    n1
  else 1], method = "constant", yleft = y[1L], yright = y[n1], 
    f = f, ties = ties)
  class(rval) <- c("stepfun", class(rval))
  attr(rval, "call") <- sys.call()
  rval
}

Answer 1

It looks like user2554330 is right about creating stepfun objects directly from two other stepfun objects, but here is a workaround in case it's useful:

x12 <- sort(unique(c(x1, x2)))
y12 <- f1(c(x12[1] - 1, x12)) + f2(c(x12[1] - 1, x12))
fs <- stepfun(x = x12, y = y12)

UPDATE1:

Or if you want to create the function from f1 and f2:

fs <- function(f1, f2) {
  xs <- sort(unique(c(get("x", envir = environment(f1)), get("x", envir = environment(f2)))))
  ys <- f1(c(xs[1] - 1, xs)) + f2(c(xs[1] - 1, xs))
  return(stepfun(x = xs, y = ys))
}
plot(fs(f1, f2))

UPDATE2:

Or for an arbitrarily long list of stepfun objects and a the ability to specify the combining function:

y3 <- c(-1, 2, 1, 0)
x3 <- c(0.5, 1, 3)
f3 <- stepfun(x = x3, y = y3)

fs <- function(lf, fun) {
  xs <- sort(unique(unlist(lapply(lf, function(x) get("x", envir = environment(x))))))
  ys <- apply(mapply(function(f) f(c(xs[1] - 1, xs)), lf), 1, FUN = fun)
  return(stepfun(x = xs, y = ys))
}
plot(fs(list(f1, f2, f3), "sum"))
plot(fs(list(f1, f2, f3), "prod"))

Answer 2

Thanks to the answers from @jblood94, @user2554330 and @rbm here I found an elegant way which I plan to use in my case. I hope that also helps others:

par(mfrow = c(2, 2))
y1 <- c(0, 1, 2, 0)
x1 <- c(1, 2, 3)
f1 <- stepfun(x = x1, y = y1)

y2 <- c(0, 1, 0)
x2 <- c(1.5, 2.5)
f2 <- stepfun(x = x2, y = y2)

plot(f1)
plot(f2)

'+.stepfun' <- function(f1, f2) {
  xs1 <- get("x", envir = environment(f1))
  xs2 <- get("x", envir = environment(f2))
  xs <- sort(unique(c(x1, x2)))
  ys <- f1(c(xs[1] - 1, xs)) + f2(c(xs[1] - 1, xs))
  return(stepfun(x = xs, y = ys))
}
f1 + f2
print(class(f1 + f2))
plot(f1 + f2, main = "Sum f1+f2")

'*.stepfun' <- function(f1, f2) {
  xs1 <- get("x", envir = environment(f1))
  xs2 <- get("x", envir = environment(f2))
  xs <- sort(unique(c(x1, x2)))
  ys <- f1(c(xs[1] - 1, xs)) * f2(c(xs[1] - 1, xs))
  return(stepfun(x = xs, y = ys))
}
f1 * f2
print(class(f1 * f2))
plot(f1 * f2, main = "Sum f1*f2")

par(mfrow = c(1, 1))

Answer 3

The fs function needs the "stepfun" class, not the results it returns. But your fs object won't work as a "stepfun" object, because R makes assumptions about those: they need to keep copies of the data that produced them, among other things. You can see what f1 keeps by looking at ls(environment(f1)). I don't know how those objects are used, but presumably they are needed.

Edited to add:

To turn fs into a "stepfun" object, you could try as.stepfun(fs). But this fails, with the error message

  Error in as.stepfun.default(fs) : 
no 'as.stepfun' method available for 'x'

The error message isn't the best (x is the internal name of the first argument, and in my opinion it would make more sense to say there's no method available for fs), but it's saying R doesn't know how to do the conversion you want.