s-news
[Top] [All Lists]

Re: moving median in S without looping

To: "'brahm@alum.mit.edu'" <brahm@alum.mit.edu>, s-news@lists.biostat.wustl.edu
Subject: Re: moving median in S without looping
From: "Warnes, Gregory R" <gregory_r_warnes@groton.pfizer.com>
Date: Tue, 21 Jan 2003 11:56:04 -0500
Cc: Markus.Loecher@scr.siemens.com
And to complete the bevy of options, here is my code (part of the R gregmisc
package) for computing running values of arbitrary functions.  It also works
fine in S-Plus.

-Greg

--- start of code ---
"running" <- function( X, fun=mean, width=min(length(X),20),
                     allow.fewer=FALSE,...)
  running2(X=X, fun=fun, width=width, allow.fewer=allow.fewer, ...)

"running2" <- function( X, Y=NULL, fun=mean, width=min(length(X),20),
                     allow.fewer=FALSE,...)
{
  n <- length(X)

  from  <-  sapply( (1:n) - width + 1, function(x) max(x,1) )
  to    <-  1:n

  elements  <- apply(cbind(from,to), 1,function(x) seq(x[1], x[2]) )

  if(is.matrix(elements))
    elements  <- as.data.frame(elements)

  if(is.null(Y))  # univariate 
    {
      funct <- function(which,what,fun,...) fun(what[which],...)
      
      Xvar <- sapply(elements, funct, what=X, fun=fun, ...)
    }
  else # bivariate
    {
      funct <- function(which,XX,YY,fun,...) fun(XX[which],YY[which], ...)
      
      Xvar <- sapply(elements, funct, XX=X, YY=Y, fun=fun, ...)
    }

  
  names(Xvar) <- paste(from,to,sep=":")

  if(!allow.fewer)
    Xvar[1:(width-1)]  <- NA
  
  return(Xvar)
}
--- end of code ---

--- start of help documentation ---
Apply a Function Over Adjacent Subsets of a Vector

Description:

     Applies a function over subsets of the vector(s) formed by taking
     a fixed number of previous points.

Usage:

     running(X, fun=mean, width=min(length(X),20), allow.fewer=FALSE,...)
     running2(X,Y, fun=mean, width=min(length(X),20), allow.fewer=FALSE,...)

Arguments:

       X: data vector 

       Y: data vector 

     fun: function to apply. Default is `mean'

   width: integer giving the number of vector elements to include in
          the subsets.  Defaults to the lesser of the length of the
          data and 20 elements.

allow.fewer: Boolean indicating whether the function should be computed
          for initial subsets with fewer than `width' points

     ...: parameters to be passed to `fun' 

Details:

     `running' applies the specified univariate function to sequential
     windows on `X'.

     `running2' applies the specified bivariate function to sequential
     windows of `X' and `Y'.

Value:

     Vector containg the results of applying the function `fun' to the
     subsets of `X' (`running') or `X' and `Y' (running2).

Author(s):

     Gregory R. Warnes Gregory_R_Warnes@groton.pfizer.com

Examples:

     running(1:20,width=5)

     plot(1:20, running(1:20,width=5))
     plot(1:20, running(1:20,width=5, allow.fewer=TRUE))

     # plot running mean and central 2 standard deviation range
     # estimated by last 40 observations
     dat <- rnorm(500, sd=1 + (1:500)/500 )
     plot(dat)
     fun <- function(x,sign=1) mean(x) + sign * sqrt(var(x))
     lines(running(dat,width=50,fun=mean,allow=TRUE),col="blue")
     lines(running(dat,width=50,fun=fun, sign=-1, allow=TRUE),col="red")
     lines(running(dat,width=50,fun=fun, sign=1, allow=TRUE),col="red")

     # plot running running correlation estimated by last 40 observations
     # along with true (local) correlation
     X <- rnorm(500, sd=1)
     Y <- X + rnorm(500, sd=(1:500)/500)
     rho <- 1/ ( 1 * sqrt(1 + (1:500)/500) )  # true 

     plot(running2(X,Y,width=20,fun=cor),col="red",type="s")
     lines(rho,type="l",col="blue")

--- end of documentation ---



LEGAL NOTICE
Unless expressly stated otherwise, this message is confidential and may be 
privileged. It is intended for the addressee(s) only. Access to this E-mail by 
anyone else is unauthorized. If you are not an addressee, any disclosure or 
copying of the contents of this E-mail or any action taken (or not taken) in 
reliance on it is unauthorized and may be unlawful. If you are not an 
addressee, please inform the sender immediately.

<Prev in Thread] Current Thread [Next in Thread>