Nothing
########################################################################
## Default grid sizes
########################################################################
default.gridsize <- function(d)
{
if (d==1) gridsize <- 401
else if (d==2) gridsize <- rep(151,d)
else if (d==3) gridsize <- rep(51, d)
else if (d>=4) gridsize <- rep(21, d)
return(gridsize)
}
default.bgridsize <- function(d)
{
if (d==1) gridsize <- 401
else if (d==2) gridsize <- rep(151,d)
else if (d==3) gridsize <- rep(31, d)
else if (d==4) gridsize <- rep(15, d)
else gridsize <- NA
return(gridsize)
}
default.bflag <- function(d, n)
{
if (d==1) thr <- 1
else if (d==2) thr <- 500
else if (d>2) thr <- 1000
bf <- n>thr
return(bf)
}
## truncate x to xmin, xmax grid
truncate.grid <- function(x, y, xmin, xmax)
{
if (is.vector(x)) { d <- 1; n <- length(x); xvec <- TRUE; x <- matrix(x, ncol=1) }
else { d <- ncol(x); n <- nrow(x); xvec <- FALSE }
xind <- rep(TRUE, n)
if (!(missing(xmax)))
{
if (d==1) msgmax <- xmax else msgmax <- paste0("c(", paste0(xmax, collapse=","), ")")
if (any(sweep(x, 2, FUN=">", xmax))) warning(paste0("Points in x greater than xmax=", msgmax, " have been excluded."))
for (i in 1:d) xind <- xind & x[,i]<=xmax[i]
}
if (!(missing(xmin)))
{
if (d==1) msgmin <- xmin else msgmin <- paste0("c(", paste0(xmin, collapse=","), ")")
if (any(sweep(x, 2, FUN="<", xmin))) warning(paste0("Points in x less than xmin=", msgmin, " have been excluded."))
for (i in 1:d) xind <- xind & x[,i]>=xmin[i]
}
if (d==1)
{
x <- x[xind,, drop=xvec]
## special case for x is 1-col matrix to force x to be
## 1-col matrix, as required for eks <= 1.0.1
}
else
{
x <- x[xind,, drop=FALSE]
}
if (!missing(y)) { y <- y[xind]; x <- list(x=x, y=y) }
return(x)
}
########################################################################
## Linear binning
## Courtesy of M Wand 2005
## Extended by T Duong to 3- and 4-dim 2006
## Extended by Gramack & Gramacki to include unconstrained b/w 2015
########################################################################
binning <- function(x, H, h, bgridsize, xmin, xmax, supp=3.7, w, gridtype="linear")
{
## default values
x <- as.matrix(x)
d <- ncol(x)
n <- nrow(x)
if (missing(bgridsize)) bgridsize <- default.gridsize(d)
if (missing(w)) w <- rep(1,n)
if (missing(h)) h <- rep(0,d)
if (!missing(H)) h <- sqrt(diag(H))
range.x <- list()
if (!missing(xmin) & !missing(xmax))
for (i in 1:d) range.x[[i]] <- c(xmin[i], xmax[i])
else if (!missing(xmin) & missing(xmax))
for (i in 1:d) range.x[[i]] <- c(xmin[i], max(x[,i]) + supp*h[i])
else if (missing(xmin) & !missing(xmax))
for (i in 1:d) range.x[[i]] <- c(min(x[,i]) - supp*h[i], xmax[i])
else
for (i in 1:d) range.x[[i]] <- c(min(x[,i]) - supp*h[i], max(x[,i]) + supp*h[i])
a <- sapply(range.x,min)
b <- sapply(range.x,max)
if (missing(gridtype)) gridtype <- rep("linear", d)
gridtype.vec <- rep("", d)
gpoints <- list()
for (id in 1:d)
{
gridtype1 <- match.arg(gridtype[i], c("linear", "sqrt", "quantile", "log"))
if (gridtype1=="linear")
gpoints[[id]] <- seq(a[id],b[id],length=bgridsize[id])
else if (gridtype1=="log")
gpoints[[id]] <- seq(exp(a[id]),exp(b[id]),length=bgridsize[id])
}
if (d==1) counts <- linbin.ks(x,gpoints[[1]], w=w)
if (d==2) counts <- linbin2D.ks(x,gpoints[[1]],gpoints[[2]], w=w)
if (d==3) counts <- linbin3D.ks(x,gpoints[[1]],gpoints[[2]],gpoints[[3]], w=w)
if (d==4) counts <- linbin4D.ks(x,gpoints[[1]],gpoints[[2]],gpoints[[3]],gpoints[[4]], w=w)
bin.counts <- list(counts=counts, eval.points=gpoints, w=w)
if (d==1) bin.counts <- lapply(bin.counts, unlist)
return(bin.counts)
}
########################################################################
## Linear binning
########################################################################
linbin.ks <- function(x, gpoints, w)
{
n <- length(x)
M <- length(gpoints)
if (missing(w)) w <- rep(1, n)
a <- gpoints[1]
b <- gpoints[M]
xi <- .C(C_massdist1d, x1=as.double(x[,1]), n=as.integer(n), a1=as.double(a), b1=as.double(b), M1=as.integer(M), weight=as.double(w), est=double(M))$est
return(xi)
}
linbin2D.ks <- function(x, gpoints1, gpoints2, w)
{
n <- nrow(x)
M1 <- length(gpoints1)
M2 <- length(gpoints2)
a1 <- gpoints1[1]
a2 <- gpoints2[1]
b1 <- gpoints1[M1]
b2 <- gpoints2[M2]
if (missing(w)) w <- rep(1, n)
## binning for interior points
out <- .C(C_massdist2d, x1=as.double(x[,1]), x2=as.double(x[,2]), n=as.integer(n), a1=as.double(a1), a2=as.double(a2), b1=as.double(b1), b2=as.double(b2), M1=as.integer(M1), M2=as.integer(M2), weight=as.double(w), est=double(M1*M2))
xi <- matrix(out$est, nrow=M1, ncol=M2)
return(xi)
}
linbin3D.ks <- function(x, gpoints1, gpoints2, gpoints3, w)
{
n <- nrow(x)
M1 <- length(gpoints1)
M2 <- length(gpoints2)
M3 <- length(gpoints3)
a1 <- gpoints1[1]
a2 <- gpoints2[1]
a3 <- gpoints3[1]
b1 <- gpoints1[M1]
b2 <- gpoints2[M2]
b3 <- gpoints3[M3]
if (missing(w)) w <- rep(1, n)
## binning for interior points
out <- .C(C_massdist3d, x1=as.double(x[,1]), x2=as.double(x[,2]), x3=as.double(x[,3]), n=as.integer(n), a1=as.double(a1), a2=as.double(a2), a3=as.double(a3), b1=as.double(b1), b2=as.double(b2), b3=as.double(b3), M1=as.integer(M1), M2=as.integer(M2), M3=as.integer(M3), weight=as.double(w), est=double(M1*M2*M3))
xi <- array(out$est, dim=c(M1,M2,M3))
return(xi)
}
linbin4D.ks <- function(x, gpoints1, gpoints2, gpoints3, gpoints4, w)
{
n <- nrow(x)
M1 <- length(gpoints1)
M2 <- length(gpoints2)
M3 <- length(gpoints3)
M4 <- length(gpoints4)
a1 <- gpoints1[1]
a2 <- gpoints2[1]
a3 <- gpoints3[1]
a4 <- gpoints4[1]
b1 <- gpoints1[M1]
b2 <- gpoints2[M2]
b3 <- gpoints3[M3]
b4 <- gpoints4[M4]
if (missing(w)) w <- rep(1, n)
## binning for interior points
out <- .C(C_massdist4d, x1=as.double(x[,1]), x2=as.double(x[,2]), x3=as.double(x[,3]), x4=as.double(x[,4]), n=as.integer(n), a1=as.double(a1), a2=as.double(a2), a3=as.double(a3), a4=as.double(a4), b1=as.double(b1), b2=as.double(b2), b3=as.double(b3), b4=as.double(b4), M1=as.integer(M1), M2=as.integer(M2), M3=as.integer(M3), M4=as.integer(M4), weight=as.double(w), est=double(M1*M2*M3*M4))
xi <- array(out$est, dim=c(M1,M2,M3,M4))
return(xi)
}
########################################################################
## Discrete convolution
########################################################################
symconv.1d <- function(keval, gcounts)
{
M <- length(gcounts)
N <- length(keval)
L <- (length(keval)+1)/2
## Smallest powers of 2 >= M + N
P <- 2^(ceiling(log2(M + N)))
## Zero-padded version of keval an gcounts
keval.zeropad <- rep(0, P)
gcounts.zeropad <- rep(0, P)
keval.zeropad[1:N] <- keval
gcounts.zeropad[L:(L+M-1)] <- gcounts
## FFTs
K <- fft(keval.zeropad)
C <- fft(gcounts.zeropad)
## Invert element-wise product of FFTs and truncate and normalise it
symconv.val <- Re(fft(K*C, inverse=TRUE)/P)[N:(N+M-1)]
return(symconv.val)
}
symconv.nd <- function(keval, gcounts, d)
{
M <- dim(gcounts)
N <- dim(keval)
L <- (dim(keval)+1)/2
## Smallest powers of 2 >= M + N
P <- 2^(ceiling(log2(M + N)))
## Zero-padded version of keval and gcounts
keval.zeropad <- array(0, dim=P)
gcounts.zeropad <- array(0, dim=P)
if (d==2)
{
keval.zeropad[1:N[1], 1:N[2]] <- keval
gcounts.zeropad[L[1]:(L[1]+M[1]-1), L[2]:(L[2]+M[2]-1)] <- gcounts
}
else if (d==3)
{
keval.zeropad[1:N[1], 1:N[2], 1:N[3]] <- keval
gcounts.zeropad[L[1]:(L[1]+M[1]-1), L[2]:(L[2]+M[2]-1), L[3]:(L[3]+M[3]-1)] <- gcounts
}
else if (d==4)
{
keval.zeropad[1:N[1], 1:N[2], 1:N[3], 1:N[4]] <- keval
gcounts.zeropad[L[1]:(L[1]+M[1]-1), L[2]:(L[2]+M[2]-1), L[3]:(L[3]+M[3]-1), L[4]:(L[4]+M[4]-1) ] <- gcounts
}
## FFTs
K <- fft(keval.zeropad)
C <- fft(gcounts.zeropad)
## Invert element-wise product of FFTs and truncate and normalise it
symconv.val <- Re(fft(K*C, inverse=TRUE)/prod(P))
if (d==2) symconv.val <- symconv.val[N[1]:(N[1]+M[1]-1), N[2]:(N[2]+M[2]-1)]
else if (d==3) symconv.val <- symconv.val[N[1]:(N[1]+M[1]-1), N[2]:(N[2]+M[2]-1), N[3]:(N[3]+M[3]-1)]
else if (d==4) symconv.val <- symconv.val[N[1]:(N[1]+M[1]-1), N[2]:(N[2]+M[2]-1), N[3]:(N[3]+M[3]-1), N[4]:(N[4]+M[4]-1)]
return(symconv.val)
}
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.