library(ks)
library(colorspace)
library(RColorBrewer)

## colour functions
pt.col <- function(pos=1, alpha=1, ...) {return(paste0(brewer.pal(12, "Paired"), format(as.hexmode(round(alpha*255,0)), width=2))[pos])}
seq.col <- function(n, alpha=1, ...){c("transparent", rev(sequential_hcl(n=n-1, h=290, power=1.1, alpha=alpha, ...)))}

## plot parameters
par(oma=c(0,0,0,0)+0.1, mgp=c(1.8,0.5,0), mar=c(2.9,2.9,0,0)+0.1, cex.axis=1.2, cex.lab=1.2)
xlab <- expression(CO[2]*" g/"*m^3)
ylab <- expression(PM[10]*" (ppm)")
xlim <- c(350, 1050)
ylim <- c(30, 340)

## prepare data
data(air)
air <- air[, c("date", "time", "co2", "pm10")]
air2 <- reshape(air, idvar="date", timevar="time", direction="wide")
air <- as.matrix(na.omit(air2[,c("co2.20:00", "pm10.20:00")]))

## Fig 7.6a
## scatter plot
fhat.air <- kde(x=air)

plot(fhat.air, cont=0, xaxs="i", yaxs="i", xlim=xlim, ylim=ylim, xlab=xlab, ylab=ylab)
points(air, cex=1, pch=16, col=pt.col(4,alpha=0.3))

## Fig 7.6b
## classical kernel density estimate
plot(fhat.air, drawlabels=FALSE, display="filled.contour", lwd=1, col.fun=seq.col, cont=seq(20,80,by=20), xaxs="i", yaxs="i", xlim=xlim, ylim=ylim, xlab=xlab, ylab=ylab)

## Fig 7.6c
## deconvolution kernel density estimate
Sigma.air <- diag(c(var(air2[,"co2.19:00"] - air2["co2.21:00"], na.rm=TRUE), var(air2[,"pm10.19:00"] - air2[,"pm10.21:00"], na.rm=TRUE)))
fhat.air.dec <- kdcde(x=air, Sigma=Sigma.air, reg=0.00021)   ## kdcde is alias for dckde

plot(fhat.air.dec, drawlabels=FALSE, display="filled.contour", lwd=1, col.fun=seq.col, cont=seq(20,80,by=20), xaxs="i", yaxs="i", xlim=xlim, ylim=ylim, xlab=xlab, ylab=ylab)