MVKSA | Daily temperature
mvstat.net
Home
  Mvksa   Daily Temp  
World Bank   Earth Quake   Stem Cell  
Foetal Cardio   Grevillea   Dumbbell  
Air Quality  
 
Daily minimum and maximum temperatures from the weather station in Badajoz, Spain from 1 January 1955 to 31 December 2015.
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))}
seq.col3d <- function(n, alpha=1, ...){rev(sequential_hcl(n=n, h=290, power=1.1, c.=c(80,70), l=c(30,50), alpha=alpha, ...))}
heat.col <- function(n, alpha=1, ...){heat_hcl(n, h=c(0, -100), l=c(75, 40), c=c(40, 80), power=1, alpha=alpha, ...)}
terrain.col <- function(n, alpha=1, ...){terrain_hcl(n, c=c(65, 0), l=c(45, 95), power=c(1/3, 1.5), alpha=alpha, ...)}
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)}

## prepare data
data(tempb)
tempb <- tempb[, c("tmin", "tmax")]

## 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 <- "Min temperature (deg C)" 
ylab <- "Max temperature (deg C)" 
xlim <- c(-5.5, 27.5)
ylim <- c(5.5, 42.5)
## Figure 1.1
## scattter plot   
plot(tempb, cex=0.7, pch=16, col=pt.col(4, alpha=0.1), asp=1, xaxs="i", yaxs="i", xlim=xlim, ylim=ylim, xlab=xlab, ylab=ylab) 

## Fig 2.1a
## histogram 
hist.tempb <- histde(tempb)

plot(hist.tempb, col.fun=seq.col, col.rect=grey(0,alpha=0.5), lty.rect=1, nbreaks=5, asp=1, xaxs="i", yaxs="i", xlim=xlim, ylim=ylim, xlab=xlab, ylab=ylab)
points(hist.tempb$eval.points[[1]][4], hist.tempb$eval.points[[2]][5], col=1, pch=16, cex=1.6)

## Fig 2.1b
## effect of anchor point placement
hist.tempb2 <- histde(tempb, adj=0.5)

plot(hist.tempb2, col.fun=seq.col, col.rect=grey(0,alpha=0.5), lty.rect=1, nbreaks=5, asp=1, xaxs="i", yaxs="i", xlim=xlim, ylim=ylim, xlab=xlab, ylab=ylab)
points(hist.tempb2$eval.points[[1]][4], hist.tempb2$eval.points[[2]][5], col=1, pch=16, cex=1.6)

## Fig 2.2a
## effect of binwidth - undersmoothing
hist.tempb.us <- histde(tempb, binw=hist.tempb$binw/3)

plot(hist.tempb.us, col.fun=seq.col, col.rect=grey(0,alpha=0.5), lty.rect=1, nbreaks=5, asp=1, xaxs="i", yaxs="i", xlim=xlim, ylim=ylim, xlab=xlab, ylab=ylab)

## Fig 2.2b
## effect of binwidth - oversmoothing
hist.tempb.os <- histde(tempb, binw=hist.tempb$binw*3)

plot(hist.tempb.os, col.fun=seq.col, col.rect=grey(0,alpha=0.5), lty.rect=1, nbreaks=5, asp=1, xaxs="i", yaxs="i", xlim=xlim, ylim=ylim, xlab=xlab, ylab=ylab)

## Fig 2.3a
## kernel functions on tempb subset
tempb.sub <- tempb[c(1608,2344,3820,5083,8184,10434,13470,13632,14044,21234),]
nsub <- nrow(tempb.sub)
H <- 10*Hpi(tempb)

for (i in 1:nsub)
{
    fhat.kern1 <- plotmixt(mus=as.numeric(tempb.sub[i,]), Sigmas=H, props=1, draw=FALSE, asp=1, xlim=xlim, ylim=ylim)
    plot(fhat.kern1, add=i>1, drawlabels=FALSE, lwd=1, col=grey(0.1, alpha=0.9), cont=seq(1,nsub,by=3)/nsub*45, asp=1, xaxs="i", yaxs="i", xlim=xlim, ylim=ylim, xlab=xlab, ylab=ylab)
}
points(tempb.sub, pch=21, bg=pt.col(4))

## Fig 2.3b
## kernel density estimate on tempb subset
fhat.tempb.sub <- kde(x=tempb.sub, H=H, xmin=c(-10.231763, -2.377353), xmax=c(29.03176, 48.55315))

plot(fhat.tempb.sub, display="filled.contour", drawlabels=FALSE, lwd=1, col.fun=seq.col, asp=1, xaxs="i", yaxs="i", xlim=xlim, ylim=ylim, xlab=xlab, ylab=ylab)
points(tempb.sub, pch=21, bg=pt.col(4))

## Fig 2.4a
## kernel density estimate - filled contour plot
fhat.tempb <- kde(x=tempb)
plot(fhat.tempb, display="filled.contour", drawlabels=FALSE, lwd=1, col.fun=seq.col, asp=1, xaxs="i", yaxs="i", xlim=xlim, ylim=ylim, xlab=xlab, ylab=ylab)

