# --------------------------------------------------- # Program: Smoothing1.S # Author: Steven M. Boker # Date: Thu Oct 28 12:37:36 EST 2004 # # Smoothing Part 1 # # --------------------------------------------------- # --------------------------------------------------- # Create a time series tX <- seq(0, 10, by=1/24) tY <- sin(tX) + rnorm(length(tX), mean=0, sd=1.5) tSeries <- data.frame(x=tX,y=tY) # --------------------------------------------------- # Plot the timeseries graphsheet(height=6.4,width=7.5) xyplot(y ~ x, data = tSeries, panel = function(x, y) { panel.xyplot(x, y) panel.abline(lm(y ~ x, data = tSeries)) }, aspect=.7, xlab="Elapsed time in days", ylab="Score") # --------------------------------------------------- # Plot the timeseries using the plot command graphsheet(height=6.4,width=7.5) plot(c(0,10), c(-4.5, 4.5), xlab="Elapsed time in days", ylab="Score", type="n") lines(tX, tY, type="l", cex=.5) # --------------------------------------------------- # Calculate a running average windowWidth <- 10 theLength <- length(tX) - (windowWidth - 1) theMeanSmooth <- rep(NA, theLength) theCenter <- rep(NA, theLength) for (i in 1:theLength) { theMeanSmooth[i] <- mean(tY[i:(i + windowWidth - 1)]) theCenter[i] <- tX[i + floor((windowWidth + 1) / 2)] } # --------------------------------------------------- # Plot the timeseries and a running average graphsheet(height=6.4,width=7.5) plot(c(0,10), c(-4.5, 4.5), xlab="Elapsed time in days", ylab="Score", type="n") lines(tX, tY, type="p", cex=.5) lines(theCenter, theMeanSmooth, type="l", col=5, lwd=2) # --------------------------------------------------- # Calculate a kernel convolution smooth using a box filter windowWidth <- 15 theKernel <- rep(1/windowWidth, windowWidth) theLength <- length(tX) - (windowWidth - 1) theMeanSmooth <- rep(NA, theLength) theCenter <- rep(NA, theLength) for (i in 1:theLength) { theMeanSmooth[i] <- sum(tY[i:(i + windowWidth - 1)] * theKernel) theCenter[i] <- tX[i + floor((windowWidth + 1) / 2)] } # --------------------------------------------------- # Plot the timeseries and a running average graphsheet(height=6.4,width=7.5) plot(c(0,10), c(-4.5, 4.5), xlab="Elapsed time in days", ylab="Score", type="n") lines(tX, tY, type="p", cex=.5) lines(theCenter, theMeanSmooth, type="l", col=5, lwd=2) # --------------------------------------------------- # A step function timeseries tStepY <- c(rep(0,50), rep(1,50), rep(-1,50), rep(0,50)) tStepX <- c(1:200) # --------------------------------------------------- # Calculate a kernel convolution smooth using a box filter windowWidth <- 15 halfWidth <- 7 theKernel <- rep(1/windowWidth, windowWidth) theLength <- length(tStepX) - (windowWidth - 1) theMeanSmooth <- rep(NA, theLength) theCenter <- rep(NA, theLength) for (i in (1+halfWidth):(theLength-halfWidth)) { theMeanSmooth[i] <- sum(tStepY[(i-halfWidth):(i+halfWidth)] * theKernel) theCenter[i] <- tStepX[i] } # --------------------------------------------------- # Plot the timeseries and the smooth graphsheet(height=6.4,width=7.5) plot(c(0,200), c(-1.25, 1.25), xlab="Samples", ylab="Score", type="n") lines(tStepX, tStepY, type="p", cex=.5) lines(theCenter, theMeanSmooth, type="l", col=5, lwd=2) # --------------------------------------------------- # Calculate repeated convolutions theMeanSmooth2 <- rep(NA, theLength) theCenter2 <- rep(NA, theLength) for (i in (1+halfWidth):(theLength-halfWidth)) { theMeanSmooth2[i] <- sum(theMeanSmooth[(i-halfWidth):(i+halfWidth)] * theKernel) theCenter2[i] <- theCenter[i] } theMeanSmooth3 <- rep(NA, theLength) theCenter3 <- rep(NA, theLength) for (i in (1+halfWidth):(theLength-halfWidth)) { theMeanSmooth3[i] <- sum(theMeanSmooth2[(i-halfWidth):(i+halfWidth)] * theKernel) theCenter3[i] <- theCenter2[i] } # --------------------------------------------------- # Plot the timeseries and the smooths graphsheet(height=6.4,width=7.5) plot(c(0,200), c(-1.25, 1.25), xlab="Samples", ylab="Score", type="n") lines(tStepX, tStepY, type="p", cex=.5) lines(theCenter, theMeanSmooth, type="l", col=5, lwd=2) lines(theCenter2, theMeanSmooth2, type="l", col=6, lwd=2) lines(theCenter3, theMeanSmooth3, type="l", col=7, lwd=2) # --------------------------------------------------- # Calculate a kernel convolution smooth using a Spencer's filter # on the Step Function data. windowWidth <- 15 halfWidth <- 7 theKernel <- c(-3,-6, -5, 3, 21, 46, 67, 74, 67, 46, 21, 3, -5, -6, -3)/320 theLength <- length(tStepX) - (windowWidth - 1) theMeanSmooth <- rep(NA, theLength) theCenter <- rep(NA, theLength) for (i in (1+halfWidth):(theLength-halfWidth)) { theMeanSmooth[i] <- sum(tStepY[(i-halfWidth):(i+halfWidth)] * theKernel) theCenter[i] <- tStepX[i] } # --------------------------------------------------- # Plot the timeseries and smooth graphsheet(height=6.4,width=7.5) plot(c(0,200), c(-1.25, 1.25), xlab="Samples", ylab="Score", type="n") lines(tStepX, tStepY, type="p", cex=.5) lines(theCenter, theMeanSmooth, type="l", col=5, lwd=2) # --------------------------------------------------- # Calculate a kernel convolution smooth using a Spencer's filter # on the noisy sine wave data. windowWidth <- 15 halfWidth <- 7 theKernel <- c(-3,-6, -5, 3, 21, 46, 67, 74, 67, 46, 21, 3, -5, -6, -3)/320 theLength <- length(tX) - (windowWidth - 1) theMeanSmooth <- rep(NA, theLength) theCenter <- rep(NA, theLength) for (i in (1+halfWidth):(theLength-halfWidth)) { theMeanSmooth[i] <- sum(tY[(i-halfWidth):(i+halfWidth)] * theKernel) theCenter[i] <- tX[i] } # --------------------------------------------------- # Plot the timeseries and smooth graphsheet(height=6.4,width=7.5) plot(c(0,10), c(-4.5, 4.5), xlab="Elapsed time in days", ylab="Score", type="n") lines(tX, tY, type="p", cex=.5) lines(theCenter, theMeanSmooth, type="l", col=5, lwd=2) # --------------------------------------------------- # Calculate a running median windowWidth <- 30 theLength <- length(tX) - (windowWidth - 1) theMedianSmooth <- rep(NA, theLength) theCenter <- rep(NA, theLength) for (i in 1:theLength) { theMedianSmooth[i] <- median(tY[i:(i + windowWidth - 1)]) theCenter[i] <- tX[i + floor((windowWidth + 1) / 2)] } # --------------------------------------------------- # Plot the timeseries and a running median graphsheet(height=6.4,width=7.5) plot(c(0,10), c(-4.5, 4.5), xlab="Elapsed time in days", ylab="Score", type="n") lines(tX, tY, type="p", cex=.5) lines(theCenter, theMedianSmooth, type="l", col=5, lwd=2) # --------------------------------------------------- # Calculate a running median using smooth() theSmooth <- smooth(tY, twice=F) # --------------------------------------------------- # Plot the timeseries and the running median graphsheet(height=6.4,width=7.5) plot(c(0,10), c(-4.5, 4.5), xlab="Elapsed time in days", ylab="Score", type="n") lines(tX, tY, type="p", cex=.5) lines(tX, theSmooth, type="l", col=5, lwd=2) # --------------------------------------------------- # Calculate a running median using smooth() theSmooth <- smooth(tY, twice=T) # --------------------------------------------------- # Plot the timeseries and the smooth graphsheet(height=6.4,width=7.5) plot(c(0,10), c(-4.5, 4.5), xlab="Elapsed time in days", ylab="Score", type="n") lines(tX, tY, type="p", cex=.5) lines(tX, theSmooth, type="l", col=5, lwd=2) # --------------------------------------------------- # Calculate a smoothing spline theSmooth <- smooth.spline(tX, tY, df=5)$y # --------------------------------------------------- # Plot the timeseries and the smooth graphsheet(height=6.4,width=7.5) plot(c(0,10), c(-4.5, 4.5), xlab="Elapsed time in days", ylab="Score", type="n") lines(tX, tY, type="p", cex=.5) lines(tX, theSmooth, type="l", col=5, lwd=2) # --------------------------------------------------- # Plot the timeseries and the Spline for several different DF graphsheet(height=6.4,width=7.5) plot(c(0,10), c(-4.5, 4.5), xlab="Elapsed time in days", ylab="Score", type="n") lines(tX, tY, type="p", cex=.5) lines(tX, smooth.spline(tX, tY, df=32)$y, type="l", col=5, lwd=2) lines(tX, smooth.spline(tX, tY, df=16)$y, type="l", col=6, lwd=2) lines(tX, smooth.spline(tX, tY, df=8)$y, type="l", col=7, lwd=2) lines(tX, smooth.spline(tX, tY, df=4)$y, type="l", col=8, lwd=2) # --------------------------------------------------- # Quit the program