# ------------------------------------------------------ # Program: VectorFields2.R # Author: Steven M. Boker # Date: Tue Nov 27 08:39:41 EST 2007 # # Vector Field plotting part 2 # # ------------------------------------------------------ # ------------------------------------------------------ # Use Multidimensional vectors to visualize the edaStudent data source("edaStudent.sdd") smoothX <- seq(0, 14, by=2) smoothY <- seq(0.125, 1, by=.125) smoothP <- matrix(NA, length(smoothX), length(smoothY)) smoothZ <- matrix(NA, length(smoothX), length(smoothY)) smoothW <- matrix(NA, length(smoothX), length(smoothY)) smoothV <- matrix(NA, length(smoothX), length(smoothY)) hx <- 5/14 hy <- 5/1 for (i in 1:length(smoothX)) { x <- smoothX[i] t1x <- abs(x - (edaStudent\$income/10000))*hx for (j in 1:length(smoothY)) { y <- smoothY[j] t1y <- abs(y - edaStudent\$pposths)*hy t2 <- rep(0, length(t1x)) t2[t1x < 1 & t1y < 1 & !is.na(t1x) & !is.na(t1y)] <- 1/2 smoothP[i,j] <- sum(t2)/(length(edaStudent\$income)) smoothZ[i,j] <- mean(edaStudent\$ztest[t1x < 1 & t1y < 1 & !is.na(edaStudent\$ztest) & !is.na(t1x) & !is.na(t1y)]) smoothW[i,j] <- mean(edaStudent\$hscgpa[t1x < 1 & t1y < 1 & !is.na(edaStudent\$hscgpa) & !is.na(t1x) & !is.na(t1y)]) smoothV[i,j] <- mean(edaStudent\$totunits[t1x < 1 & t1y < 1 & !is.na(edaStudent\$totunits) & !is.na(t1x) & !is.na(t1y)]) } } smoothP <- 100 * (smoothP/sum(smoothP)) smoothStudentFrame <- data.frame( income=rep(smoothX,length(smoothY)), pposths=rep(smoothY,each=length(smoothX)), density=c(smoothP), ztest=(c(smoothZ) - mean(c(smoothZ), na.rm=TRUE)) / sqrt(var(c(smoothZ), na.rm=TRUE)), hscgpa=(c(smoothW) - mean(c(smoothW), na.rm=TRUE)) / sqrt(var(c(smoothW), na.rm=TRUE)), totunits=(c(smoothV) - mean(c(smoothV), na.rm=TRUE)) / sqrt(var(c(smoothV), na.rm=TRUE)) ) summary(smoothStudentFrame) # ------------------------------------------------------ # Plot a whisker field of the fitted slopes pdf("VF2StudentMVF.pdf", height=5, width=5) plot(c(0,14), c(0,1.0), type="n" , xlab="Income (in 10K)", ylab="Probability of Post HS") scale <- .2 yScale <- 1/7 for (i in 1:length(smoothStudentFrame\$income)) { tx <- smoothStudentFrame\$income[i] ty <- smoothStudentFrame\$pposths[i] td <- smoothStudentFrame\$density[i] ta <- smoothStudentFrame\$ztest[i] txinc <- (scale * td) / sqrt( 1 + (ta^2)) tyinc <- ta * txinc * yScale lines(c(tx), c(ty), type="p", cex=0.5) lines(c(tx, tx+txinc), c(ty, ty+tyinc), type="l") tx1 <- tx+txinc ty1 <- ty+tyinc ta <- smoothStudentFrame\$hscgpa[i] txinc <- (scale * td) / sqrt( 1 + (ta^2)) tyinc <- ta * txinc * yScale lines(c(tx1, tx1+txinc), c(ty1, ty1+tyinc), type="l") tx2 <- tx1+txinc ty2 <- ty1+tyinc ta <- smoothStudentFrame\$totunits[i] txinc <- (scale * td) / sqrt( 1 + (ta^2)) tyinc <- ta * txinc * yScale lines(c(tx2, tx2+txinc), c(ty2, ty2+tyinc), type="l") text(9,0.02, "o-Test-GPA-Units") } dev.off() # ------------------------------------------------------ # Use Multidimensional vectors to visualize the Iris data iris <- data.frame(read.table("iris.dat", header=T)) summary(iris) smoothX <- seq(4, 8, by=.5) smoothY <- seq(2, 4.5, by=.25) smoothP <- matrix(NA, length(smoothX), length(smoothY)) smoothZ <- matrix(NA, length(smoothX), length(smoothY)) smoothW <- matrix(NA, length(smoothX), length(smoothY)) hx <- 5/4 hy <- 5/2.5 for (i in 1:length(smoothX)) { x <- smoothX[i] t1x <- abs(x - (iris\$sepal.length))*hx for (j in 1:length(smoothY)) { y <- smoothY[j] t1y <- abs(y - iris\$sepal.width)*hy t2 <- rep(0, length(t1x)) t2[t1x < 1 & t1y < 1 & !is.na(t1x) & !is.na(t1y)] <- 1/2 smoothP[i,j] <- sum(t2)/(length(iris\$sepal.length)) smoothZ[i,j] <- mean(iris\$petal.length[t1x < 1 & t1y < 1 & !is.na(iris\$petal.length) & !is.na(t1x) & !is.na(t1y)]) smoothW[i,j] <- mean(iris\$petal.width[t1x < 1 & t1y < 1 & !is.na(iris\$petal.width) & !is.na(t1x) & !is.na(t1y)]) } } smoothP <- 100 * (smoothP/sum(smoothP)) smoothIrisFrame <- data.frame( sepal.length=rep(smoothX,length(smoothY)), sepal.width=rep(smoothY,each=length(smoothX)), density=c(smoothP), petal.length=(c(smoothZ) - mean(c(smoothZ), na.rm=TRUE)) / sqrt(var(c(smoothZ), na.rm=TRUE)), petal.width=(c(smoothW) - mean(c(smoothW), na.rm=TRUE)) / sqrt(var(c(smoothW), na.rm=TRUE))) summary(smoothIrisFrame) # ------------------------------------------------------ # Plot a whisker field of the fitted slopes pdf("VF2IrisMVF.pdf", height=5, width=5) plot(c(4, 8), c(1.75, 4.5), type="n" , xlab="Sepal Length", ylab="Sepal Width") scale <- .1 for (i in 1:length(smoothIrisFrame\$sepal.length)) { tx <- smoothIrisFrame\$sepal.length[i] ty <- smoothIrisFrame\$sepal.width[i] td <- smoothIrisFrame\$density[i] ta <- smoothIrisFrame\$petal.length[i] txinc <- (scale * td) / sqrt( 1 + (ta^2)) tyinc <- ta * txinc lines(c(tx), c(ty), type="p") lines(c(tx, tx+txinc), c(ty, ty+tyinc), type="l") tx1 <- tx+txinc ty1 <- ty+tyinc ta <- smoothIrisFrame\$petal.width[i] txinc <- (scale * td) / sqrt( 1 + (ta^2)) tyinc <- ta * txinc lines(c(tx1, tx1+txinc), c(ty1, ty1+tyinc), type="l") text(6,1.8, "o-PetalLength-PetalWidth") } dev.off() # ------------------------------------------------------ # Create and plot four coupled time series TimeSeries1 <- rep(NA, 201) TimeSeries1[1] <- 90 for (i in 2:201) { TimeSeries1[i] <- .98 * TimeSeries1[i-1] } TimeSeries1 <- 100 - TimeSeries1 + rnorm(201, mean=0, sd=10) TimeSeries2 <- .3 * TimeSeries1 * sin(c(1:201)/6) + rnorm(201, mean=0, sd=3) + TimeSeries1 TimeSeries3 <- TimeSeries1[seq(201,1,by=-1)] + rnorm(201, mean=0, sd=3) TimeSeries4 <- .3 * TimeSeries1 * sin(2+c(1:201)/6) + rnorm(201, mean=0, sd=3) + TimeSeries3 TimeSeriesAge <- c(0:200) / 3 pdf("VF2TSPlot1.pdf", height=5, width=6) plot(c(0,70), c(0,130), type="n" , xlab="Age", ylab="Total Score") lines(TimeSeriesAge, TimeSeries1, type="l", lty=1, col=1) lines(TimeSeriesAge, TimeSeries2, type="l", lty=1, col=2) lines(TimeSeriesAge, TimeSeries3, type="l", lty=1, col=4) lines(TimeSeriesAge, TimeSeries4, type="l", lty=1, col=5) dev.off() # ------------------------------------------------------ # Use Multidimensional vectors to visualize the coupled time series smoothX <- seq(min(TimeSeries1), max(TimeSeries1), by=10) smoothY <- seq(min(TimeSeries2), max(TimeSeries1), by=10) smoothP <- matrix(NA, length(smoothX), length(smoothY)) smoothZ <- matrix(NA, length(smoothX), length(smoothY)) smoothW <- matrix(NA, length(smoothX), length(smoothY)) smoothV <- matrix(NA, length(smoothX), length(smoothY)) hx <- 5/(max(TimeSeries1) - min(TimeSeries1)) hy <- 5/(max(TimeSeries2) - min(TimeSeries2)) for (i in 1:length(smoothX)) { x <- smoothX[i] t1x <- abs(x - (TimeSeriesAge))*hx for (j in 1:length(smoothY)) { y <- smoothY[j] t1y <- abs(y - TimeSeries1)*hy t2 <- rep(0, length(t1x)) t2[t1x < 1 & t1y < 1 & !is.na(t1x) & !is.na(t1y)] <- 1/2 smoothP[i,j] <- sum(t2)/(length(TimeSeriesAge)) smoothZ[i,j] <- mean(TimeSeries2[t1x < 1 & t1y < 1 & !is.na(TimeSeries2) & !is.na(t1x) & !is.na(t1y)]) smoothW[i,j] <- mean(TimeSeries3[t1x < 1 & t1y < 1 & !is.na(TimeSeries3) & !is.na(t1x) & !is.na(t1y)]) smoothV[i,j] <- mean(TimeSeries4[t1x < 1 & t1y < 1 & !is.na(TimeSeries4) & !is.na(t1x) & !is.na(t1y)]) } } smoothP <- 100 * (smoothP/sum(smoothP)) smoothTSFrame <- data.frame( age=rep(smoothX,length(smoothY)), x1=rep(smoothY,each=length(smoothX)), density=c(smoothP), x2=(c(smoothZ) - mean(c(smoothZ), na.rm=TRUE)) / sqrt(var(c(smoothZ), na.rm=TRUE)), x3=(c(smoothW) - mean(c(smoothW), na.rm=TRUE)) / sqrt(var(c(smoothW), na.rm=TRUE)), x4=(c(smoothV) - mean(c(smoothV), na.rm=TRUE)) / sqrt(var(c(smoothV), na.rm=TRUE))) summary(smoothTSFrame) # ------------------------------------------------------ # Plot a whisker field of the fitted slopes pdf("VF2TSWhiskerField1.pdf", height=5, width=5) plot(c(0, 70), c(0, 130), type="n" , xlab="Age", ylab="Time Series 1") scale <- 2 for (i in 1:length(smoothTSFrame\$age)) { tx <- smoothTSFrame\$age[i] ty <- smoothTSFrame\$x1[i] td <- smoothTSFrame\$density[i] ta <- smoothTSFrame\$x2[i] txinc <- (scale * td) / sqrt( 1 + (ta^2)) tyinc <- ta * txinc lines(c(tx), c(ty), type="p") lines(c(tx, tx+txinc), c(ty, ty+tyinc), type="l") tx1 <- tx+txinc ty1 <- ty+tyinc ta <- smoothTSFrame\$x3[i] txinc <- (scale * td) / sqrt( 1 + (ta^2)) tyinc <- ta * txinc lines(c(tx1, tx1+txinc), c(ty1, ty1+tyinc), type="l") tx2 <- tx1+txinc ty2 <- ty1+tyinc ta <- smoothTSFrame\$x4[i] txinc <- (scale * td) / sqrt( 1 + (ta^2)) tyinc <- ta * txinc lines(c(tx2, tx2+txinc), c(ty2, ty2+tyinc), type="l") text(35,125, "o-X2-X3-X4") } dev.off()