# --------------------------------------------------- # Program: Smoothing1.R # Author: Steven M. Boker # Date: Thu Oct 12 08:25:05 EDT 2006 # # Smoothing Part 1 # # --------------------------------------------------- # ---------------------------------------------- # Load required libraries library(lattice) library(stats) # --------------------------------------------------- # 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 pdf("Smooth1TimeSeries1.pdf", height=5, width=6) xyplot(y ~ x, data = tSeries, panel = function(x, y) { panel.xyplot(x, y, type="l") panel.abline(lm(y ~ x, data = tSeries)) }, aspect=.7, xlab = list(label="Elapsed Time (days)",cex=1.75), ylab=list(label="Score",cex=1.75), scales=list(cex=1.25)) dev.off() # --------------------------------------------------- # Plot the timeseries using the plot command pdf("Smooth1TimeSeries2.pdf", height=5, width=6) plot(c(0,10), c(-4.5, 4.5), xlab="Elapsed Time (days)", ylab="Score", type="n") lines(tX, tY, type="l", cex=.5) dev.off() # --------------------------------------------------- # Calculate a running average windowWidth <- 2 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 pdf("Smooth1Mean1.pdf", height=5, width=6) plot(c(0,10), c(-4.5, 4.5), xlab="Elapsed Time (days)", ylab="Score", type="n") lines(tX, tY, type="p", cex=.5) lines(theCenter, theMeanSmooth, type="l", col=5, lwd=2) dev.off() # --------------------------------------------------- # 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 pdf("Smooth1Mean2.pdf", height=5, width=6) plot(c(0,10), c(-4.5, 4.5), xlab="Elapsed Time (days)", ylab="Score", type="n") lines(tX, tY, type="p", cex=.5) lines(theCenter, theMeanSmooth, type="l", col=5, lwd=2) dev.off() # --------------------------------------------------- # Calculate a running average with a wider window windowWidth <- 30 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 the running average pdf("Smooth1Mean3.pdf", height=5, width=6) plot(c(0,10), c(-4.5, 4.5), xlab="Elapsed Time (days)", ylab="Score", type="n") lines(tX, tY, type="p", cex=.5) lines(theCenter, theMeanSmooth, type="l", col=5, lwd=2) dev.off() # --------------------------------------------------- # 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 pdf("Smooth1BoxFilter1.pdf", height=5, width=6) plot(c(0,10), c(-4.5, 4.5), xlab="Elapsed Time (days)", ylab="Score", type="n") lines(tX, tY, type="p", cex=.5) lines(theCenter, theMeanSmooth, type="l", col=5, lwd=2) dev.off() # --------------------------------------------------- # 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 pdf("Smooth1BoxFilter2.pdf", height=5, width=6) 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) dev.off() # --------------------------------------------------- # 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] } theMeanSmooth4 <- rep(NA, theLength) theCenter4 <- rep(NA, theLength) for (i in (1+halfWidth):(theLength-halfWidth)) { theMeanSmooth4[i] <- sum(theMeanSmooth3[(i-halfWidth):(i+halfWidth)] * theKernel) theCenter4[i] <- theCenter3[i] } theMeanSmooth5 <- rep(NA, theLength) theCenter5 <- rep(NA, theLength) for (i in (1+halfWidth):(theLength-halfWidth)) { theMeanSmooth5[i] <- sum(theMeanSmooth4[(i-halfWidth):(i+halfWidth)] * theKernel) theCenter5[i] <- theCenter4[i] } # --------------------------------------------------- # Plot the timeseries and the smooths pdf("Smooth1BoxFilter2a.pdf", height=5, width=6) plot(c(0,200), c(-1.25, 1.25), xlab="Samples", ylab="Score", type="n") lines(tStepX, tStepY, type="p", cex=.5) dev.off() # --------------------------------------------------- # Plot the timeseries and the smooths pdf("Smooth1BoxFilter2b.pdf", height=5, width=6) 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) dev.off() # --------------------------------------------------- # Plot the timeseries and the smooths pdf("Smooth1BoxFilter2c.pdf", height=5, width=6) 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) dev.off() # --------------------------------------------------- # Plot the timeseries and the smooths pdf("Smooth1BoxFilter2d.pdf", height=5, width=6) 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) dev.off() # --------------------------------------------------- # Plot the timeseries and the smooths pdf("Smooth1BoxFilter2e.pdf", height=5, width=6) 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) lines(theCenter4, theMeanSmooth4, type="l", col=4, lwd=2) dev.off() # --------------------------------------------------- # Plot the timeseries and the smooths pdf("Smooth1BoxFilter2f.pdf", height=5, width=6) 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) lines(theCenter4, theMeanSmooth4, type="l", col=4, lwd=2) lines(theCenter5, theMeanSmooth5, type="l", col=3, lwd=2) dev.off() # --------------------------------------------------- # 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] } 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] } theMeanSmooth4 <- rep(NA, theLength) theCenter4 <- rep(NA, theLength) for (i in (1+halfWidth):(theLength-halfWidth)) { theMeanSmooth4[i] <- sum(theMeanSmooth3[(i-halfWidth):(i+halfWidth)] * theKernel) theCenter4[i] <- theCenter3[i] } theMeanSmooth5 <- rep(NA, theLength) theCenter5 <- rep(NA, theLength) for (i in (1+halfWidth):(theLength-halfWidth)) { theMeanSmooth5[i] <- sum(theMeanSmooth4[(i-halfWidth):(i+halfWidth)] * theKernel) theCenter5[i] <- theCenter4[i] } # --------------------------------------------------- # Plot the timeseries and smooth pdf("Smooth1SpencerFilter1a.pdf", height=5, width=6) plot(c(0,200), c(-1.25, 1.25), xlab="Samples", ylab="Score", type="n") lines(tStepX, tStepY, type="p", cex=.5) dev.off() pdf("Smooth1SpencerFilter1b.pdf", height=5, width=6) 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) dev.off() pdf("Smooth1SpencerFilter1c.pdf", height=5, width=6) 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) dev.off() pdf("Smooth1SpencerFilter1d.pdf", height=5, width=6) 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) dev.off() pdf("Smooth1SpencerFilter1e.pdf", height=5, width=6) 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) lines(theCenter4, theMeanSmooth4, type="l", col=4, lwd=2) dev.off() pdf("Smooth1SpencerFilter1f.pdf", height=5, width=6) 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) lines(theCenter4, theMeanSmooth4, type="l", col=4, lwd=2) lines(theCenter5, theMeanSmooth5, type="l", col=3, lwd=2) dev.off() # --------------------------------------------------- # 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] } 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] } theMeanSmooth4 <- rep(NA, theLength) theCenter4 <- rep(NA, theLength) for (i in (1+halfWidth):(theLength-halfWidth)) { theMeanSmooth4[i] <- sum(theMeanSmooth3[(i-halfWidth):(i+halfWidth)] * theKernel) theCenter4[i] <- theCenter3[i] } theMeanSmooth5 <- rep(NA, theLength) theCenter5 <- rep(NA, theLength) for (i in (1+halfWidth):(theLength-halfWidth)) { theMeanSmooth5[i] <- sum(theMeanSmooth4[(i-halfWidth):(i+halfWidth)] * theKernel) theCenter5[i] <- theCenter4[i] } # --------------------------------------------------- # Plot the timeseries and smooth pdf("Smooth1SpencerFilter2a.pdf", height=5, width=6) plot(c(0,10), c(-4.5, 4.5), xlab="Elapsed Time (days)", ylab="Score", type="n") lines(tX, tY, type="p", cex=.5) dev.off() pdf("Smooth1SpencerFilter2b.pdf", height=5, width=6) plot(c(0,10), c(-4.5, 4.5), xlab="Elapsed Time (days)", ylab="Score", type="n") lines(tX, tY, type="p", cex=.5) lines(theCenter, theMeanSmooth, type="l", col=5, lwd=2) dev.off() pdf("Smooth1SpencerFilter2c.pdf", height=5, width=6) plot(c(0,10), c(-4.5, 4.5), xlab="Elapsed Time (days)", ylab="Score", type="n") lines(tX, tY, type="p", cex=.5) lines(theCenter, theMeanSmooth, type="l", col=5, lwd=2) lines(theCenter2, theMeanSmooth2, type="l", col=6, lwd=2) dev.off() pdf("Smooth1SpencerFilter2d.pdf", height=5, width=6) plot(c(0,10), c(-4.5, 4.5), xlab="Elapsed Time (days)", ylab="Score", type="n") lines(tX, tY, 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) dev.off() pdf("Smooth1SpencerFilter2e.pdf", height=5, width=6) plot(c(0,10), c(-4.5, 4.5), xlab="Elapsed Time (days)", ylab="Score", type="n") lines(tX, tY, 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) lines(theCenter4, theMeanSmooth4, type="l", col=4, lwd=2) dev.off() pdf("Smooth1SpencerFilter2f.pdf", height=5, width=6) plot(c(0,10), c(-4.5, 4.5), xlab="Elapsed Time (days)", ylab="Score", type="n") lines(tX, tY, 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) lines(theCenter4, theMeanSmooth4, type="l", col=4, lwd=2) lines(theCenter5, theMeanSmooth5, type="l", col=3, lwd=2) dev.off() # --------------------------------------------------- # 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 pdf("Smooth1Median1.pdf", height=5, width=6) plot(c(0,10), c(-4.5, 4.5), xlab="Elapsed Time (days)", ylab="Score", type="n") lines(tX, tY, type="p", cex=.5) lines(theCenter, theMedianSmooth, type="l", col=5, lwd=2) dev.off() # --------------------------------------------------- # Calculate a running median using smooth() theSmooth <- smooth(tY, twice=F) # --------------------------------------------------- # Plot the timeseries and the running median pdf("Smooth1Median2.pdf", height=5, width=6) plot(c(0,10), c(-4.5, 4.5), xlab="Elapsed Time (days)", ylab="Score", type="n") lines(tX, tY, type="p", cex=.5) lines(tX, theSmooth, type="l", col=5, lwd=2) dev.off() # --------------------------------------------------- # Calculate a running median using smooth() theSmooth <- smooth(tY, twice=T) # --------------------------------------------------- # Plot the timeseries and the smooth pdf("Smooth1Median3.pdf", height=5, width=6) plot(c(0,10), c(-4.5, 4.5), xlab="Elapsed Time (days)", ylab="Score", type="n") lines(tX, tY, type="p", cex=.5) lines(tX, theSmooth, type="l", col=5, lwd=2) dev.off() # --------------------------------------------------- # Calculate a smoothing spline theSmooth <- smooth.spline(tX, tY, df=5)$y # --------------------------------------------------- # Plot the timeseries and the smooth pdf("Smooth1Spline1.pdf", height=5, width=6) plot(c(0,10), c(-4.5, 4.5), xlab="Elapsed Time (days)", ylab="Score", type="n") lines(tX, tY, type="p", cex=.5) lines(tX, theSmooth, type="l", col=5, lwd=2) dev.off() # --------------------------------------------------- # Plot the timeseries and the Spline for several different DF pdf("Smooth1Spline2a.pdf", height=5, width=6) plot(c(0,10), c(-4.5, 4.5), xlab="Elapsed Time (days)", ylab="Score", type="n") lines(tX, tY, type="p", cex=.5) dev.off() pdf("Smooth1Spline2b.pdf", height=5, width=6) plot(c(0,10), c(-4.5, 4.5), xlab="Elapsed Time (days)", ylab="Score", type="n") lines(tX, tY, type="p", cex=.5) lines(tX, smooth.spline(tX, tY, df=2)$y, type="l", col=5, lwd=2) dev.off() pdf("Smooth1Spline2c.pdf", height=5, width=6) plot(c(0,10), c(-4.5, 4.5), xlab="Elapsed Time (days)", ylab="Score", type="n") lines(tX, tY, type="p", cex=.5) lines(tX, smooth.spline(tX, tY, df=2)$y, type="l", col=5, lwd=2) lines(tX, smooth.spline(tX, tY, df=4)$y, type="l", col=6, lwd=2) dev.off() pdf("Smooth1Spline2d.pdf", height=5, width=6) plot(c(0,10), c(-4.5, 4.5), xlab="Elapsed Time (days)", ylab="Score", type="n") lines(tX, tY, type="p", cex=.5) lines(tX, smooth.spline(tX, tY, df=2)$y, type="l", col=5, lwd=2) lines(tX, smooth.spline(tX, tY, df=4)$y, type="l", col=6, lwd=2) lines(tX, smooth.spline(tX, tY, df=8)$y, type="l", col=7, lwd=2) dev.off() pdf("Smooth1Spline2e.pdf", height=5, width=6) plot(c(0,10), c(-4.5, 4.5), xlab="Elapsed Time (days)", ylab="Score", type="n") lines(tX, tY, type="p", cex=.5) lines(tX, smooth.spline(tX, tY, df=2)$y, type="l", col=5, lwd=2) lines(tX, smooth.spline(tX, tY, df=4)$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=16)$y, type="l", col=4, lwd=2) dev.off() pdf("Smooth1Spline2f.pdf", height=5, width=6) plot(c(0,10), c(-4.5, 4.5), xlab="Elapsed Time (days)", ylab="Score", type="n") lines(tX, tY, type="p", cex=.5) lines(tX, smooth.spline(tX, tY, df=2)$y, type="l", col=5, lwd=2) lines(tX, smooth.spline(tX, tY, df=4)$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=16)$y, type="l", col=4, lwd=2) lines(tX, smooth.spline(tX, tY, df=32)$y, type="l", col=3, lwd=2) dev.off() # --------------------------------------------------- # Quit the program