4

I want to avoid the following loop:

for(i in 1:2){
vectVal[i] = myFunc(M[,,i],S[,,i],phi2, sig2)
}

by using the apply function.

The problem is that the arguments passed to the apply function contain arrays (--> M and S) and scalars (--> phi2 and sig2).

I tried the following:

apply(M,3,myFunc, S = S, phi2 = phi2, sig2 = sig2)

which resulted in an error message because S is an array and not a matrix as required in myFunc (see below):

Here is a reproducible code:

M = array(data = c(
  0.5, 0.7, 0.45,
  0.5, 0.3, 0.45,
  0.5, 0.7, 0.3,
  0.5, 0.3, 0.7,
  0.5, 0.7, 0.45,
  0.5, 0.3, 0.55),
  dim = c(3,2,2),
)

S = array(data = c(
   0.7723229, -0.2149794, -0.2159068,
  -0.2149794,  0.7723229, -0.2083123,
  -0.2159068, -0.2083123,  0.7723229,
   0.7723229, -0.2149794, -0.2149794,
  -0.2149794,  0.7723229, -0.1783025,
  -0.2149794, -0.1783025,  0.7723229,
   0.7723229, -0.2149794, -0.2176665,
  -0.2149794,  0.7723229, -0.2111496,
  -0.2176665, -0.2111496,  0.7723229),
   dim = c(3,3,2)
)

phi2 = 0.5
sig2 = 0.3

myFunc = function(M, S, phi2, sig2){
  valMult = M[,1]%*%diag(S)
  valEnd = valMult + phi2 - sig2
return(valEnd)
}

vectVal = vector(length = 2)

for(i in 1:2){
  vectVal[i] = myFunc(M[,,i],S[,,i],phi2, sig2)
}

vectVal

Does someone has an idea?

user3118627
  • 41
  • 3
  • +1 for a example data and code, but this runs without error for me. Perhaps try it again in a clean R session. – Thomas Dec 19 '13 at 10:27
  • 1
    @Thomas it runs ok, they want an alternative to the `for` loop. Perhaps off-topic and belongs on code review. Also, `S` has 6 elements in the 3rd dimension. Do you really only want to use the first 2? – Simon O'Hanlon Dec 19 '13 at 10:28
  • @ Thomas: exactly, the for loop works, but not the apply-version – user3118627 Dec 19 '13 at 10:33
  • 1
    The exact equivalent of your for loop with *apply would be `sapply(1:2,function(x)myFunc(M[,,x],S[,,x],phi2,sig2))` but I don't really see the improvement. – plannapus Dec 19 '13 at 10:35
  • @SimonO101: My mistake, the 3rd dimension should have only 2 elements (I corrected the code). – user3118627 Dec 19 '13 at 10:36
  • @plannapus: Since the 3rd dimension of M and S in the original code is much bigger than 2, I thought the apply-version would require less time – user3118627 Dec 19 '13 at 10:39
  • @user3118627 I'm not so sure. You can try benchmarking say the `mapply` answer I give below but I am not sure it will be faster than a vanilla for loop. Try it. – Simon O'Hanlon Dec 19 '13 at 10:40
  • @user3118627 Take a look at this one : http://stackoverflow.com/questions/2275896/is-rs-apply-family-more-than-syntactic-sugar – Joris Meys Dec 19 '13 at 10:47

1 Answers1

3

One (not particularly efficient) way would be to use plyr to split your arrays into lists (each element of the lists are the third dimension of your arrays). You could then use mapply to run your function like so:

require( plyr)
ml <- alply( M , 3 )
sl <- alply( S , 3 )
mapply( myFunc , ml , sl , phi2 , sig2 )
#       1        2 
#1.474333 1.358484

Update:

A more vectorised alternative (but still not as fast as for and %*% [see @JorisMeys comment below]) is to get the diag of S and then use colSums and matrix multiplication like so to achieve the same result:

s <- apply(S,3,diag)
colSums( M[,1,] * s ) + phi2 - sig2
# [1] 1.474333 1.358484

Update, update:

@JorisMeys has written a vectorised extractor function for getting the diagonal elements of 3D square arrays. Check this out.

Community
  • 1
  • 1
Simon O'Hanlon
  • 58,647
  • 14
  • 142
  • 184
  • This is about 90 times slower than the for-loop. Although very useful, plyr is not the best choice when performance is an issue. – Joris Meys Dec 19 '13 at 10:45
  • mostly plyr, but mapply doesn't offer a speedup compared to the loop either. – Joris Meys Dec 19 '13 at 10:50
  • Take a look at the `rbenchmark` package. That will tell you that the loop in the question outperforms the rest, for the simple reason that it's rather impossible to beat `%*%`. OP will have to look at eg Rcpp if he wants a speedup. – Joris Meys Dec 19 '13 at 11:09
  • @JorisMeys yeah, I was being lazy not checking. Oh darn. – Simon O'Hanlon Dec 19 '13 at 11:12