FAQ
Is it possible to compute the variance of every column in a matrix
without a loop?
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !) To: r-help-request at stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._

## Search Discussions

•  at Mar 17, 2002 at 11:15 am ⇧

Francisco J Molina <fjmolina@lbl.gov> writes:

Is it possible to compute the variance of every column in a matrix
without a loop?
apply(m,2,var)

--
O__ ---- Peter Dalgaard Blegdamsvej 3
c/ /'_ --- Dept. of Biostatistics 2200 Cph. N
(*) \(*) -- University of Copenhagen Denmark Ph: (+45) 35327918
~~~~~~~~~~ - (p.dalgaard at biostat.ku.dk) FAX: (+45) 35327907
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !) To: r-help-request at stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
•  at Mar 17, 2002 at 11:45 am ⇧

Francisco J Molina wrote:
Is it possible to compute the variance of every column in a matrix
without a loop?
apply(X, 2, var)

Uwe Ligges
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !) To: r-help-request at stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
•  at Mar 17, 2002 at 11:57 am ⇧

On 03/16/02 19:43, Francisco J Molina wrote:
Is it possible to compute the variance of every column in a matrix
without a loop?
If the matrix is m1,

apply(m1,2,var)

or, with missing data,

apply(m1,2,var,na.rm=T)

In general, apply() and related functions are well worth
understanding.
--
Jonathan Baron, Professor of Psychology, University of Pennsylvania
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !) To: r-help-request at stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
•  at Mar 17, 2002 at 10:50 pm ⇧

On Sat, 16 Mar 2002, Francisco J Molina wrote:

Is it possible to compute the variance of every column in a matrix
without a loop?
No, but you can make it look as if there's no loop.

You can hide the for() loop inside a function with
apply(m, 2, var)
if the reason you want to avoid a loop is to make your code look simpler.

In 1.5.0 you will be able to do
colMeans(m*m)-colMeans(m)^2
if the reason is to have the loop take place in compiled code and be
faster.

-thomas

-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !) To: r-help-request at stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
•  at Mar 18, 2002 at 7:40 am ⇧

Thomas Lumley wrote:
On Sat, 16 Mar 2002, Francisco J Molina wrote:

Is it possible to compute the variance of every column in a matrix
without a loop?
No, but you can make it look as if there's no loop.

You can hide the for() loop inside a function with
apply(m, 2, var)
if the reason you want to avoid a loop is to make your code look simpler.

In 1.5.0 you will be able to do
colMeans(m*m)-colMeans(m)^2

..., but only for well conditioned problems, because using Steiner's
theorem is risky from the computational point of view. Try this extreme
example several times:

m <- rnorm(10, 1e8, 1e-7) # large mean, small variance
mean(m^2) - mean(m)^2 # typical results: -4, -2, 0, 2, 4

BTW: Of course this is a computational problem that is *not* related to
R!

if the reason is to have the loop take place in compiled code and be
faster.
Uwe Ligges
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !) To: r-help-request at stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
•  at Mar 18, 2002 at 3:28 pm ⇧

Francisco J Molina wrote:
Is it possible to compute the variance of every column in a matrix
without a loop?
1) apply(m, 2, var) (as several people said - simple, but not optimally fast)

2) colVars(m) from package "colSums", available at:
<http://cran.r-project.org/src/contrib/Devel/colSums_1.1.tar.gz>

3) In 1.5.0, I believe colSums and colMeans will be available (as fast C
routines). I don't know if Brian D. Ripley plans to include colVars too,
but if not the definition is simple enough:
-------------------------------------------------------------------------------
colVars <- function(x, na.rmúLSE, dims=1, unbiased=TRUE, SumSquaresúLSE,
twopassúLSE) {
if (SumSquares) return(colSums(x^2, na.rm, dims))
N <- colSums(!is.na(x), FALSE, dims)
Nm1 <- if (unbiased) N-1 else N
if (twopass) {x <- if (dims==length(dim(x))) x - mean(x, na.rm=na.rm) else
sweep(x, (dims+1):length(dim(x)), colMeans(x,na.rm,dims))}
(colSums(x^2, na.rm, dims) - colSums(x, na.rm, dims)^2/N) / Nm1
}
-------------------------------------------------------------------------------

With twopass=TRUE, this should act very much like the S-Plus function by the
same name. (twopass=TRUE is needed for numerical accuracy if mu >> sigma.)

For a matrix (as opposed to a >2-dimensional array), with unbiasedúLSE, you'd
get the same result as Thomas Lumley's suggestion:
R> colMeans(m*m)-colMeans(m)^2
However, this is *not* the same as apply(m, 2, var) because (sample) variances
are normally calculated with unbiased=TRUE (i.e. with N-1 in the denominator).
IF you have no NA's, then I'd modify Thomas's suggestion slightly:
R> N <- nrow(m)
R> N/(N-1) * (colMeans(m*m)-colMeans(m)^2)
--
-- David Brahm (brahm at alum.mit.edu)
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !) To: r-help-request at stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._

## Related Discussions

Discussion Overview
 group r-help categories r posted Mar 17, '02 at 3:43a active Mar 18, '02 at 3:28p posts 7 users 6 website r-project.org irc #r

### 6 users in discussion

Content

People

Support

Translate

site design / logo © 2017 Grokbase