## Fig 2.4b
## kernel density estimate - perspective plot
plot(fhat.tempb, display="persp", phi=20, theta=20, col.fun=function(n){seq.col(n, alpha=1)[-1]}, border=grey(0, alpha=0.2), thin=3, xlab="Min temp.", ylab="Max temp.", zlab="")

## Fig 2.5a
## filled contour plot - 9 probability contour levels
plot(fhat.tempb, display="filled.contour", drawlabels=FALSE, lwd=1, col.fun=seq.col, cont=seq(10,90,by=10), asp=1, xaxs="i", yaxs="i", xlim=xlim, ylim=ylim, xlab=xlab, ylab=ylab)

## Fig 2.5b
## filled contour plot - 9 linear contour levels
ct.lev <- seq(0.0009, 0.0059, length=9)
	  
plot(fhat.tempb, display="filled.contour", drawlabels=FALSE, lwd=1, abs.cont=ct.lev, col.fun=seq.col, asp=1, xaxs="i", yaxs="i", xlim=xlim, ylim=ylim, xlab=xlab, ylab=ylab)

## Fig 2.6a
## filled contour plot - 50% probability contour level
plot(fhat.tempb, display="filled.contour", drawlabels=FALSE, lwd=1, cont=50, col=seq.col(9)[c(1,6)], asp=1, xaxs="i", yaxs="i", xlim=xlim, ylim=ylim, xlab=xlab, ylab=ylab)

## Fig 2.6b
## filled contour plot - 95% probability contour level
plot(fhat.tempb, display="filled.contour", drawlabels=FALSE, lwd=1, cont=95, col=seq.col(9)[c(1,2)], asp=1, xaxs="i", yaxs="i", xlim=xlim, ylim=ylim, xlab=xlab, ylab=ylab)

## Fig 2.7a
## filled contour plot - heat colour scheme 
plot(fhat.tempb, display="filled.contour", drawlabels=FALSE, lwd=1, col.fun=heat.col, asp=1, xaxs="i", yaxs="i", xlim=xlim, ylim=ylim, xlab=xlab, ylab=ylab)

## Fig 2.7b
## filled contour plot - terrain colour scheme 
plot(fhat.tempb, display="filled.contour", drawlabels=FALSE, lwd=1 col.fun=terrain.col, asp=1, xaxs="i", yaxs="i", xlim=xlim, ylim=ylim, xlab=xlab, ylab=ylab)

## Fig 2.8a
## unconstrained bandwidth matrix
plot(fhat.tempb, display="filled.contour", drawlabels=FALSE, lwd=1, asp=1, col.fun=seq.col, xaxs="i", yaxs="i", xlim=xlim, ylim=ylim, xlab=xlab, ylab=ylab)
plotmixt(mus=c(26,38), Sigmas=fhat.tempb$H, props=1, drawlabels=FALSE, add=TRUE)

## Fig 2.8b
## diagonal bandwidth matrix
plot(fhat.tempb.diag, display="filled.contour", drawlabels=FALSE, lwd=1 col.fun=seq.col, asp=1, xaxs="i", yaxs="i", xlim=xlim, ylim=ylim, ylab=ylab)
plotmixt(mus=c(26,38), Sigmas=fhat.tempb.diag$H, props=1, drawlabels=FALSE, add=TRUE)

## Fig 5.1a
## kernel density gradient estimate 
fhat1.tempb <- kdde(x=tempb, deriv.order=1)

## contour plot - partial derivative f^(1,0)
plot(fhat1.tempb, display="filled.contour", drawlabels=FALSE, lwd=1, col.fun=diverge.col, which.deriv.ind=1, asp=1, xaxs="i", yaxs="i", xlim=xlim, ylim=ylim, xlab=xlab, ylab=ylab)
box()

## Fig 5.1b
## contour plot - partial derivative f^(0,1)
plot(fhat1.tempb, display="filled.contour", drawlabels=FALSE, lwd=1, col.fun=diverge.col, which.deriv.ind=2, asp=1, xaxs="i", yaxs="i", xlim=xlim, ylim=ylim, xlab=xlab, ylab=ylab)
box()

## Fig 5.1c
## quiver plot
plot(fhat1.tempb, display="quiver", transf=1/12, thin=6, lwd=2, col.fun=diverge.col, asp=1, xaxs="i", yaxs="i", xlim=xlim, ylim=ylim, xlab=xlab, ylab=ylab)
box()

## Fig 5.2a
## kernel density Hessian estimate
fhat2.tempb <- kdde(x=tempb, deriv.order=2)
                                                  
## partial derivative f^(2,0)
plot(fhat2.tempb, display="filled.contour", drawlabels=FALSE, lwd=1, col.fun=diverge.col, which.deriv.ind=1, asp=1, xaxs="i", yaxs="i", xlim=xlim, ylim=ylim, xlab=xlab, ylab=ylab)
box()

## Fig 5.2b
## partial derivative f^(1,1)
plot(fhat2.tempb, display="filled.contour", drawlabels=FALSE, lwd=1, col.fun=diverge.col, which.deriv.ind=2, asp=1, xaxs="i", yaxs="i", xlim=xlim, ylim=ylim, xlab=xlab, ylab=ylab)
box()

