# -----------------------------------------------------
# Program: ThreeD1.S  
#  Author: Steven M. Boker
#    Date: Thu Nov 4 13:49:16 EST 2004
#
#  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 two histograms of the X variable

graphsheet(height=6.4,width=7.5)
hist(tX, nclass=10, probability=T,
        xlab = "X Score", 
        ylab = "Probability")
 
graphsheet(height=6.4,width=7.5)
hist(tX, nclass=5, probability=T,
        xlab = "X Score", 
        ylab = "Probability")
 
# -----------------------------------------------------
#  Calculate a smoothed version using a box kernel of width 10

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) 
}


graphsheet(height=6.4,width=7.5)
plot(c(-20,40),c(0,0.05), 
        xlab = "X Score", 
        ylab = "Probability", type="n")
lines(smoothX, smoothF, type="l")

# -----------------------------------------------------
#  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)

graphsheet(height=6.4,width=7.5)
plot(c(-20,40),c(0,0.05), 
        xlab = "X Score", 
        ylab = "Probability", type="n")
lines(smoothF2, type="l", col=6)

# -----------------------------------------------------
#  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)

graphsheet(height=6.4,width=7.5)
plot(c(-20,40),c(0,0.07), 
        xlab = "X Score", 
        ylab = "Probability", type="n")
lines(smoothF2, type="l", col=6)
lines(smoothF3, type="l", col=7)
text(30,.045,"rectangular", col=6)
text(30,.040,"gaussian", col=7)


# -----------------------------------------------------
#  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

graphsheet(height=6.4,width=7.5)
contour(smoothX, smoothY, smoothF,
        xlab = "X Score", 
        ylab = "Y Score")

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

graphsheet(height=6.4,width=7.5)
image(smoothX, smoothY, smoothF, 
        xlab = "X Score", 
        ylab = "Y Score")
contour(smoothX, smoothY, smoothF, add=T)

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

graphsheet(height=6.4,width=7.5)
persp(smoothX, smoothY, smoothF, 
        xlab = "X Score", 
        ylab = "Y Score",
	zlab = "Probability")


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

q()