# -----------------------------------------------------
# Program: ThreeD1.S  
#  Author: Steven M. Boker
#    Date: Tue Nov 7 10:35:32 EST 2006
#
#  Probability Density Estimation, Density Plots
#  Three-D Plots Part 1
#
# -----------------------------------------------------

# -----------------------------------------------------
#  Create a mixed normal distribution and regression relation

tX <- c(rnorm(150, mean=10, sd=10), rnorm(50, mean=0, sd=5))
tY <- -10 + .5 * tX + rnorm(200,mean=0,sd=5)

tData <- data.frame(x=tX,y=tY)
summary(tData)


# -----------------------------------------------------
#  Plot histograms of the X variable

pdf("ThreeD1Hist1.pdf", height=5, width=6)
hist(tX, nclass=10, probability=T, main="", col=4,
        xlab = "X Score", 
        ylab = "Probability")
dev.off()

pdf("ThreeD1Hist2.pdf", height=5, width=6)
hist(tX, nclass=5, probability=T, main="", col=4,
        xlab = "X Score", 
        ylab = "Probability")
dev.off()
 
pdf("ThreeD1Hist3.pdf", height=5, width=6)
hist(tX, nclass=3, probability=T, main="", col=4,
        xlab = "X Score", 
        ylab = "Probability")
dev.off()
 
# -----------------------------------------------------
#  Calculate a smoothed version using a box kernel

smoothX <- seq(-20, 40, by=1)
smoothF <- rep(NA, length(smoothX))
h <- 10
for (i in 1:length(smoothX)) {
    x <- smoothX[i]
    t1 <- abs(x - tX)/h
    t2 <- rep(0, length(t1))
    t2[t1 < 1] <- 1/2
    smoothF[i] <- sum(t2)/(length(tX)*h) 
}


pdf("ThreeD1PDF1.pdf", height=5, width=6)
plot(c(-20,40),c(0,0.05), 
        xlab = "X Score", 
        ylab = "Probability", type="n")
lines(smoothX, smoothF, type="l")
dev.off()

smoothX <- seq(-20, 40, by=1)
smoothF <- rep(NA, length(smoothX))
h <- 5
for (i in 1:length(smoothX)) {
    x <- smoothX[i]
    t1 <- abs(x - tX)/h
    t2 <- rep(0, length(t1))
    t2[t1 < 1] <- 1/2
    smoothF[i] <- sum(t2)/(length(tX)*h) 
}


pdf("ThreeD1PDF1b.pdf", height=5, width=6)
plot(c(-20,40),c(0,0.05), 
        xlab = "X Score", 
        ylab = "Probability", type="n")
lines(smoothX, smoothF, type="l")
dev.off()

smoothX <- seq(-20, 40, by=1)
smoothF <- rep(NA, length(smoothX))
h <- 3
for (i in 1:length(smoothX)) {
    x <- smoothX[i]
    t1 <- abs(x - tX)/h
    t2 <- rep(0, length(t1))
    t2[t1 < 1] <- 1/2
    smoothF[i] <- sum(t2)/(length(tX)*h) 
}


pdf("ThreeD1PDF1c.pdf", height=5, width=6)
plot(c(-20,40),c(0,0.05), 
        xlab = "X Score", 
        ylab = "Probability", type="n")
lines(smoothX, smoothF, type="l")
dev.off()

# -----------------------------------------------------
#  Calculate a smoothed version using the density function
#  and a Rectangular kernel

smoothF2 <- density(tX, n=60, width=5, 
                    window="rectangular",
                    from=-20, to=40)

pdf("ThreeD1PDF2.pdf", height=5, width=6)
plot(c(-20,40),c(0,0.05), 
        xlab = "X Score", 
        ylab = "Probability", type="n")
lines(smoothF2, type="l", col=4, lwd=2)
dev.off()

# -----------------------------------------------------
#  Calculate a smoothed version using the density function
#  and a Gaussian and Rectangular kernel

smoothF2 <- density(tX, n=60, width=5, window="rectangular",
                    from=-20, to=40)

smoothF3 <- density(tX, n=60, width=5, window="gaussian",
                    from=-20, to=40)

pdf("ThreeD1PDF3.pdf", height=5, width=6)
plot(c(-20,40),c(0,0.07), 
        xlab = "X Score", 
        ylab = "Probability", type="n")
lines(smoothF2, type="l", col=2, lwd=2)
lines(smoothF3, type="l", col=4, lwd=2)
text(28,.055,"rectangular", col=2, cex=1.5)
text(28,.050,"gaussian", col=4, cex=1.5)
dev.off()


# -----------------------------------------------------
#  Calculate a smoothed bivariate density function using a
#  Rectangular kernel

smoothX <- seq(-20, 40, by=2)
smoothY <- seq(-20, 10, by=2)
smoothF <- matrix(NA, length(smoothX), length(smoothY))
h <- 5
for (i in 1:length(smoothX)) {
   x <- smoothX[i]
   t1x <- abs(x - tX)/h
   for (j in 1:length(smoothY)) {
      y <- smoothY[j]
      t1y <- abs(y - tY)/h
      t2 <- rep(0, length(t1x))
      t2[t1x < 1 & t1y < 1] <- 1/2
      smoothF[i,j] <- sum(t2)/(length(tX)*h) 
   }
}
smoothF <- smoothF/sum(smoothF)

# -----------------------------------------------------
#  Plot a line contour graph

pdf("ThreeD1Contour1.pdf", height=5, width=6)
contour(smoothX, smoothY, smoothF,
        xlab = "X Score", 
        ylab = "Y Score")
dev.off()

# -----------------------------------------------------
#  Plot an image density graph and add a line contour graph

pdf("ThreeD1Contour2.pdf", height=5, width=6)
image(smoothX, smoothY, smoothF, col = heat.colors(12),
        xlab = "X Score", 
        ylab = "Y Score")
contour(smoothX, smoothY, smoothF, add=T)
dev.off()

pdf("ThreeD1Contour2b.pdf", height=5, width=6)
image(smoothX, smoothY, smoothF, col = topo.colors(12),
        xlab = "X Score", 
        ylab = "Y Score")
contour(smoothX, smoothY, smoothF, add=T)
dev.off()

# -----------------------------------------------------
#  Plot a 3-d Perspective plot of the 2-D density function

pdf("ThreeD1Persp1.pdf", height=5, width=6)
persp(smoothX, smoothY, smoothF, 
      theta = 45, 
      phi = 25, 
      r = 3, 
      d = 1,
      shade = .25,
      col = 4,
      xlab = "X Score", 
      ylab = "Y Score",
      zlab = "Probability")
dev.off()

for(i in seq(0,360,by=20)) {
    tfilename <- paste("ThreeD1Persp2-", i, ".pdf", sep="")
    pdf(tfilename, height=5, width=6)
    persp(smoothX, smoothY, smoothF, 
          theta = i, 
          phi = 25, 
          r = 3, 
          d = 1,
          shade = .25,
          col = 4,
          xlab = "X Score", 
          ylab = "Y Score",
          zlab = "Probability")
    dev.off()
}

# -----------------------------------------------------
# Quit the program