# Stochastik IV, Übungsblatt 2

# Aufgabe 5

library(MASS)

# Einlesen der Daten und Einbinden in Suchpfad

data.olives <- 
read.table("D:/Univ Augsburg/Lehre/SS07/VO_MultiVar/Datensaetze/Olives.txt",
head = TRUE, sep = "\t", quote = "")

attach(data.olives)

# Teilaufgabe (a), Scatterplot der Variablen palmitoleic und oleic 
# (eingefärbt nach den Regionen North, Sardinia und South)

plot(palmitoleic, oleic, main = "Scatterplot \"oleic versus palmitoleic\"",
col = as.numeric(Region), cex = 1.5, pch = 19)

legend(c(181, 263), c(7980, 8350), c("North", "Sardinia", "South"), 
       fill = c(1, 2, 3))

# Teilaufgabe (b), 2-dimensionale Dichteschätzer

# Bandbreite bw.nrd0 (Silverman rule)

par(mfrow = c(1,2))
plot(palmitoleic, oleic, col = "lightblue")
dens.est.nrd0 <- kde2d(palmitoleic, oleic, h = c(bw.nrd0(palmitoleic), 
                                                 bw.nrd0(oleic)), n = 30)
contour(dens.est.nrd0, nlevels = 5, method = "simple", add = TRUE)
title("Scatter and Contour Plots (bw.nrd0)")
persp(dens.est.nrd0, main = "Perspective Plot (bw.nrd0)",
xlab = "palmitoleic", ylab = "oleic", zlab = "density", 
phi = 30, theta = 0, d = 5, col = "lightblue")

# Bandbreite bw.nrd (Scott rule)

windows()

par(mfrow = c(1,2))
plot(palmitoleic, oleic, col = "lightblue")
dens.est.nrd <- kde2d(palmitoleic, oleic, h = c(bw.nrd(palmitoleic), 
                                                bw.nrd(oleic)), n = 30)
contour(dens.est.nrd, nlevels = 5, method = "simple", add = TRUE)
title("Scatter and Contour Plots (bw.nrd)")
persp(dens.est.nrd, main = "Perspective Plot (bw.nrd)",
xlab = "palmitoleic", ylab = "oleic", zlab = "density", 
phi = 30, theta = 0, d = 5, col = "lightblue")

# Bandbreite bw.ucv (unbiased cross-validation)

windows()

par(mfrow = c(1,2))
plot(palmitoleic, oleic, col = "lightblue")
dens.est.ucv <- kde2d(palmitoleic, oleic, h = c(bw.ucv(palmitoleic), 
                                                bw.ucv(oleic)), n = 30)
contour(dens.est.ucv, nlevels = 5, method = "simple", add = TRUE)
title("Scatter and Contour Plots (bw.ucv)")
persp(dens.est.ucv, main = "Perspective Plot (bw.ucv)",
xlab = "palmitoleic", ylab = "oleic", zlab = "density", 
phi = 30, theta = 0, d = 5, col = "lightblue")

# Bandbreite bw.bcv (biased cross-validation)

windows()

par(mfrow = c(1,2))
plot(palmitoleic, oleic, col = "lightblue")
dens.est.bcv <- kde2d(palmitoleic, oleic, h = c(bw.bcv(palmitoleic), 
                                                bw.bcv(oleic)), n = 30)
contour(dens.est.bcv, nlevels = 5, method = "simple", add = TRUE)
title("Scatter and Contour Plots (bw.bcv)")
persp(dens.est.bcv, main = "Perspective Plot (bw.bcv)",
xlab = "palmitoleic", ylab = "oleic", zlab = "density", 
phi = 30, theta = 0, d = 5, col = "lightblue")

# Bandbreite bw.SJ.ste (Sheather & Jones rule, ste)

windows()

par(mfrow = c(1,2))
plot(palmitoleic, oleic, col = "lightblue")
dens.est.SJ.ste <- kde2d(palmitoleic, oleic, 
                         h = c(bw.SJ(palmitoleic, method = "ste"), 
                               bw.SJ(oleic, method = "ste")), n = 30)
