What are the rules for threading a function over a vector in R?

user_2382517 0
Viewed 117

I have some code which I call with two vectors of different length, lets call them A and B. However, I wrote the function having in mind a single element of A with the expectation that it will be automatically threaded over A. To be concrete,

A <- rnorm(5)
B <- rnorm(30)

foo <- function(x,B){
    sum( cos(x*B) ) # calculate sum_i cos(x*B[i])
}

sum( exp(foo(A,B)) ) # expecting this to calculate the exponent for each A[j] and add over j

I need to get

Σ_j exp( Σ_i cos(A[j]*B[i])

and not

Σ_ij exp(cos(A[j]*B[i])) OR  exp(cos(Σ_ij A[j]*B[i]))

I suspect that the last R expression is ambiguous, since the declaration of foo does not know B is always a vector. What are the formal rules and am I right to worry about the ambiguity?

lapply r
2 Answers

If we want to loop over the 'A', then use sapply , and apply the foo on each of the elements of 'A' with anonymous function call and get the sum of the output vector

sum(exp(sapply(A, function(x) foo(x, B))))

In the OP's example with the expression foo(A, B), the product A*B is computed first, and since the lengths of A and B are unequal, the recycling rule takes priority. There is no error message coming out, just because by pure luck the vector length of one is a multiple of the other.

user_37322710

You can also Vectorize the x input. I think this is what you were expecting. At the end of the day, this will work it's way down to an mappy() implementation which is a multivariate sapply, so probably best to just do it yourself as with the solution from akrun.

foo2 <- Vectorize(foo, "x")

sum(exp(foo2(A, B)))

The "formal rules" as you put them is quite simply how R does help("Arithmetic").

The binary operators return vectors containing the result of the element by element operations. If involving a zero-length vector the result has length zero. Otherwise, the elements of shorter vectors are recycled as necessary (with a warning when they are recycled only fractionally). The operators are + for addition, - for subtraction, * for multiplication, / for division and ^ for exponentiation.

So when you use x*B, it is doing element-wise multiplication. Nothing changes when you pass A into the function instead of x.

Simply go through your lines one at a time.

  1. x*B will be a vector of length max(length(x, B)). When they are not of the same length, R will recycle elements of the shorter vector (i.e., repeat them).
  2. cos(x*B) will be a vector of the same length as step (1), but now the cosine of that value.
  3. sum( cos(x*B) ) will sum that vector, returning a single number.
  4. foo(A,B) does steps (1) through (3), but with your defined A and B. Note that in your example A is recycled 6 times to get to the length of B. In other words, what you entered as A is being used as rep(A, 6) in the multiplication step. Nothing about a function definition in R says that foo(A,B) should be repeated for each element of vector A. So it behaves literally as you wrote it, basically swapping in A for x in the function code.
  5. exp(foo(A,B)) will take the result from foo from step 3 (which is a scalar) and raise it to an exponent.
  6. sum( exp(foo(A,B)) ) does nothing, since step (5) is a scalar, there is nothing to sum.
user_109174790
Related