# --------------------------------------------------------------------- # Program: Smoothing2.S # Author: Steven M. Boker # Date: Tue Nov 6 09:29:06 EST 2007 # # 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) # --------------------------------------------------------------------- # Demonstrate the effect of the bicubic weighting function for LOWESS. 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 trueX <- c(-30:30) trueY <- -10 + .2 * trueX trueY[trueX>0] <- -10 + 1.2 * trueX[trueX>0] pdf("Smooth2Scatter0.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() # --------------------------------------------------------------------- # A linear model might make us happy tLM <- lm (y~x, data=tData) summary(tLM) 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) abline(tLM) dev.off() pdf("Smooth2Scatter2.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(trueX, trueY, type="l", col="red") abline(tLM) 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) lines(trueX, trueY, type="l", col="red") dev.off() # --------------------------------------------------------------------- # Plot the weights when there is a gap in the X values tXGap <- c(1:30, 70:100) q <- 50 for (x in seq(0,100,by=10)) { tD <- abs(tXGap - x) tXwin <- tXGap[rank(tD) <= q] tDX <- max(abs(tXwin - x)) tWt <- (1 - abs((tXwin-x)/tDX)^3)^3 tFilename <- paste("Smooth2BicubicXGap", 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="red") dev.off() } # --------------------------------------------------------------------- # Drop the middle third of the data tSelect <- tX < -10 | tX > 10 pdf("Smooth2ScatterGap0.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) dev.off() # --------------------------------------------------------------------- # A linear model might make us happy tLM <- lm (y[tSelect]~x[tSelect], data=tData) summary(tLM) pdf("Smooth2ScatterGap1.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) abline(tLM) dev.off() pdf("Smooth2ScatterGap2.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) abline(tLM) lines(trueX, trueY, type="l", col="red") dev.off() # --------------------------------------------------------------------- # Plot the data and the smooth theSmooth <- lowess(tX[tSelect], tY[tSelect]) 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) lines(trueX, trueY, type="l", col="red") 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) trueX <- seq(3,17,by=.1) trueY <- 1.5*sin(trueX) # --------------------------------------------------------------------- # Plot the simulated data pdf("Smooth2ScatterCycle2.pdf", height=5, width=6) 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() pdf("Smooth2ScatterCycle3.pdf", height=5, width=6) plot(c(0, 20), c(-10, 10), xlab="Time in days", ylab="Y Score", type="n") lines(tX, tY, type="p", cex=.5) lines(trueX, trueY, type="l", col="red", lwd=2) dev.off() # --------------------------------------------------------------------- # Calculate a Lowess smooth theSmooth <- lowess(tX, tY, f=1) # --------------------------------------------------------------------- # Plot the timeseries and the smooth pdf("Smooth2LowessCycle3.pdf", height=5, width=6) 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="blue", lwd=2) lines(trueX, trueY, type="l", col="red", lwd=2) dev.off() for (theF in c(1,3/4, 1/2, 1/4, 1/8)) { theSmooth <- lowess(tX, tY, f=theF) tFilename <- paste("Smooth2LowessCycle3F", theF, ".pdf", sep="") pdf(tFilename, height=5, width=6) 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="blue", lwd=2) lines(trueX, trueY, type="l", col="red", 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("Smooth2Loess3c.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="blue", lwd=2) lines(trueX, trueY, type="l", col="red", lwd=2) dev.off() # --------------------------------------------------------------------- # Quit the program