contour(dens.est.SJ.ste, nlevels = 5, method = "simple", add = TRUE)
title("Scatter and Contour Plots (bw.SJ.ste)")
persp(dens.est.SJ.ste, main = "Perspective Plot (bw.SJ.ste)",
xlab = "palmitoleic", ylab = "oleic", zlab = "density", 
phi = 30, theta = 0, d = 5, col = "lightblue")

# Bandbreite bw.SJ.dpi (Sheather & Jones rule, dpi)

windows()

par(mfrow = c(1,2))
plot(palmitoleic, oleic, col = "lightblue")
dens.est.SJ.dpi <- kde2d(palmitoleic, oleic, 
                         h = c(bw.SJ(palmitoleic, method = "dpi"), 
                               bw.SJ(oleic, method = "dpi")), n = 30)
contour(dens.est.SJ.dpi, nlevels = 5, method = "simple", add = TRUE)
title("Scatter and Contour Plots (bw.SJ.dpi)")
persp(dens.est.SJ.dpi, main = "Perspective Plot (bw.SJ.dpi)",
xlab = "palmitoleic", ylab = "oleic", zlab = "density", 
phi = 30, theta = 0, d = 5, col = "lightblue")

# Bandbreite bandwidth.nrd (MASS rule)

windows()

par(mfrow = c(1,2))
plot(palmitoleic, oleic, col = "lightblue")
dens.est.bandwidth.nrd <- kde2d(palmitoleic, oleic, 
                                h = c(bandwidth.nrd(palmitoleic), 
                                      bandwidth.nrd(oleic)), n = 30)
contour(dens.est.bandwidth.nrd, nlevels = 5, method = "simple", add = TRUE)
title("Scatter and Contour Plots (bandwidth.nrd)")
persp(dens.est.bandwidth.nrd, main = "Perspective Plot (bandwidth.nrd)",
xlab = "palmitoleic", ylab = "oleic", zlab = "density", 
phi = 30, theta = 0, d = 5, col = "lightblue")

# Bandbreite ucv (MASS rule)

windows()

par(mfrow = c(1,2))
plot(palmitoleic, oleic, col = "lightblue")
dens.est.ucv <- kde2d(palmitoleic, oleic, 
                      h = c(ucv(palmitoleic), 
                            ucv(oleic)), n = 30)
contour(dens.est.ucv, nlevels = 5, method = "simple", add = TRUE)
title("Scatter and Contour Plots (ucv)")
persp(dens.est.ucv, main = "Perspective Plot (ucv)",
xlab = "palmitoleic", ylab = "oleic", zlab = "density", 
phi = 30, theta = 0, d = 5, col = "lightblue")

# Bandbreite bcv (MASS rule)

windows()

par(mfrow = c(1,2))
plot(palmitoleic, oleic, col = "lightblue")
dens.est.bcv <- kde2d(palmitoleic, oleic, 
                      h = c(bcv(palmitoleic), 
                            bcv(oleic)), n = 30)
contour(dens.est.bcv, nlevels = 5, method = "simple", add = TRUE)
title("Scatter and Contour Plots (bcv)")
persp(dens.est.bcv, main = "Perspective Plot (bcv)",
xlab = "palmitoleic", ylab = "oleic", zlab = "density", 
phi = 30, theta = 0, d = 5, col = "lightblue")

# Bandbreite width.SJ-ste (MASS rule)

windows()

par(mfrow = c(1,2))
plot(palmitoleic, oleic, col = "lightblue")
dens.est.width.SJ.ste <- kde2d(palmitoleic, oleic, 
                               h = c(width.SJ(palmitoleic, method = "ste"), 
                                     width.SJ(oleic, method = "ste")), n = 30)
contour(dens.est.width.SJ.ste, nlevels = 5, method = "simple", add = TRUE)
title("Scatter and Contour Plots (width.SJ.ste)")
persp(dens.est.width.SJ.ste, main = "Perspective Plot (width.SJ.ste)",
xlab = "palmitoleic", ylab = "oleic", zlab = "density", 
phi = 30, theta = 0, d = 5, col = "lightblue")

