R help - Mar 2002 - compute variance of every column in a matrix without a loop

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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._

Francisco J Molina <FJMolina at 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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._

Francisco J Molina wrote:
> 
> Is it possible to compute the variance of every column in a matrix
> without a loop?
What about
 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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._

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
Home page: http://www.sas.upenn.edu/~baron
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._

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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._

Francisco J Molina <FJMolina at lbl.gov> 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=FALSE, dims=1, unbiased=TRUE, SumSquares=FALSE,
                    twopass=FALSE) {
  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=FALSE,
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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._