# --------------------------------------------------------------------- # Program: Smoothing2.S # Author: Steven M. Boker # Date: Tue Oct 24 08:05:23 EDT 2006 # # Smoothing Part 2 # # --------------------------------------------------------------------- # --------------------------------------------------------------------- # 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) # --------------------------------------------------------------------- # Calculate a smoothing spline theSmooth <- smooth.spline(tX, tY, df=5)$y # --------------------------------------------------------------------- # Plot the data and the smooth pdf("Smooth2Spline1.pdf", height=5, width=6) 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) dev.off() # --------------------------------------------------------------------- # Plot the timeseries and the Spline for several different DF pdf("Smooth2Spline2.pdf", height=5, width=6) 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) lines(tX, sin(tX), type="l", col=1, lwd=1) dev.off() # --------------------------------------------------------------------- # Demonstrate the effect of the bicubic weighting function. tX <- c(1:100) q <-100 x <- 50 tD <- abs(tX - x) tXwin <- tX[rank(tD) <= q] tDX <- max(abs(tXwin - x)) tWt <- (1 - abs((tXwin-x)/tDX)^3)^3 pdf("Smooth2Bicubic1.pdf", height=5, width=6) plot(c(0,100), c(0, 1), xlab="X", ylab="Weight", type="n") lines(tXwin, tWt, type="p", cex=.5, col=6) lines(c(x,x), c(0,1), type="l", lty=2, col=3) dev.off() q <-100 x <- 100 tD <- abs(tX - x) tXwin <- tX[rank(tD) <= q] tDX <- max(abs(tXwin - x)) tWt <- (1 - abs((tXwin-x)/tDX)^3)^3 pdf("Smooth2Bicubic2.pdf", height=5, width=6) plot(c(0,100), c(0, 1), xlab="X", ylab="Weight", type="n") lines(tXwin, tWt, type="p", cex=.5, col=6) lines(c(x,x), c(0,1), type="l", lty=2, col=3) dev.off() q <- 50 for (x in seq(0,100,by=10)) { tD <- abs(tX - x) tXwin <- tX[rank(tD) <= q] tDX <- max(abs(tXwin - x)) tWt <- (1 - abs((tXwin-x)/tDX)^3)^3 tFilename <- paste("Smooth2BicubicX", x, ".pdf", sep="") pdf(tFilename, height=5, width=6) plot(c(0,100), c(0, 1), xlab="X", ylab="Weight", type="n") lines(tXwin, tWt, type="p", cex=.5, col=3) lines(c(x,x), c(0,1), type="l", lty=2, col=3) dev.off() } # --------------------------------------------------------------------- # Create a bivariate distribution with a piecewise linear relationship tX <- rnorm(200, mean=0, sd=10) tY <- rep(NA, 200) tY[tX <= 0] <- -10 + .2 * tX[tX <= 0] + rnorm(200,mean=0,sd=5)[tX <= 0] tY[tX > 0] <- -10 + 1.2 * tX[tX > 0] + rnorm(200,mean=0,sd=5)[tX > 0] tData <- data.frame(x=tX,y=tY) summary(tData) # --------------------------------------------------------------------- # Plot the simulated data pdf("Smooth2Scatter1.pdf", height=5, width=5) plot(c(-30, 30), c(-30, 30), xlab="X Score", ylab="Y Score", type="n") lines(tX, tY, type="p", cex=.5) dev.off() # --------------------------------------------------------------------- # Calculate a Lowess smooth theSmooth <- lowess(tX, tY) # --------------------------------------------------------------------- # Plot the data and the smooth pdf("Smooth2Lowess1.pdf", height=5, width=5) plot(c(-30, 30), c(-30, 30), xlab="X Score", ylab="Y Score", type="n") lines(tX, tY, type="p", cex=.5) lines(theSmooth, type="l", col=5, lwd=2) dev.off() # --------------------------------------------------------------------- # Drop the middle third of the data and recalculate a Lowess smooth tSelect <- tX < -10 | tX > 10 theSmooth <- lowess(tX[tSelect], tY[tSelect]) # --------------------------------------------------------------------- # Plot the data and the smooth pdf("Smooth2Lowess2.pdf", height=5, width=5) plot(c(-30, 30), c(-30, 30), xlab="X Score", ylab="Y Score", type="n") lines(tX[tSelect], tY[tSelect], type="p", cex=.5) lines(theSmooth, type="l", col=5, lwd=2) dev.off() # --------------------------------------------------------------------- # Calculate a Loess smooth tLoess <- loess(tY ~ tX) fitX <- seq(min(tX), max(tX), 1) fitY <- predict(tLoess, data.frame(tX = fitX), se = TRUE)$fit # --------------------------------------------------------------------- # Plot the data and the smooth pdf("Smooth2Lowess1c.pdf", height=5, width=5) plot(c(-30, 30), c(-30, 30), xlab="X Score", ylab="Y Score", type="n") lines(tX, tY, type="p", cex=.5) lines(fitX, fitY, type="l", col=5, lwd=2) dev.off() # --------------------------------------------------------------------- # Create a bivariate distribution with a cyclic relationship tX <- rnorm(200, mean=10, sd=3) tY <- 1.5*sin(tX) + rnorm(200,mean=0,sd=2) tData <- data.frame(x=tX,y=tY) summary(tData) # --------------------------------------------------------------------- # Plot the simulated data pdf("Smooth2Scatter2.pdf", height=5, width=5) plot(c(0, 20), c(-10, 10), xlab="Time in days", ylab="Y Score", type="n") lines(tX, tY, type="p", cex=.5) dev.off() # --------------------------------------------------------------------- # Calculate a Lowess smooth theSmooth <- lowess(tX, tY, f=1) # --------------------------------------------------------------------- # Plot the timeseries and the smooth pdf("Smooth2Lowess3.pdf", height=5, width=5) plot(c(0, 20), c(-10, 10), xlab="X Score", ylab="Y Score", type="n") lines(tX, tY, type="p", cex=.5) lines(theSmooth, type="l", col=5, lwd=2) dev.off() # --------------------------------------------------------------------- # Calculate a Loess smooth tLoess <- loess(tY ~ tX, span=1/2) fitX <- seq(min(tX), max(tX), 1) fitY <- predict(tLoess, data.frame(tX = fitX), se = TRUE)$fit # --------------------------------------------------------------------- # Plot the data and the smooth pdf("Smooth2Lowess3c.pdf", height=5, width=5) plot(c(0, 20), c(-10, 10), xlab="X Score", ylab="Y Score", type="n") lines(tX, tY, type="p", cex=.5) lines(fitX, fitY, type="l", col=5, lwd=2) dev.off() # --------------------------------------------------------------------- # Quit the program