r - Efficient apply or mapply for multiple matrix arguments by row -

-

i have 2 matrices want apply function to, rows:

matrixa            gsm83009  gsm83037  gsm83002  gsm83029  gsm83041 100001_at  5.873321  5.416164  3.512227  6.064150  3.713696 100005_at  5.807870  6.810829  6.105804  6.644000  6.142413 100006_at  2.757023  4.144046  1.622930  1.831877  3.694880  matrixb           gsm82939 gsm82940 gsm82974 gsm82975 100001_at 3.673556 2.372952 3.228049 3.555816 100005_at 6.916954 6.909533 6.928252 7.003377 100006_at 4.277985 4.856986 3.670161 4.075533

i've found several similar questions, not whole lot of answers: mapply matrices, multi matrix row-wise mapply?. code have splits matrices row lists, having split makes rather slow , not faster loop, considering have 9000 rows in each matrix:

scores <- mapply(t.test.stat, split(matrixa, row(matrixa)), split(matrixb, row(matrixb)))

the function simple, finding t-value:

t.test.stat <- function(x, y) {     return( (mean(x) - mean(y)) / sqrt(var(x)/length(x) + var(y)/length(y)) ) }

splitting matrices isn't biggest contributor evaluation time.

set.seed(21) matrixa <- matrix(rnorm(5 * 9000), nrow = 9000) matrixb <- matrix(rnorm(4 * 9000), nrow = 9000)  system.time( scores <- mapply(t.test.stat,     split(matrixa, row(matrixa)), split(matrixb, row(matrixb))) ) #    user  system elapsed  #    1.57    0.00    1.58  sma <- split(matrixa, row(matrixa)) smb <- split(matrixb, row(matrixb)) system.time( scores <- mapply(t.test.stat, sma, smb) ) #    user  system elapsed  #    1.14    0.00    1.14

look @ output rprof see of time is--not surprisingly--spent evaluating t.test.stat (mean, var, etc.). basically, there's quite bit of overhead function calls.

rprof() scores <- mapply(t.test.stat, sma, smb) rprof(null) summaryrprof()

you may able find faster generalized solutions, none approach speed of vectorized solution below.

since function simple, can take advantage of vectorized rowmeans function instantaneously (though it's bit messy):

system.time({ nca <- ncol(matrixa) ncb <- ncol(matrixb) ans <- (rowmeans(matrixa)-rowmeans(matrixb)) /   sqrt( rowmeans((matrixa-rowmeans(matrixa))^2)*(nca/(nca-1))/nca +         rowmeans((matrixb-rowmeans(matrixb))^2)*(ncb/(ncb-1))/ncb ) }) #    user  system elapsed  #      0       0       0  head(ans) # [1]  0.8272511 -1.0965269  0.9862844 -0.6026452 -0.2477661  1.1896181

update
here's "cleaner" version using rowvars function:

rowvars <- function(x, na.rm=false, dims=1l) {   rowmeans((x-rowmeans(x, na.rm, dims))^2, na.rm, dims)*(ncol(x)/(ncol(x)-1)) } ans <- (rowmeans(matrixa)-rowmeans(matrixb)) /   sqrt( rowvars(matrixa)/ncol(matrixa) + rowvars(matrixb)/ncol(matrixb) )

Comments

Post a Comment

Popular posts from this blog

php - What is the difference between $_SERVER['PATH_INFO'] and $_SERVER['ORIG_PATH_INFO']? -

-
what’s difference between $_server['path_info'] , $_server['orig_path_info'] ? how use them? when run print_r($_server) , path_info , orig_path_info not present in array. why not? how can enable them? i have read php manual on them, still don’t understand them. the path_info variable present if invoke php script this: http://www.example.com/phpinfo.php/hello_there it's /hello_there part after .php script. if don't invoke url that, there won't $_server["path_info"] environment variable. the porig_ prefix uncommon. path_info standard cgi-environment variable, , should never prefixed. did read that? (there issues around php3/php4 if invoked php interpreter via cgi-bin/ - hardly has such setups today.) for reference: http://www.ietf.org/rfc/rfc3875
Read more

fortran - Function return type mismatch -

-
i'm attempting recode old c++ program in fortran make use of lapack (i'm aware c++ have lapack++, i'm having lot of trouble installing it, gave up). i didn't have problems compilation, when had variables declared real . when started coding section of program required lapack, found parameters passed dsyev need double precision . tried change double precision (including changing hard coded numbers double precision counterparts, i.e. 0.0 -> 0.0d0) when try compile, following error of functions , subroutines i've written: error: return type mismatch of function <function> @ (1) (real(4)/real(8)) i'm not sure coming from, in program has been changed double precision. for example, i've declared following: double precision :: alpha(3),d(3),zeta1,zeta2 double precision :: a1(3),a2(3),d1(3),d2(3) double precision :: pi pi = 3.14159265359d0 alpha = (/0.109818d0, 0.405771d0, 2.22766d0/) d = (/0.444635d0, 0.535328d0, 0.154329d0 /) 10 i=1,3 ...
Read more

queue - mq_receive: message too long -

-
i implementing communication between 2 processes using queue. problem when call function mq_receive, error: message long. i have done following: struct mq_attr attr; long size = attr.mq_msgsize; .... // initializing queue "/gateway" int rc = mq_receive(gateway, buffer, size, &prio); if print size value, size=1, while when print same size program (got same mechanism), not long integer ( -1217186280 )... how can solve error?....so while size = 1, believe it's right "message long" why 1? p.s. have tried put : int rc = mq_receive(gateway, buffer, sizeof(buffer), &prio); but no result. it seems need read docs more carefully. when call mq_receive should pass size of destination buffer. this size must greater mq_msgsize attribute of queue . in addition, seems have error in queue attributes initialisation makes proper mq_receive call impossible. here standard message queue session: fill mq_attr struct ( doc ): struct mq_a...
Read more