mvstat.net MVKSA Book GPS Trajectory Contact
Daily
Temp World
Bank Earth
Quake Stem
Cell Foetal
Cardio Grevillea
  Dumb
Bell Air
Quality

World Bank development indicators for 2011. 
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)

## load data
data(worldbank)

## Fig 1.2
## scatter plot matrix
wb <- worldbank[,-1]
wb[,2] <- wb[,2]/1000
var.lab <- c(expression(CO[2]*"/capita"), "GDP/capita", "GDP growth", "Inflation", "Internet", "Agriculture")

pairs(wb, upper.panel=NULL, pch=16, col=pt.col(4,alpha=0.4), labels=var.lab, gap=0, oma=c(1.5,1.5,1.5,1.5)+0.1, mgp=c(0.5,0.5,0), mar=c(0,0,0,0))

## Fig 4.1a
## GDP-inflation scatter plot
xlab <- "GDP growth (%)"
ylab <- "Inflation (%)"
xlim <- c(-70 ,25)
ylim <- c(-25, 70)
zlim <- c(0, 0.022)
wb1 <- as.matrix(na.omit(worldbank[,4:5]))

## kernel density estimate
fhat.wb1 <- kde(x=wb1)

plot(fhat.wb1, cont=0, asp=1, xlim=xlim, ylim=ylim, xaxs="i", yaxs="i", xlab=xlab, ylab=ylab)
points(wb1, cex=1.4, pch=16, col=pt.col(4,alpha=0.5))
plotmixt(mus=c(-60,-10), Sigmas=fhat.wb1$H, props=1, drawlabels=FALSE, add=TRUE)

## Fig 4.1b
## fixed bandwidth kernel density estimate
plot(fhat.wb1, display="persp", thin=5, zlim=zlim, xlab=xlab, ylab=ylab, zlab="", box=TRUE, phi=20, col.fun=function(n){seq.col(n, alpha=1)}, border=grey(0, alpha=0.3))

## Fig 4.1c
## bandwidth varies with estimation point
fhat.wb1.ball <- kde.balloon(x=wb1, verbose=TRUE)
ep.ind1 <- c(11,31,51,71,91,111,131)
ep.ind2 <- c(11,31,51,71,91,111,131)
ep.sub <- as.matrix(expand.grid( fhat.wb1.ball$eval.points[[1]][ep.ind1], fhat.wb1.ball$eval.points[[2]][ep.ind2]))
H.pi.sub <- fhat.wb1.ball$H [as.matrix(expand.grid(ep.ind1, ep.ind2))]

plot(fhat.wb1, cont=0, asp=1, xaxs="i", yaxs="i", xlim=xlim, ylim=ylim, xlab=xlab, ylab=ylab)
for (i in 1:length(H.pi.sub))
{
    if (!is.na(H.pi.sub[i]) & H.pi.sub[i]>0)
    {
        plotmixt(mus=ep.sub[i,], Sigmas=H.pi.sub[i]*diag(2), props=1, drawlabels=FALSE, add=TRUE)
         points(ep.sub[i,1], ep.sub[i,2], cex=1, pch=21, col=grey(0,alpha=0.5), bg=pt.col(4))
    }
}

## Fig 4.1d
## balloon kernel density estimate
plot(fhat.wb1.ball, display="persp", thin=5, zlim=zlim, xlab=xlab, ylab=ylab, zlab="", box=TRUE, phi=20, col.fun=function(n){seq.col(n, alpha=1)}, border=grey(0, alpha=0.3))

## Fig 4.1e
## bandwidth varies with data point
fhat.wb1.sp <- kde.sp(x=wb1)
wb1.ind <- c(15, 17, 53, 62, 65 ,92 ,96, 107, 121, 129, 175)
H.ab.sub <- fhat.wb1.sp$H[wb1.ind]

plot(fhat.wb1, cont=0, asp=1, xaxs="i", yaxs="i", xlim=xlim, ylim=ylim, xlab=xlab, ylab=ylab)
for (i in 1:length(wb1.ind))
{
    if (H.ab.sub[i]>40) plotmixt(mus=wb1[wb1.ind[i],], Sigmas=H.ab.sub[i]*diag(2), props=1, drawlabels=FALSE, add=TRUE, xlim=c(-40,30), ylim=c(-20,80))
    else plotmixt(mus=wb1[wb1.ind[i],], Sigmas=H.ab.sub[i]*diag(2), props=1, drawlabels=FALSE, add=TRUE)
}
points(wb1[wb1.ind,], cex=1, pch=21, col=grey(0,alpha=0.5), bg=pt.col(4))

## Fig 4.1f
## sample point kernel density estimate
par(mar=c(0,0,0,0))
plot(fhat.wb1.sp, display="persp", thin=5, zlim=zlim, xlab=xlab, ylab=ylab, zlab="", box=TRUE, phi=20, col.fun=function(n){seq.col(n, alpha=1)}, border=grey(0, alpha=0.3))

## Fig 4.2a
## CO2-GDP scatter plot 
xlab <- expression(CO[2]*" emissions per capita ('000 Kg)")
ylab <- "GDP per capita (current '000 USD)"
xlim <- c(-5,55)
ylim <- c(-20,165)
wb2 <- as.matrix(na.omit(worldbank[,2:3]))
wb2[,2] <- wb2[,2]/1000
fhat.wb2 <- kde(x=wb2)

plot(fhat.wb2, cont=0, xaxs="i", yaxs="i", xlim=xlim, ylim=ylim, xlab=xlab, ylab=ylab)
points(wb2, cex=1.4, pch=16, col=pt.col(4,alpha=0.5))
rect(0, 0, 200, 200, border=grey(0.5, alpha=0.5), lty=2)

