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, ...)))}
seq.col2 <- function(n, alpha=1, ...) {c("transparent",tail(rev(diverge.col(2*n-1, alpha=alpha, ...)), n-1))}
diverge.col <- function(n, alpha=1, ...){cols <- paste0(brewer.pal(n, "PuOr"), as.hexmode(round(alpha*255,0))); cols[n %/% 2+1] <- grey(1, alpha=0); return(cols)}
## plot parameters
par(oma=c(0,0,0,0), mgp=c(0,0,0), mar=c(0,0,0,0), cex.axis=1.2, cex.lab=1.2)
xlim <- c(-5.6,5.6)
ylim <- c(-5.6,5.6)
## dumbbell density random sample
mus <- rbind(c(-2,2), c(0,0), c(2,-2))
Sigmas <- rbind(diag(2), 0.8*invvech(c(1, -0.9, 1)), diag(2))
cwt <- 3/11
props <- c((1-cwt)/2, cwt, (1-cwt)/2)
samp <- 1000
set.seed(81928192)
x <- rmvnorm.mixt(n=samp, mus=mus, Sigmas=Sigmas, props=props)
## Fig 2.10a ## target dumbbell density target.db <- plotmixt(mus=mus, Sigmas=Sigmas, props=props, drawlabels=FALSE, xlim=xlim, ylim=ylim, draw=FALSE) plot(target.db, drawlabels=FALSE, axes=FALSE, lwd=1, display="filled.contour", col.fun=seq.col, asp=1, xaxs="i", yaxs="i", xlim=c(-4,4), ylim=c(-4,4), xlab="", ylab="") box()
## Fig 2.10b ## kernel density estimate with unconstrained bandwidth fhat.db <- kde(x=x, xmin=c(xlim[1], ylim[1]), xmax=c(xlim[2], ylim[2])) db.single <- plotmixt(mus=c(3,3), Sigmas=fhat.db$H, props=1, drawlabels=FALSE, xlim=xlim, ylim=ylim, draw=FALSE) plot(fhat.db, drawlabels=FALSE, axes=FALSE, lwd=1, display="filled.contour", col.fun=seq.col, asp=1, xaxs="i", yaxs="i", xlim=c(-4,4), ylim=c(-4,4), xlab="", ylab="") plot(db.single, add=TRUE, drawlabels=FALSE, lwd=1) box()
## Fig 2.10c ## kernel density estimate with diagonal bandwidth fhat.db.diag <- kde(x=x, H=Hpi.diag(x), xmin=c(xlim[1], ylim[1]), xmax=c(xlim[2], ylim[2])) db.single.diag <- plotmixt(mus=c(3,3), Sigmas=fhat.db.diag$H, props=1, drawlabels=FALSE, xlim=xlim, ylim=ylim, draw=FALSE) plot(fhat.db.diag, drawlabels=FALSE, axes=FALSE, lwd=1, display="filled.contour", col.fun=seq.col, xaxs="i", yaxs="i", xlim=c(-4,4), ylim=c(-4,4), xlab="", ylab="") plot(db.single.diag, add=TRUE, drawlabels=FALSE, lwd=1) box()
## Fig 5.4a ## target dumbbell density gradient target1.db <- plotmixt(mus=mus, Sigmas=Sigmas, props=props, drawlabels=FALSE, xlim=xlim, ylim=ylim, draw=FALSE, deriv.order=1) plot(target1.db, display="quiver", thin=6, lwd=2, transf=1/12, scale=0.8, asp=1, xaxs="i", yaxs="i", xlim=c(-4,4), ylim=c(-4,4), xlab="", ylab="") box()
## Fig 5.4b ## kernel density gradient estimate with unconstrained bandwidth fhat1.db <- kdde(x=x, deriv.order=1, xmin=c(xlim[1], ylim[1]), xmax=c(xlim[2], ylim[2])) plot(fhat1.db, display="quiver", thin=6, lwd=2, transf=1/12, scale=0.8, asp=1, xaxs="i", yaxs="i", xlim=c(-4,4), ylim=c(-4,4), xlab="", ylab="") box()
## Fig 5.4c ## kernel density gradient estimate with diagonal bandwidth fhat1.db.diag <- kdde(x=x, H=Hpi.diag(x, deriv.order=1), binned=TRUE, deriv.order=1, xmin=c(xlim[1], ylim[1]), xmax=c(xlim[2], ylim[2])) png.create(file="fig5-4c-dumbbell-kdde.png") plot(fhat1.db.diag, display="quiver", thin=6, lwd=2, transf=1/12, scale=0.8, asp=1, xaxs="i", yaxs="i", xlim=c(-4,4), ylim=c(-4,4), xlab="", ylab="") box()
## Fig 5.4d ## target dumbbell density Hessian target2.db <- plotmixt(mus=mus, Sigmas=Sigmas, props=props, drawlabels=FALSE, xlim=xlim, ylim=ylim, draw=FALSE, deriv.order=2) target2.db.curv <- kcurv(target2.db) png.create(file="fig5-4d-dumbbell-kdde.png") plot(target2.db.curv, display="filled.contour", drawlabels=FALSE, lwd=1, col.fun=seq.col2, asp=1, xaxs="i", yaxs="i", xlim=c(-4,4), ylim=c(-4,4), xlab="", ylab="") box()
## Fig 5.4e ## kernel density Hessian estimate with unconstrained bandwidth fhat2.db <- kdde(x=x, deriv.order=2, xmin=c(xlim[1], ylim[1]), xmax=c(xlim[2], ylim[2])) fhat2.db.curv <- kcurv(fhat2.db) plot(fhat2.db.curv, display="filled.contour", drawlabels=FALSE, lwd=1, col.fun=seq.col2, asp=1, xaxs="i", yaxs="i", xlim=c(-4,4), ylim=c(-4,4), xlab="", ylab="") box()
## Fig 5.4f ## kernel density Hessian estimate with diagonal bandwidth fhat2.db.diag <- kdde(x=x, H=Hpi.diag(x, deriv.order=2), binned=TRUE, deriv.order=2, xmin=c(xlim[1], ylim[1]), xmax=c(xlim[2], ylim[2])) fhat2.db.curv.diag <- kcurv(fhat2.db.diag) plot(fhat2.db.curv.diag, display="filled.contour", drawlabels=FALSE, lwd=1, col.fun=seq.col2, asp=1, xaxs="i", yaxs="i", xlim=c(-4,4), ylim=c(-4,4), xlab="", ylab="") box()