## Fig 5.2c
## partial derivative f^(0,2)
plot(fhat2.tempb, display="filled.contour", drawlabels=FALSE, lwd=1, col.fun=diverge.col, which.deriv.ind=4, asp=1, xaxs="i", yaxs="i", xlim=xlim, ylim=ylim, xlab=xlab, ylab=ylab)
box()

## Fig 5.2d
## kernel summary curvature estimate
fhat2.tempb.curv <- kcurv(fhat2.tempb)

plot(fhat2.tempb.curv, display="filled.contour", drawlabels=FALSE, lwd=1, col.fun=seq.col2, asp=1, xaxs="i", yaxs="i", xlim=xlim, ylim=ylim, xlab=xlab, ylab=ylab)
box()

## Fig 6.1a
## modal region - 25% contour of density estimate
plot(fhat.tempb, display="filled.contour", drawlabels=FALSE, lwd=1, asp=1, cont=25, col.fun=seq.col, xaxs="i", yaxs="i", xlim=xlim, ylim=ylim, xlab=xlab, ylab=ylab)

## Fig 6.1b
## modal region - 25% contour of summary curvature 
plot(fhat2.tempb.curv, display="filled.contour", drawlabels=FALSE, cont=25, col.fun=seq.col2, lwd=1, asp=1, xaxs="i", yaxs="i", xlim=xlim, ylim=ylim, xlab=xlab, ylab=ylab)
box()                                               

## Fig 6.3a
## Devroye-Wise support estimate
supp.tempb.dw <- dwsupp(tempb, H=fhat.tempb$H, xmin=c(-20,-10), xmax=c(35,55))
chull.tempb.dw <- ksupp(supp.tempb.dw, cont=99, convex.hull=TRUE)

plot(fhat.tempb, cont=0, asp=1, xaxs="i", yaxs="i", xlim=xlim, ylim=ylim, xlab=xlab, ylab=ylab)
points(tempb, cex=0.4, pch=16, col=pt.col(4,alpha=0.1))
plot(supp.tempb.dw, add=TRUE, drawlabels=FALSE, lwd=2, col=1, abs.cont=0.5, asp=1)
polygon(chull.tempb.dw, lwd=2, lty=3, border=grey(0,alpha=0.5))

## Fig 6.3b
## kernel support estimate
chull.tempb9995 <- ksupp(fhat.tempb, cont=99.95, convex.hull=TRUE)
chull.tempb <- fhat.tempb$x[chull(fhat.tempb$x),]

plot(fhat.tempb, cont=0, asp=1, xaxs="i", yaxs="i", xlim=xlim, ylim=ylim, xlab=xlab, ylab=ylab)
points(tempb, cex=0.4, pch=16, col=pt.col(4,alpha=0.1))
plot(fhat.tempb, add=TRUE, drawlabels=FALSE, lwd=2, col=1, abs.cont=contourLevels(fhat.tempb, cont=c(95, 99.5, 99.95)))
polygon(chull.tempb9995, lwd=2, lty=2, border=seq.col(3)[3])

## Fig 6.7
## kernel mean shift clustering
## warning: kms may take hours to execute on tempb
ms.tempb <- kms(tempb, verbose=TRUE, min.clust.size=0.05*nrow(tempb))

plot(fhat.tempb, cont=0, asp=1, xaxs="i", yaxs="i", xlim=xlim, ylim=ylim, ylab=ylab, xlab=xlab)
plot(ms.tempb, pch=16, col=pt.col(c(4,10,2), alpha=0.1), add=TRUE)
points(ms.tempb$mode, pch=21, col=1, bg=pt.col(7))

## Fig 6.11
## kernel feature significance
fs.tempb <- kfs(tempb)

plot(fs.tempb, display="filled.contour", col=pt.col(7), lwd=2, asp=1, xaxs="i", yaxs="i", xlim=xlim, ylim=ylim, xlab=xlab, ylab=ylab)
plot(fhat2.tempb.curv, cont=25, add=TRUE, drawlabels=FALSE, col=pt.col(10), lwd=5, lty="12")

## Fig 8.2a
## coarse estimation grid size
fhat.tempb.gs1 <- kde(x=tempb, bgridsize=rep(11,2))

plot(fhat.tempb, display="slice", drawlabels=FALSE, col=pt.col(10,0.5), lwd=3, asp=1, xaxs="i", yaxs="i", xlim=xlim, ylim=ylim, xlab=xlab, ylab=ylab)
plot(fhat.tempb.gs1, display="slice", drawlabels=FALSE, add=TRUE, col=1,lwd=3, lty=3)

## Fig 8.2b
## fine estimation grid 
fhat.tempb.gs2 <- kde(x=tempb, bgridsize=rep(301,2))

plot(fhat.tempb, display="slice", drawlabels=FALSE, col=pt.col(10,0.5), lwd=3, asp=1, xaxs="i", yaxs="i", xlim=xlim, ylim=ylim, xlab=xlab, ylab=ylab)
plot(fhat.tempb.gs2, display="slice", drawlabels=FALSE, add=TRUE, col=1,lwd=3, lty=3)