[R] average matrices across a list

[Ben Bolker]

Gregory Ryslik <rsaber <at> comcast.net> writes:


> mymats <- vector('list', 5)
> set.seed(246)
> 
> # Generate a list of five 3 x 3 matrices
> for(i in 1:5) mymats[[i]] <- matrix(sample(1:9), nrow = 3)
> 
> mymats[[5]][1,1]<-NA
> mymats[[4]][2,2]<-NA
> mymats
> 
> matrixadder<-function(u,v){
> 	na.u<-is.na(u)
> 	na.v<-is.na(v)	
> 	ifelse(na.u & na.v, NA, ifelse(na.u, 0, u)+ ifelse(na.v,0,v))	
> }
> 
> Reduce('matrixadder',mymats)


  I was going to suggest that my solution would be faster,
but it turns out to be slower (!) -- 11 seconds vs 6 seconds
for 10,000 replications.

  It's cleverer, but I don't know if that's really a virtue.

library(abind)
apply(do.call(abind,c(mymats,list(along=3))),c(1,2),sum,na.rm=TRUE)

system.time(replicate(10000,
                      apply(do.call(abind,c(mymats,list(along=3))),
                            c(1,2),sum,na.rm=TRUE)))