## Fig 4.2b
## fixed kernel density estimate
plot(fhat.wb2, drawlabels=FALSE, display="filled.contour", lwd=1, cont=seq(10,90,by=10), col.fun=function(n, alpha=1, ...){c("transparent", rev(sequential_hcl(n=n-1, h=290, power=0.6, alpha=alpha, ...)))}, xaxs="i", yaxs="i", xlim=xlim, ylim=ylim, xlab=xlab, ylab=ylab)
rect(0, 0, 200, 200, border=grey(0.5, alpha=0.5), lty=2)

## Fig 4.2c
## transformation kernel density estimate
fhat.wb2.trans <- kde(x=wb2, adj.positive=c(0,0), positive=TRUE)

plot(fhat.wb2.trans, drawlabels=FALSE, display="filled.contour", lwd=1, cont=seq(10,90,by=10), col.fun=function(n, alpha=1, ...){c("transparent", rev(sequential_hcl(n=n-1, h=290, power=0.6, alpha=alpha, ...)))}, xaxs="i", yaxs="i", xlim=xlim, ylim=ylim, xlab=xlab, ylab=ylab)
rect(0, 0, 200, 200, border=grey(0.5, alpha=0.5), lty=2)

## Fig 4.3a 
## Internet-agriculture scatter plot
xlab <- "Internet users (per 100)"
ylab <- "Agriculture value added (% of GDP)"
xlim <- c(-10, 110)
ylim <- c(-10, 110)
wb3 <- as.matrix(na.omit(worldbank[,6:7]))
fhat.wb3 <- kde(x=wb3)

plot(fhat.wb3, cont=0, asp=1, xaxs="i", yaxs="i", xlim=xlim, ylim=ylim, xlab=xlab, ylab=ylab)
points(wb3, cex=1.4, pch=16, col=pt.col(4,alpha=0.5))
rect(0, 0, 100, 100, border=grey(0.5, alpha=0.5), lty=2)
plotmixt(mus=c(80, 80), Sigmas=fhat.wb3$H, props=1, xlim=c(0,100), ylim=c(0,100), drawlabels=FALSE, add=TRUE)	  

## Fig 4.3b
## fixed kernel density estimate
plot(fhat.wb3, display="filled.contour", col.fun=function(n){seq.col(n, alpha=1)}, border=1, lwd=1, drawlabels=FALSE, xaxs="i", yaxs="i", cont=seq(10,90,by=20), asp=1, xaxs="i", yaxs="i", xlim=xlim, ylim=ylim, xlab=xlab, ylab=ylab)
rect(0, 0, 100, 100, border=grey(0.5, alpha=0.5), lty=2)

## Fig 4.3c
## beta boundary kernels 
fhat.wb3.beta <- kde.boundary(x=wb3, xmin=c(0,0), xmax=c(100,100), boundary.kernel="beta")

bb.sub <- c(1:5,35,150)
for (i in bb.sub)
{
    bb1 <- ks:::dmvbeta.prod.kernel2.2d(x=wb3[i,], xmax=c(100,100), hs=2*sqrt(diag(fhat.wb3.beta$H)))
    plot(bb1, display="slice", lwd=1, drawlabels=FALSE, add=i>1, asp=1, xaxs="i", yaxs="i", xlim=xlim, ylim=ylim, xlab=xlab, ylab=ylab)
}
points(wb3[bb.sub,], cex=1, pch=21, col=grey(0,alpha=0.5), bg=pt.col(4))
rect(0, 0, 100, 100, border=grey(0.5, alpha=0.5), lty=2)

## Fig 4.3d
## beta boundary kernel density estimate
plot(fhat.wb3.beta, display="filled.contour", col.fun=function(n){seq.col(n, alpha=1)}, border=1, lwd=1, drawlabels=FALSE, xaxs="i", yaxs="i", cont=seq(10,90,by=20), asp=1, xaxs="i", yaxs="i", xlim=xlim, ylim=ylim, xlab=xlab, ylab=ylab)
rect(0, 0, 100, 100, border=grey(0.5, alpha=0.5), lty=2)

## Fig 4.3e
## linear boundary kernels
fhat.wb3.LB <- kde.boundary(x=wb3, xmin=c(0,0), xmax=c(100,100), boundary.kernel="linear")

bb.sub <- c(1:5,35,150)
for (i in bb.sub)
{
    bb1 <- ks:::dmvnorm.LB.kernel.2d(x=wb3[i,], xmax=c(100,100), H=fhat.wb3.LB$H)
    plot(bb1, display="slice", lwd=1, drawlabels=FALSE, add=i>1, asp=1, xaxs="i", yaxs="i", xlim=xlim, ylim=ylim, xlab=xlab, ylab=ylab)
}
points(wb3[bb.sub,], cex=1, pch=21, col=grey(0,alpha=0.5), bg=pt.col(4))
rect(0, 0, 100, 100, border=grey(0.5, alpha=0.5), lty=2)

## Fig 4.3f
## linear boundary kernel density estimate
plot(fhat.wb3.LB, display="filled.contour", col.fun=function(n){seq.col(n, alpha=1)}, border=1, lwd=1, drawlabels=FALSE, cont=seq(10,90,by=20), asp=1, xaxs="i", yaxs="i", xlim=xlim, ylim=ylim, xlab=xlab, ylab=ylab)
rect(0, 0, 100, 100, border=grey(0.5, alpha=0.5), lty=2)