# Bandbreite width.SJ-dpi (MASS rule)

windows()

par(mfrow = c(1,2))
plot(palmitoleic, oleic, col = "lightblue")
dens.est.width.SJ.dpi <- kde2d(palmitoleic, oleic, 
                               h = c(width.SJ(palmitoleic, method = "dpi"), 
                                     width.SJ(oleic, method = "dpi")), n = 30)
contour(dens.est.width.SJ.dpi, nlevels = 5, method = "simple", add = TRUE)
title("Scatter and Contour Plots (width.SJ.dpi)")
persp(dens.est.width.SJ.dpi, main = "Perspective Plot (width.SJ.dpi)",
xlab = "palmitoleic", ylab = "oleic", zlab = "density", 
phi = 30, theta = 0, d = 5, col = "lightblue")

# BEMERKUNGEN:
#
# 1. Anhand obiger Graphiken gesehen ist eine Bandbreitenwahl basierend auf 
# der MASS-Implementation bandwidth.nrd() (4 mal der Basisfunktion bw.nrd();
# siehe 2.a und 2.b unten) für eine sinnvolle Dichteschätzung empfehlenswert!
#
# 2. Bemerkung: bandwidth.nrd() ist evtl. speziell fuer den bivariaten 
# Dichteschätzer formuliert?!
#
# 2.a. Die Funktion bandwidth.nrd() ist definiert als
#
# function(x)
# {
# 	r <- quantile(x, c(0.25, 0.75))
#    	h <- (r[2] - r[1])/1.34
#   	4 * 1.06 * min(sqrt(var(x)), h) * length(x)^(-1/5)
# }
#
# 2.b. im Gegensatz zur Definition von bw.nrd()
#
# function (x) 
# {
#   	if (length(x) < 2) 
#     	stop("need at least 2 data points")
#    	r <- quantile(x, c(0.25, 0.75))
#    	h <- (r[2] - r[1])/1.34
#    	1.06 * min(sqrt(var(x)), h) * length(x)^(-1/5)
# }

# Teilaufgabe (c), Randdichten als Projektionen auf die Achsen

oldpar <- par(no.readonly = TRUE)
par(mar = c(5.1, 4.1, 6.1, 5.1))

# Scatter- und Contour-Plot

plot(palmitoleic, oleic, col = as.numeric(Region), cex = 1.25, pch = 19)
dens.est.bandwidth.nrd <- kde2d(palmitoleic, oleic, 
                                h = c(bandwidth.nrd(palmitoleic), 
                                      bandwidth.nrd(oleic)), n = 30)
contour(dens.est.bandwidth.nrd, nlevels = 5, method = "simple", lwd = 1.5,
        add = TRUE)

# 1D Rug- / Scatter-Plot

rug(palmitoleic, ticksize = 0.0225, col = "lightblue")
rug(oleic, side = 2, ticksize = 0.0225, col = "lightpink")
box()

# Projektion auf palmitoleic-Achse

dens.palm <- density(palmitoleic, bw = bandwidth.nrd(palmitoleic), from = 0)
dens.val.palm <- pretty(range(dens.palm$y)) 
axis(side = 4, at = 6500 + 400000 * dens.val.palm, lab = dens.val.palm) 
lines(dens.palm$x, 6500 + 400000 * dens.palm$y, lwd = 2, col = "lightblue")
mtext("density (palmitoleic)", side = 4, line = 3)

# Projektion auf oleic-Achse

dens.oleic <- density(oleic, bw = bandwidth.nrd(oleic), from = 0)
dens.val.oleic <- pretty(range(dens.oleic$y)) 
axis(side = 3, at = 50 + 350000 * dens.val.oleic, lab = dens.val.oleic) 
lines(50 + 350000 * dens.oleic$y, dens.oleic$x, lwd = 2, col = "lightpink")
mtext("density (oleic)", side = 3, line = 3)

legend(200, 8500, 
       c("North", "Sardinia", "South", "rug & density (palmitoleic)",
         "rug & density (oleic)"), 
       fill = c(1, 2, 3, "lightblue", "lightpink"), bty = "n")