必ず前後を含めたメディアンで置き換えてしまう関数。

median.filter.univariate <- function(sequence, width.of.wings) {
  seq.o <- sequence # o:original
  seq.n <- seq.o    # n:new
  n.seq <- length(sequence)
  if( floor(width.of.wings)==width.of.wings ) {
    if( floor(width.of.wings/2)*2==width.of.wings ) {
      stop(message="the width of wings should be odd.")
    } else {
      n.wing <- floor(width.of.wings/2)
    }
  } else {
    stop(message="the width of wings should be integer (and also odd.")
  }
  for( i in c(1:n.seq) ) {
    wing <- c((i-n.wing):(i+n.wing))
    wing <- wing[wing>=1]      # Left edge of the sequence
    wing <- wing[wing<=n.seq]  # Right edge of the sequence
    seq.n[i] <- median(seq.o[wing])
  }
  return(seq.n)
}

閾値を超えたジャンプがあったときだけ、前後を含めたメディアンで置き換える関数。

median.filter.univariate.with.threshold <- function(sequence, 
                                                    width.of.wings,
                                                    filter.threshold) {
  seq.o <- sequence    # o:original
  seq.n <- seq.o       # n:new
  seq.d <- diff(sequence) # d:differential
  n.seq <- length(sequence)
  if( floor(width.of.wings)==width.of.wings ) {
    if( floor(width.of.wings/2)*2==width.of.wings ) {
      stop(message="the width of wings should be odd.")
    } else {
      n.wing <- floor(width.of.wings/2)
    }
  } else {
    stop(message="the width of wings should be integer (and also odd.")
  }
  for( i in c(2:n.seq) ) {
    if( abs(seq.d[i-1])>=filter.threshold ) {
      wing <- c((i-n.wing):(i+n.wing))
      wing <- wing[wing>=1]      # Left edge of the sequence
      wing <- wing[wing<=n.seq]  # Right edge of the sequence
      seq.n[i] <- median(seq.o[wing])
    }
  }
  return(seq.n)
}

デモ用のコード。

test.filter <- function(n) {
  X <- rnorm(n)
  plot(X)
  lines(median.filter.univariate(X,3), col="black")
  lines(median.filter.univariate.with.threshold(X,3,1.5), col="black", lty=2)
  lines(median.filter.univariate(X,5), col="blue")
  lines(median.filter.univariate.with.threshold(X,5,1.5), col="blue", lty=2)
  lines(median.filter.univariate(X,7), col="green")
  lines(median.filter.univariate.with.threshold(X,7,1.5), col="green", lty=2)
}

test.filter(100)