Re: R equivalent of Splus rowVars function

[David Brahm <brahm@alum.mit.edu>]

Mark Leeds <mleeds@mlp.com> wrote (to S-News):
> does anyone know the R equivalent of the SPlus rowVars function ?

Andy Liaw <andy_liaw@merck.com> replied:
> More seriously, I seem to recall David Brahms at one time had created an R
> package with these dimensional summary statistics, using C code.  (And I
> pointed him to the `two-pass' algorithm for variance.)

Here are the functions that didn't make it into R's base package, which should
be very similar to the S-Plus functions of the same names.  The "twopass"
argument determines whether Andy's two-pass algorithm (Chan Golub & LeVegue) is
used (it's slower but more accurate).  In real life I set the "twopass" default
to FALSE, because in finance noise is always bigger than signal.

I am cc'ing to R-help, as this is really an R question.
-- 
                              -- David Brahm (brahm@alum.mit.edu)


colVars <- function(x, na.rm=FALSE, dims=1, unbiased=TRUE, SumSquares=FALSE,
                    twopass=TRUE) {
  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
}

rowVars <- function(x, na.rm=FALSE, dims=1, unbiased=TRUE, SumSquares=FALSE,
                    twopass=TRUE) {
  if (SumSquares) return(rowSums(x^2, na.rm, dims))
  N <- rowSums(!is.na(x), FALSE, dims)
  Nm1 <- if (unbiased) N-1 else N
  if (twopass) {x <- if (dims==0) x - mean(x, na.rm=na.rm) else
                     sweep(x, 1:dims, rowMeans(x,na.rm,dims))}
  (rowSums(x^2, na.rm, dims) - rowSums(x, na.rm, dims)^2/N) / Nm1
}

colStdevs <- function(x, ...) sqrt(colVars(x, ...))

rowStdevs <- function(x, ...) sqrt(rowVars(x, ...))