Nothing
###############################################################################
## Multivariate kernel density derivative estimate
###############################################################################
kdde <- function(x, H, h, deriv.order=0, gridsize, gridtype, xmin, xmax, supp=3.7, eval.points, binned, bgridsize, positive=FALSE, adj.positive, w, deriv.vec=TRUE, verbose=FALSE)
{
## default values
r <- deriv.order
ksd <- ks.defaults(x=x, w=w, binned=binned, bgridsize=bgridsize, gridsize=gridsize)
d <- ksd$d; n <- ksd$n; w <- ksd$w
binned <- ksd$binned
gridsize <- ksd$gridsize
bgridsize <- ksd$bgridsize
## clip data to xmin,xmax grid for binned estimation
grid.clip <- binned
if (grid.clip)
{
if (!missing(xmax)) xmax <- xmax[1:d]
if (!missing(xmin)) xmin <- xmin[1:d]
if (positive & missing(xmin)) { xmin <- rep(0,d) }
xt <- truncate.grid(x=x, y=w, xmin=xmin, xmax=xmax)
x <- xt$x; w <- xt$y; n <- length(w)
}
## default bandwidths
if (d==1 & missing(h))
{
if (positive) x1 <- log(x) else x1 <- x
h <- hpi(x=x1, nstage=2, binned=default.bflag(d=d, n=n), deriv.order=r)
}
if (d>1 & missing(H))
{
if (positive) x1 <- log(x) else x1 <- x
if ((r>0) & (d>2)) nstage <- 1 else nstage <- 2
H <- Hpi(x=x1, nstage=nstage, binned=default.bflag(d=d, n=n), deriv.order=r, verbose=verbose)
}
## compute binned estimator
if (binned)
{
if (positive & is.vector(x))
{
y <- log(x)
fhat <- kdde.binned(x=y, H=H, h=h, deriv.order=r, bgridsize=bgridsize, xmin=xmin, xmax=xmax, w=w)
fhat$estimate <- fhat$estimate/exp(fhat$eval.points)
fhat$eval.points <- exp(fhat$eval.points)
fhat$x <- x
}
else
fhat <- kdde.binned(x=x, H=H, h=h, deriv.order=r, bgridsize=bgridsize, xmin=xmin, xmax=xmax, w=w, deriv.vec=deriv.vec, verbose=verbose)
if (!missing(eval.points))
{
fhat$estimate <- predict(fhat, x=eval.points)
fhat$eval.points <- eval.points
}
}
else
{
## compute exact (non-binned) estimator
## 1-dimensional
if (d==1)
{
if (!missing(H) & !missing(h)) stop("Both H and h are both specified")
if (missing(h)) h <- sqrt(H)
if (missing(eval.points))
fhat <- kdde.grid.1d(x=x, h=h, gridsize=gridsize, supp=supp, xmin=xmin, xmax=xmax, gridtype=gridtype, w=w, deriv.order=r)
else
fhat <- kdde.points.1d(x=x, h=h, eval.points=eval.points, w=w, deriv.order=r)
}
## multi-dimensional
else
{
if (is.data.frame(x)) x <- as.matrix(x)
if (missing(eval.points))
{
if (d==2)
fhat <- kdde.grid.2d(x=x, H=H, gridsize=gridsize, supp=supp, xmin=xmin, xmax=xmax, gridtype=gridtype, w=w, deriv.order=r, deriv.vec=deriv.vec, verbose=verbose)
else if (d==3)
fhat <- kdde.grid.3d(x=x, H=H, gridsize=gridsize, supp=supp, xmin=xmin, xmax=xmax, gridtype=gridtype, w=w, deriv.order=r, deriv.vec=deriv.vec, verbose=verbose)
else
stop("Need to specify eval.points for more than 3 dimensions")
}
else
fhat <- kdde.points(x=x, H=H, eval.points=eval.points, w=w, deriv.order=r, deriv.vec=deriv.vec)
}
}
fhat$binned <- binned
fhat$names <- parse.name(x)
fhat$type <- "kdde"
class(fhat) <- "kdde"
return(fhat)
}
###############################################################################
## Multivariate binned kernel density derivative estimate
###############################################################################
kdde.binned <- function(x, H, h, deriv.order, bgridsize, xmin, xmax, bin.par, w, deriv.vec=TRUE, deriv.index, verbose=FALSE)
{
r <- deriv.order
if (length(r)>1) stop("deriv.order should be a non-negative integer")
## linear binning
if (missing(bin.par))
{
if (is.vector(x)) { d <- 1 }
else { d <- ncol(x) }
if (d==1)
if (missing(H)) { H <- as.matrix(h^2) }
else { h <- sqrt(H); H <- as.matrix(H) }
if (d==1) Hd <- H else Hd <- diag(diag(H))
bin.par <- binning(x=x, H=Hd, h=h, bgridsize=bgridsize, xmin=xmin, xmax=xmax, supp=3.7+max(r), w=w)
}
else
{
if (!is.list(bin.par$eval.points)) { d <- 1; bgridsize <- length(bin.par$eval.points) }
else { d <- length(bin.par$eval.points); bgridsize <- sapply(bin.par$eval.points, length) }
w <- bin.par$w
if (d==1)
if (missing(H)) H <- as.matrix(h^2)
else { h <- sqrt(H); H <- as.matrix(H) }
}
if (d==1)
{
fhat <- kdde.binned.1d(h=h, deriv.order=r, bin.par=bin.par)
eval.points <- fhat$eval.points
est <- fhat$estimate
}
else
{
ind.mat <- dmvnorm.deriv(x=rep(0,d), mu=rep(0,d), Sigma=H, deriv.order=r, only.index=TRUE, deriv.vec=deriv.vec)
fhat.grid <- kdde.binned.nd(H=H, deriv.order=r, bin.par=bin.par, verbose=verbose, deriv.vec=deriv.vec)
}
if (missing(x)) x <- NULL
if (d==1)
{
if (r==0) fhat <- list(x=x, eval.points=unlist(eval.points), estimate=est, h=h, H=h^2, gridtype="linear", gridded=TRUE, binned=TRUE, names=NULL, w=w)
else
fhat <- list(x=x, eval.points=unlist(eval.points), estimate=est, h=h, H=h^2, gridtype="linear", gridded=TRUE, binned=TRUE, names=NULL, w=w, deriv.order=r, deriv.ind=r)
}
else
{
if (r==0)
fhat <- list(x=x, eval.points=fhat.grid$eval.points, estimate=fhat.grid$estimate[[1]], H=H, gridtype="linear", gridded=TRUE, binned=TRUE, names=NULL, w=w)
else
fhat <- list(x=x, eval.points=fhat.grid$eval.points, estimate=fhat.grid$estimate, H=H, gridtype="linear", gridded=TRUE, binned=TRUE, names=NULL, w=w, deriv.order=r, deriv.ind=ind.mat)
}
class(fhat) <- "kdde"
return(fhat)
}
kdde.binned.1d <- function(h, deriv.order, bin.par)
{
r <- deriv.order
n <- sum(bin.par$counts)
a <- min(bin.par$eval.points)
b <- max(bin.par$eval.points)
M <- length(bin.par$eval.points)
L <- min(ceiling((4+r)*h*(M-1)/(b-a)), M-1)
delta <- (b-a)/(M-1)
N <- 2*L-1
grid1 <- seq(-(L-1), L-1)
keval <- dnorm.deriv(x=delta*grid1, mu=0, sigma=h, deriv.order=r)/n
est <- symconv.1d(keval, bin.par$counts)
return(list(eval.points=bin.par$eval.points, estimate=est))
}
kdde.binned.nd <- function(H, deriv.order, bin.par, verbose=FALSE, deriv.vec=TRUE)
{
d <- ncol(H)
r <- deriv.order
n <- sum(bin.par$counts)
a <- sapply(bin.par$eval.points, min)
b <- sapply(bin.par$eval.points, max)
M <- sapply(bin.par$eval.points, length)
L <- pmin(ceiling((4+r)*max(sqrt(abs(diag(H))))*(M-1)/(b-a)), M-1)
delta <- (b-a)/(M-1)
if (min(L)<=0) warning(paste("Binning grid too coarse for current (small) bandwidth: consider increasing grid size for dimensions", toString(which(pmin(L)<=1))))
N <- 2*L-1
if(d==2)
{
grid1 <- seq(-(L[1]-1), L[1]-1)
grid2 <- seq(-(L[2]-1), L[2]-1)
xgrid <- expand.grid(delta[1]*grid1, delta[2]*grid2)
}
else if (d==3)
{
grid1 <- seq(-(L[1]-1), L[1]-1)
grid2 <- seq(-(L[2]-1), L[2]-1)
grid3 <- seq(-(L[3]-1), L[3]-1)
xgrid <- expand.grid(delta[1]*grid1, delta[2]*grid2, delta[3]*grid3)
}
else if (d==4)
{
grid1 <- seq(-(L[1]-1), L[1]-1)
grid2 <- seq(-(L[2]-1), L[2]-1)
grid3 <- seq(-(L[3]-1), L[3]-1)
grid4 <- seq(-(L[4]-1), L[4]-1)
xgrid <- expand.grid(delta[1]*grid1, delta[2]*grid2, delta[3]*grid3, delta[4]*grid4)
}
deriv.index <- dmvnorm.deriv(x=rep(0,d), mu=rep(0,d), Sigma=H, deriv.order=r, add.index=TRUE, only.index=TRUE, deriv.vec=TRUE)
deriv.index.minimal <- dmvnorm.deriv(x=rep(0,d), mu=rep(0,d), Sigma=H, deriv.order=r, add.index=TRUE, only.index=TRUE, deriv.vec=FALSE)
if (verbose) pb <- txtProgressBar()
if (r==0)
{
n.seq <- block.indices(1, nrow(xgrid), d=d, r=r, diff=FALSE)
est.list <- vector(1, mode="list")
est.list[[1]] <- array(0, dim=dim(bin.par$counts))
}
else if (r>0)
{
n.deriv <- nrow(deriv.index)
n.deriv.minimal <- nrow(deriv.index.minimal)
if (deriv.vec) n.est.list <- n.deriv else n.est.list <- n.deriv.minimal
est.list <- vector(n.est.list, mode="list")
for (j in 1:n.est.list) est.list[[j]] <- array(0, dim=dim(bin.par$counts))
if (d^r >= 3^7) n.seq <- block.indices(1, nrow(xgrid), d=d, r=r, diff=FALSE, block.limit=1e4)
else n.seq <- block.indices(1, nrow(xgrid), d=d, r=r, diff=FALSE, block.limit=1e5)
}
for (i in 1:(length(n.seq)-1))
{
if (verbose) setTxtProgressBar(pb, i/(length(n.seq)-1))
keval <- dmvnorm.deriv(x=xgrid[n.seq[i]:(n.seq[i+1]-1),], mu=rep(0,d), Sigma=H, deriv.order=r, add.index=TRUE, deriv.vec=FALSE)$deriv/n
if (r==0) keval <- as.matrix(keval, ncol=1)
else if (is.vector(keval)) keval <- as.matrix(t(keval), nrow=1)
est <- list()
## loop over only unique partial derivative indices
nderiv <- nrow(deriv.index.minimal)
if (!(is.null(nderiv)))
for (s in 1:nderiv)
{
if (deriv.vec) deriv.rep.index <- which.mat(deriv.index.minimal[s,], deriv.index)
else deriv.rep.index <- s
kevals <- array(keval[,s], dim=N)
est.temp <- symconv.nd(kevals, bin.par$counts, d=d)
for (s2 in 1:length(deriv.rep.index)) est[[deriv.rep.index[s2]]] <- est.temp
}
else
{
kevals <- lapply(as.data.frame(keval), function(x) { array(x, dim=N) })
est <- lapply(kevals, function(x) { symconv.nd(x, bin.par$counts, d=d) })
}
if (r==0) est.list[[1]] <- est.list[[1]] + est[[1]]
else if (r>0) for (j in 1:n.est.list) est.list[[j]] <- est.list[[j]] + est[[j]]
}
if (verbose) close(pb)
fhatr <- list(eval.points=bin.par$eval.points, estimate=est.list, deriv.order=r)
return(fhatr)
}
#############################################################################
## Univariate kernel density derivative estimate on a grid
#############################################################################
kdde.grid.1d <- function(x, h, gridsize, supp=3.7, positive=FALSE, adj.positive, xmin, xmax, gridtype, w, deriv.order=0)
{
r <- deriv.order
if (r==0)
fhatr <- kde(x=x, h=h, gridsize=gridsize, supp=supp, positive=positive, adj.positive=adj.positive, xmin=xmin, xmax=xmax, gridtype=gridtype, w=w)
else
{
if (missing(xmin)) xmin <- min(x) - h*supp
if (missing(xmax)) xmax <- max(x) + h*supp
if (missing(gridtype)) gridtype <- "linear"
y <- x
gridtype1 <- match.arg(gridtype, c("linear", "sqrt"))
if (gridtype1=="linear")
gridy <- seq(xmin, xmax, length=gridsize)
else if (gridtype1=="sqrt")
{
gridy.temp <- seq(sign(xmin)*sqrt(abs(xmin)), sign(xmax)*sqrt(abs(xmax)), length=gridsize)
gridy <- sign(gridy.temp) * gridy.temp^2
}
gridtype.vec <- gridtype1
n <- length(y)
est <- dnorm.deriv.mixt(x=gridy, mus=y, sigmas=rep(h, n), props=w/n, deriv.order=r)
fhatr <- list(x=y, eval.points=gridy, estimate=est, h=h, H=h^2, gridtype=gridtype.vec, gridded=TRUE, binned=FALSE, names=NULL, w=w, deriv.order=r, deriv.ind=deriv.order)
class(fhatr) <- "kdde"
}
return(fhatr)
}
##############################################################################
## Bivariate kernel density derivative estimate on a grid
## Computes all mixed partial derivatives for a given deriv.order
##############################################################################
kdde.grid.2d <- function(x, H, gridsize, supp, gridx=NULL, grid.pts=NULL, xmin, xmax, gridtype, w, deriv.order=0, deriv.vec=TRUE, verbose=FALSE)
{
d <- 2
r <- deriv.order
if (r==0)
fhatr <- kde(x=x, H=H, gridsize=gridsize, supp=supp, xmin=xmin, xmax=xmax, gridtype=gridtype, w=w, verbose=verbose)
else
{
## initialise grid
n <- nrow(x)
if (is.null(gridx))
gridx <- make.grid.ks(x, matrix.sqrt(H), tol=supp, gridsize=gridsize, xmin=xmin, xmax=xmax, gridtype=gridtype)
suppx <- make.supp(x, matrix.sqrt(H), tol=supp)
if (is.null(grid.pts))
grid.pts <- find.gridpts(gridx, suppx)
nderiv <- d^r
fhat.grid <- list()
for (k in 1:nderiv)
fhat.grid[[k]] <- matrix(0, nrow=length(gridx[[1]]), ncol=length(gridx[[2]]))
if (verbose) pb <- txtProgressBar()
for (i in 1:n)
{
## compute evaluation points
eval.x <- gridx[[1]][grid.pts$xmin[i,1]:grid.pts$xmax[i,1]]
eval.y <- gridx[[2]][grid.pts$xmin[i,2]:grid.pts$xmax[i,2]]
eval.x.ind <- c(grid.pts$xmin[i,1]:grid.pts$xmax[i,1])
eval.y.ind <- c(grid.pts$xmin[i,2]:grid.pts$xmax[i,2])
eval.x.len <- length(eval.x)
eval.pts <- expand.grid(eval.x, eval.y)
## Create list of matrices for different partial derivatives
fhat <- dmvnorm.deriv(x=eval.pts, mu=x[i,], Sigma=H, deriv.order=r)
## place vector of density estimate values `fhat' onto grid 'fhat.grid'
for (k in 1:nderiv)
for (j in 1:length(eval.y))
fhat.grid[[k]][eval.x.ind, eval.y.ind[j]] <- fhat.grid[[k]][eval.x.ind, eval.y.ind[j]] + w[i]*fhat[((j-1) * eval.x.len + 1):(j * eval.x.len),k]
if (verbose) setTxtProgressBar(pb, i/n)
}
if (verbose) close(pb)
for (k in 1:nderiv) fhat.grid[[k]] <- fhat.grid[[k]]/n
gridx1 <- list(gridx[[1]], gridx[[2]])
ind.mat <- dmvnorm.deriv(x=rep(0,d), mu=rep(0,d), Sigma=H, deriv.order=r, only.index=TRUE)
if (!deriv.vec)
{
fhat.grid.vech <- list()
deriv.ind <- unique(ind.mat)
for (i in 1:nrow(deriv.ind))
{
which.deriv <- which.mat(deriv.ind[i,], ind.mat)[1]
fhat.grid.vech[[i]] <- fhat.grid[[which.deriv]]
}
ind.mat <- deriv.ind
fhat.grid <- fhat.grid.vech
}
fhatr <- list(x=x, eval.points=gridx1, estimate=fhat.grid, H=H, gridtype=gridx$gridtype, gridded=TRUE, binned=FALSE, names=NULL, w=w, deriv.order=deriv.order, deriv.ind=ind.mat)
}
return(fhatr)
}
##############################################################################
## Trivariate kernel density derivative estimate on a grid
## Computes all mixed partial derivatives for a given deriv.order
##############################################################################
kdde.grid.3d <- function(x, H, gridsize, supp, gridx=NULL, grid.pts=NULL, xmin, xmax, gridtype, w, deriv.order=0, deriv.vec=TRUE, verbose=FALSE)
{
d <- 3
r <- deriv.order
if (r==0)
fhatr <- kde(x=x, H=H, gridsize=gridsize, supp=supp, xmin=xmin, xmax=xmax, gridtype=gridtype, w=w, verbose=verbose)
else
{
## initialise grid
n <- nrow(x)
if (is.null(gridx))
gridx <- make.grid.ks(x, matrix.sqrt(H), tol=supp, gridsize=gridsize, xmin=xmin, xmax=xmax, gridtype=gridtype)
suppx <- make.supp(x, matrix.sqrt(H), tol=supp)
if (is.null(grid.pts))
grid.pts <- find.gridpts(gridx, suppx)
nderiv <- d^r
fhat.grid <- list()
for (k in 1:nderiv)
fhat.grid[[k]] <- array(0, dim=gridsize)
if (verbose) pb <- txtProgressBar()
for (i in 1:n)
{
## compute evaluation points
eval.x <- gridx[[1]][grid.pts$xmin[i,1]:grid.pts$xmax[i,1]]
eval.y <- gridx[[2]][grid.pts$xmin[i,2]:grid.pts$xmax[i,2]]
eval.z <- gridx[[3]][grid.pts$xmin[i,3]:grid.pts$xmax[i,3]]
eval.x.ind <- c(grid.pts$xmin[i,1]:grid.pts$xmax[i,1])
eval.y.ind <- c(grid.pts$xmin[i,2]:grid.pts$xmax[i,2])
eval.z.ind <- c(grid.pts$xmin[i,3]:grid.pts$xmax[i,3])
eval.x.len <- length(eval.x)
eval.pts <- expand.grid(eval.x, eval.y)
## create list of matrices for different partial derivatives
## place vector of density estimate values `fhat' onto grid 'fhat.grid'
for (ell in 1:nderiv)
for (k in 1:length(eval.z))
{
fhat <- w[i]*dmvnorm.deriv(cbind(eval.pts, eval.z[k]), x[i,], H, deriv.order=r)
for (j in 1:length(eval.y))
fhat.grid[[ell]][eval.x.ind,eval.y.ind[j], eval.z.ind[k]] <-
fhat.grid[[ell]][eval.x.ind, eval.y.ind[j], eval.z.ind[k]] +
fhat[((j-1) * eval.x.len + 1):(j * eval.x.len), ell]
}
if (verbose) setTxtProgressBar(pb, i/n)
}
if (verbose) close(pb)
for (k in 1:nderiv) fhat.grid[[k]] <- fhat.grid[[k]]/n
gridx1 <- list(gridx[[1]], gridx[[2]], gridx[[3]])
ind.mat <- dmvnorm.deriv(x=rep(0,d), mu=rep(0,d), Sigma=H, deriv.order=r, only.index=TRUE)
if (!deriv.vec)
{
fhat.grid.vech <- list()
deriv.ind <- unique(ind.mat)
for (i in 1:nrow(deriv.ind))
{
which.deriv <- which.mat(deriv.ind[i,], ind.mat)[1]
fhat.grid.vech[[i]] <- fhat.grid[[which.deriv]]
}
ind.mat <- deriv.ind
fhat.grid <- fhat.grid.vech
}
fhatr <- list(x=x, eval.points=gridx1, estimate=fhat.grid, H=H, gridtype=gridx$gridtype, gridded=TRUE, binned=FALSE, names=NULL, w=w, deriv.order=deriv.order, deriv.ind=ind.mat)
}
return(fhatr)
}
#############################################################################
## Multivariate kernel density estimate using normal kernels,
## evaluated at each sample point
#############################################################################
kdde.points.1d <- function(x, h, eval.points, w, deriv.order=0)
{
r <- deriv.order
n <- length(x)
fhat <- dnorm.deriv.mixt(x=eval.points, mus=x, sigmas=rep(h,n), props=w/n, deriv.order=r)
return(list(x=x, eval.points=eval.points, estimate=fhat, h=h, H=h^2, gridded=FALSE, binned=FALSE, names=NULL, w=w, deriv.order=r, deriv.ind=r))
}
kdde.points <- function(x, H, eval.points, w, deriv.order=0, deriv.vec=TRUE)
{
n <- nrow(x)
Hs <- replicate(n, H, simplify=FALSE)
Hs <- do.call(rbind, Hs)
r <- deriv.order
fhat <- dmvnorm.deriv.mixt(x=eval.points, mus=x, Sigmas=Hs, props=w/n, deriv.order=r, deriv.vec=deriv.vec, add.index=TRUE)
return(list(x=x, eval.points=eval.points, estimate=fhat$deriv, H=H, gridded=FALSE, binned=FALSE, names=NULL, w=w, deriv.order=r, deriv.ind=fhat$deriv.ind))
}
#############################################################################
## S3 methods for KDDE objects
#############################################################################
## plot method
plot.kdde <- function(x, ...)
{
fhat <- x
if (is.null(fhat$deriv.order))
{
class(fhat) <- "kde"
plot(fhat, ...)
}
else
{
if (is.vector(fhat$x))
{
plotkdde.1d(fhat, ...)
invisible()
}
else
{
d <- ncol(fhat$x)
if (d==2)
{
opr <- options()$preferRaster; if (!is.null(opr)) if (!opr) options("preferRaster"=TRUE)
plotret <- plotkdde.2d(fhat, ...)
if (!is.null(opr)) options("preferRaster"=opr)
invisible(plotret)
}
else if (d==3)
{
plotkdde.3d(fhat, ...)
invisible()
}
else
stop ("Plot function only available for 1, 2 or 3-d data")
}
}
}
plotkdde.1d <- function(fhat, xlab, ylab="Density derivative function", cont=50, abs.cont, ...)
{
if (missing(xlab)) xlab <- fhat$names[1]
if (missing(abs.cont))
{
abs.cont <- as.matrix(contourLevels(fhat, approx=TRUE, cont=cont), ncol = length(cont))
abs.cont <- c(abs.cont[1, ], rev(abs.cont[2, ]))
}
class(fhat) <- "kde"
plot(fhat, xlab=xlab, ylab=ylab, abs.cont=abs.cont, ...)
}
plotkdde.2d <- function(fhat, which.deriv.ind=1, cont=c(25,50,75), abs.cont, display="slice", xlab, ylab, zlab="Density derivative function", col, col.fun, alpha=1, kdde.flag=TRUE, thin=3, transf=1, neg.grad=FALSE, ...)
{
disp <- match.arg(display, c("slice", "persp", "image", "filled.contour", "filled.contour2", "quiver"))
if (disp=="filled.contour2") disp <- "filled.contour"
if (missing(col.fun)) col.fun <- function(n) { hcl.colors(n, palette="Blue-Red", alpha=alpha) }
if (missing(xlab)) xlab <- fhat$names[1]
if (missing(ylab)) ylab <- fhat$names[2]
if (disp=="slice" | disp=="filled.contour")
{
if (missing(abs.cont))
{
abs.cont <- as.matrix(contourLevels(fhat, approx=TRUE, cont=cont, which.deriv.ind=which.deriv.ind), ncol=length(cont))
abs.cont <- c(abs.cont[1,], rev(abs.cont[2,]))
}
if (missing(col))
{
if (disp=="slice") col <- col.fun(n=length(abs.cont))
else if (disp=="filled.contour")
{
col <- col.fun(n=length(abs.cont)+1)
col[median(1:length(col))] <- "transparent"
}
}
}
if (disp=="quiver")
{
if (fhat$deriv.order==1)
plotquiver(fhat=fhat, thin=thin, transf=transf, neg.grad=neg.grad, col=col, xlab=xlab, ylab=ylab, alpha=alpha, ...)
else warning("Quiver plot requires gradient estimate.")
}
else
{
fhat.temp <- fhat
fhat.temp$deriv.ind <- fhat.temp$deriv.ind[which.deriv.ind,]
fhat.temp$estimate <- fhat.temp$estimate[[which.deriv.ind]]
fhat <- fhat.temp
class(fhat) <- "kde"
if (disp=="persp") plot(fhat, display=display, abs.cont=abs.cont, xlab=xlab, ylab=ylab, zlab=zlab, col.fun=col.fun, kdde.flag=kdde.flag, col=col, thin=thin, alpha=alpha, ...)
else plot(fhat, display=display, abs.cont=abs.cont, xlab=xlab, ylab=ylab, zlab=zlab, col.fun=col.fun, kdde.flag=kdde.flag, col=col, alpha=alpha, ...)
}
}
plotkdde.3d <- function(fhat, which.deriv.ind=1, display="plot3D", cont=c(25,50,75), abs.cont, colors, col, col.fun, xlab, ylab, zlab, ...)
{
if (missing(col.fun)) col.fun <- function(n) { hcl.colors(n, palette="Blue-Red") }
if (missing(xlab)) xlab <- fhat$names[1]
if (missing(ylab)) ylab <- fhat$names[2]
if (missing(zlab)) ylab <- fhat$names[3]
if (missing(abs.cont))
{
abs.cont <- as.matrix(contourLevels(fhat, approx=TRUE, cont=cont, which.deriv.ind=which.deriv.ind), ncol=length(cont))
abs.cont <- c(abs.cont[1,], rev(abs.cont[2,]))
}
fhat.temp <- fhat
fhat.temp$deriv.ind <- fhat.temp$deriv.ind[which.deriv.ind,]
fhat.temp$estimate <- fhat.temp$estimate[[which.deriv.ind]]
fhat <- fhat.temp
class(fhat) <- "kde"
if (missing(col))
{
col <- col.fun(n=length(abs.cont)+1)
nc <- length(col)
col <- col[-median(1:nc)]
}
colors <- col
plot(fhat, display=display, abs.cont=abs.cont, colors=col, col=col, xlab=xlab, ylab=ylab, zlab=zlab, ...)
}
######################################################################
## Quiver plot
######################################################################
plotquiver <- function(fhat, thin=5, transf=1, neg.grad=FALSE, xlab, ylab, col, add=FALSE, scale, length=0.1, alpha=1, ...)
{
if (!requireNamespace("pracma", quietly=TRUE)) stop("Install the pracma package as it is required.", call.=FALSE)
if (missing(col)) col <- 1
col <- transparency.col(col, alpha=alpha)
ev <- fhat$eval.points
est <- fhat$estimate
if (transf!=0)
{
est[[1]] <- sign(est[[1]])*abs(est[[1]])^(transf)
est[[2]] <- sign(est[[2]])*abs(est[[2]])^(transf)
}
thin1.ind <- seq(1, length(ev[[1]]), by=thin)
thin2.ind <- seq(1, length(ev[[2]]), by=thin)
evx <- ev[[1]][thin1.ind]
evy <- ev[[2]][thin2.ind]
fx <- est[[1]][thin1.ind, thin2.ind]
fy <- est[[2]][thin1.ind, thin2.ind]
if (neg.grad) { fx <- -fx; fy <- -fy }
if (missing(xlab)) xlab <- fhat$names[1]
if (missing(ylab)) ylab <- fhat$names[2]
if (!add) plot(fhat, abs.cont=c(0,0), col="transparent", xlab=xlab, ylab=ylab)
grid.xy <- pracma::meshgrid(evy, evx)
x0 <- grid.xy$Y
y0 <- grid.xy$X
## default scale factor - arrows should exceed a bin
arrow.len <- sqrt(fx^2 + fy^2)
bin.diag <- min(diff(evx), diff(evy))
if (missing(scale)) scale <- bin.diag/max(arrow.len)
## remove `zero-length' arrows i.e. length < 1e-3 inches
## adapted from https://stackoverflow.com/questions/52689959/how-small-is-a-zero-length-arrow/52690054
units <- par(c('usr', 'pin'))
x2in <- with(units, pin[1L]/diff(usr[1:2]))
y2in <- with(units, pin[2L]/diff(usr[3:4]))
arrow2in <- sqrt((x2in*scale*fx)^2 + (y2in*scale*fy)^2)
arrowint <- sort(c(min(arrow2in)-0.1*abs(max(arrow2in)), seq(1e-3, length, length=5), max(arrow2in)+0.1*abs(max(arrow2in))))
nsmall.f <- cut(arrow2in, arrowint, labels=FALSE)
nsf.lab <- unique(nsmall.f[nsmall.f>1])
for (i in nsf.lab)
{
nsf <- which(nsmall.f==i)
if (length(nsf)>0) pracma::quiver(x=x0[nsf], y=y0[nsf], u=fx[nsf], v=fy[nsf], col=col, scale=scale, length=arrowint[i], ...)
}
}
## contourLevels method
contourLevels.kdde <- function(x, prob, cont, nlevels=5, approx=TRUE, which.deriv.ind=1, ...)
{
fhat <- x
if (is.vector(fhat$x))
{
d <- 1; n <- length(fhat$x)
if (!is.null(fhat$deriv.order))
{
fhat.temp <- fhat
fhat.temp$deriv.ind <-fhat.temp$deriv.ind[which.deriv.ind]
fhat <- fhat.temp
}
}
else
{
d <- ncol(fhat$x); n <-nrow(fhat$x)
if (!is.matrix(fhat$x)) fhat$x <- as.matrix(fhat$x)
if (!is.null(fhat$deriv.order))
{
fhat.temp <- fhat
fhat.temp$estimate <- fhat.temp$estimate[[which.deriv.ind]]
fhat.temp$deriv.ind <-fhat.temp$deriv.ind[which.deriv.ind,]
fhat <- fhat.temp
}
}
if (is.null(x$w)) w <- rep(1, n)
else w <- x$w
if (is.null(fhat$gridded))
{
if (d==1) fhat$gridded <- fhat$binned
else fhat$gridded <- is.list(fhat$eval.points)
}
if (missing(prob) & missing(cont))
hts <- pretty(fhat$estimate, n=nlevels)
else
{
if (approx & fhat$gridded)
dobs <- predict.kde(fhat, x=fhat$x)
else
dobs <- kdde(x=fhat$x, H=fhat$H, eval.points=fhat$x, w=w, deriv.order=fhat$deriv.order)$estimate[,which.deriv.ind]
if (is.null(fhat$deriv.order))
{
if (!missing(prob) & missing(cont)) hts <- quantile(dobs[dobs>=0], prob=prob)
if (missing(prob) & !missing(cont)) hts <- quantile(dobs[dobs>=0], prob=(100-cont)/100)
}
else
{
if (!missing(prob) & missing(cont)) hts <- rbind(-quantile(abs(dobs[dobs<0]), prob=prob), quantile(dobs[dobs>=0], prob=prob))
if (missing(prob) & !missing(cont)) hts <- rbind(-quantile(abs(dobs[dobs<0]), prob=(100-cont)/100), quantile(dobs[dobs>=0], prob=(100-cont)/100))
}
}
return(hts)
}
## predict method for KDDE
predict.kdde <- function(object, ..., x)
{
fhat <- object
if (is.vector(fhat$H)) d <- 1 else d <- ncol(fhat$H)
if (d==1) n <- length(x)
else
{
if (is.vector(x)) x <- matrix(x, nrow=1)
else x <- as.matrix(x)
n <- nrow(x)
}
if (!is.null(fhat$deriv.ind))
{
if (is.vector(fhat$deriv.ind)) pk.mat <- predict.kde(fhat, x=x, ...)
else
{
nd <- nrow(fhat$deriv.ind)
pk.mat <- matrix(0, ncol=nd, nrow=n)
for (i in 1:nd)
{
fhat.temp <- fhat
fhat.temp$estimate <- fhat$estimate[[i]]
pk.mat[,i] <- predict.kde(fhat.temp, x=x, ...)
}
}
}
else
pk.mat <- predict.kde(fhat, x=x, ...)
return(drop(pk.mat))
}
######################################################################
## Summary kernel curvature
######################################################################
kcurv <- function(fhat, compute.cont=TRUE)
{
fhat.curv <- fhat
if (is.vector(fhat$H)) d <- 1 else d <- ncol(fhat$H)
if (fhat$deriv.order!=2) stop("Requires output from kdde(, deriv.order=2).")
if (d==1)
{
Hessian.det <- fhat$estimate
local.mode <- fhat$estimate <0
fhat.curv$estimate <- local.mode*abs(Hessian.det)
}
else if (d>1)
{
fhat.est <- sapply(fhat$estimate, as.vector)
Hessian.det <- sapply(seq(1,nrow(fhat.est)), function(i) { det(invvec(fhat.est[i,])) })
Hessian.eigen <- lapply(lapply(seq(1,nrow(fhat.est)), function(i) { invvec(fhat.est[i,]) }), eigen, only.values=TRUE)
Hessian.eigen <- t(sapply(Hessian.eigen, getElement, "values"))
local.mode <- apply(Hessian.eigen <= 0, 1, all)
fhat.curv$estimate <- local.mode*array(abs(Hessian.det), dim=dim(fhat$estimate[[1]]))
}
fhat.curv$deriv.order <- NULL
fhat.curv$deriv.ind <- NULL
if (compute.cont)
{
fhat.temp <- fhat.curv
fhat.temp$x <- fhat.curv$x[predict(fhat.curv, x=fhat.curv$x)>0,]
fhat.temp$estimate <- fhat.temp$estimate
fhat.curv$cont <- contourLevels(fhat.temp, cont=1:99)
}
class(fhat.curv) <- "kde"
fhat.curv$type <- "kcurv"
return(fhat.curv)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.