[Rd] outer not applying a constant function
Martin Maechler
maechler at stat.math.ethz.ch
Mon Mar 20 14:36:42 CET 2017
>>>>> Gebhardt, Albrecht <Albrecht.Gebhardt at aau.at>
>>>>> on Sun, 19 Mar 2017 09:14:56 +0000 writes:
> Hi,
> the function outer can not apply a constant function as in the last line of the following example:
>> xg <- 1:4
>> yg <- 1:4
>> fxyg <- outer(xg, yg, function(x,y) x*y)
>> fconstg <- outer(xg, yg, function(x,y) 1.0)
> Error in outer(xg, yg, function(x, y) 1) :
> dims [product 16] do not match the length of object [1]
> Of course there are simpler ways to construct a constant matrix, that is not my point.
> It happens for me in the context of generating matrices of partial derivatives, and if on of these partial derivatives happens to be constant it fails.
> So e.g this works:
> library(Deriv)
> f <- function(x,y) (x-1.5)*(y-1)*(x-1.8)+(y-1.9)^2*(x-1.1)^3
> fx <- Deriv(f,"x")
> fy <- Deriv(f,"y")
> fxy <- Deriv(Deriv(f,"y"),"x")
> fxx <- Deriv(Deriv(f,"x"),"x")
> fyy <- Deriv(Deriv(f,"y"),"y")
> fg <- outer(xg,yg,f)
> fxg <- outer(xg,yg,fx)
> fyg <- outer(xg,yg,fy)
> fxyg <- outer(xg,yg,fxy)
> fxxg <- outer(xg,yg,fxx)
> fyyg <- outer(xg,yg,fyy)
> And with
> f <- function(x,y) x+y
> it stops working. Of course I can manually fix this for that special case, but thats not my point. I simply thought "outer" should be able to handle constant functions.
?outer clearly states that FUN needs to be vectorized
but function(x,y) 1 is not.
It is easy to solve by wrapping the function in Vectorize(.):
> x <- 1:3; y <- 1:4
> outer(x,y, function(x,y) 1)
Error in dim(robj) <- c(dX, dY) :
dims [product 12] do not match the length of object [1]
> outer(x,y, Vectorize(function(x,y) 1))
[,1] [,2] [,3] [,4]
[1,] 1 1 1 1
[2,] 1 1 1 1
[3,] 1 1 1 1
----------------
So, your "should" above must be read in the sense
"It really would be convenient here and
correspond to other "recycling" behavior of R"
and I agree with that, having experienced the same inconvenience
as you several times in the past.
outer() being a nice R-level function (i.e., no C speed up)
makes it easy to improve:
Adding something like the line
if(length(robj) == 1L) robj <- rep.int(robj, dX*dY)
before dim(robj) <- c(dX, dY) [which gave the error]
would solve the issue and not cost much (in the cases it is unneeded).
Or is this a bad idea?
More information about the R-devel
mailing list