Re: moving median in S without looping

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

Markus Loecher <Markus.Loecher@scr.siemens.com> wrote:
> does anyone know how to compute a moving median in Splus without using
> explicit looping ?

I wrote a C program and an interface to R that computes running medians (for N
odd) of matrix columns (also works for a single vector).  You'll have to modify
the interface for S-Plus.  You can find the full R package (also includes a
running mean, faster than "filter") at
  <http://cran.r-project.org/src/contrib/Devel/runStat_1.1.tar.gz>


-- (begin R code) --
runMedian <- function(m, N=5) {
  if (!(N %% 2)) stop("N should be odd.")
  d <- od <- dim(m)
  if (!(nd <- length(d))) d <- c(length(m), 1)               # Vector -> matrix
  if (nd > 2) d <- c(d[1], prod(d[-1]))                      # Array  -> matrix
  if (d[1] <= N) {m[] <- NA; return(m)}
  m <- array(.C("runMedian", m=as.real(m), as.integer(d), as.integer(N))$m, d)
  m[1:N, ] <- NA
  if (nd==2) m else if (!nd) m[ ,1] else array(m, od)     # Matrix/Vector/Array
}
-- (end R code) --


-- (begin C code) --
void runMedian(double *m, int *d, int *N) {
  int i, j, k, a, R=(*N+1)/2, r[*N];
  double old, new, *x;

  for (j=1; j <= d[1]; j++) {                           /* Loop over columns */
    x = m + j*d[0] - 1 - *N;                         /* x[*N] = m[nrow(m),j] */

    for (i=0; i < *N; i++) {                  /* Calculate initial ranks "r" */
      r[i] = 1;                               /* r higher if larger or later */
      for (k=0;   k < i;  k++) if (x[i] >= x[k]) r[i]++;
      for (k=i+1; k < *N; k++) if (x[i] >  x[k]) r[i]++;
      if (r[i]==R) x[*N] = x[i];            /* If rank=R, this is the median */
    }

    for (i=d[0]-1; i > *N; i--) {                      /* Now x[*N] = m[i,j] */
      x--;  old=x[*N];  new=x[0];  a=*N;              /* a = rank of new guy */
      for (k=*N-1; k > 0; k--) {          /* Recalculate each element's rank */
        r[k] = r[k-1];                   /* Shift previous iteration's ranks */
        if (x[k] >= new) {r[k]++; a--;}     /* Are we adding a lower number? */
        if (x[k] >  old) {r[k]--;}       /* Are we dropping a higher number? */
        if (r[k]==R) x[*N] = x[k];          /* If rank=R, this is the median */
      }
      r[0] = a;
      if (a==R) x[*N] = new;                   /* Is the new guy the median? */
    }
  }
}
-- (end C code) --

                              -- David Brahm (brahm@alum.mit.edu)