R/prelim.R

In ks: Kernel Smoothing

##########################################################################
## Default values for parameter choices for ks objects
##########################################################################

ks.defaults <- function(x, w, binned, bgridsize, gridsize)
{
    ## dimensions of x
    if (is.vector(x)) { d <- 1; n <- length(x) }
    else { d <- ncol(x); n <- nrow(x) }

    ## default uniform weights
    if (missing(w)) w <- rep(1,n)
    else if (!missing(w))
    {
        if (!(identical(all.equal(sum(w), n), TRUE)))
        {
            warning("Weights don't sum to sample size - they have been scaled accordingly\n")
            w <- w*n/sum(w)
        }
    }

    ## default binning flag
    if (missing(binned)) binned <- default.bflag(d=d, n=n)

    ## default grid sizes
    if (missing(bgridsize))
    {
        if (missing(gridsize)) bgridsize <- default.bgridsize(d)
        else bgridsize <- gridsize
    }
    if (missing(gridsize)) gridsize <- default.gridsize(d)
    if (length(gridsize)==1) gridsize <- rep(gridsize, d)
    if (length(bgridsize)==1) bgridsize <- rep(bgridsize, d)
    ksd <- list(d=d, n=n, w=w, binned=binned, bgridsize=bgridsize, gridsize=gridsize)

    return(ksd)
}

##########################################################################
## Parse variable name
##########################################################################

parse.name <- function(x)
{
    if (is.vector(x))
    {
        d <- 1
        x.names <- deparse(substitute(x))
    }
    else
    {  
        d <- ncol(x)
        x.names <- colnames(x)
        if (is.null(x.names))
        {
            x.names <- strsplit(deparse(substitute(x)), "\\[")[[1]][1]
            x.names <- paste(x.names, "[, ", 1:d,"]",sep="") 
        }
    }

    return(x.names)
}

##########################################################################
## Basic vectors and matrices and their operations
##########################################################################
  
## vec operator 

vec <- function(x, byrow=FALSE)
{
    if (is.vector(x)) return (x)
    if (byrow) x <- t(x)
    d <- ncol(x)
    vecx <- vector()
    for (j in 1:d) vecx <- c(vecx, x[,j])

    return(vecx)           
}

## vech operator

vech <- function(x)
{
    if (is.vector(x))
    {
    if (length(x)==1)
        return (x)
    else
        stop("vech undefined for vectors")
    }
    else if (is.matrix(x))
    {  
        d <- ncol(x)
        if (d!=nrow(x)) stop("vech only defined for square matrices")

        vechx <- vector()
        for (j in 1:d) vechx <- c(vechx, x[j:d,j])
        return(vechx)
    }
}

## inverse vec operator 
    
invvec <- function(x, ncol, nrow, byrow=FALSE)
{
    if (length(x)==1) return(x)

    d <- sqrt(length(x))
    if (missing(ncol) | missing(nrow))
    {
        ncol <- d; nrow <- d
        if (round(d) != d)
            stop("Need to specify nrow and ncol for non-square matrices")
    }

    invvecx <- matrix(0, nrow = nrow, ncol = ncol)
    if (byrow)
    for (j in 1:nrow)
        invvecx[j,] <- x[c(1:ncol) + (j-1)*ncol]
    else
        for (j in 1:ncol)
            invvecx[,j] <- x[c(1:nrow) + (j-1)*nrow]

    return(invvecx)
}

## inverse vech operator 

invvech <- function(x)
{
  if (length(x)==1) return(x)
  
  d <- (-1 + sqrt(8*length(x) + 1))/2
  if (round(d) != d)
    stop("Number of elements in x will not form a square matrix")
  invvechx <- matrix(0, nrow=d, ncol=d)

  for (j in 1:d)
        invvechx[j:d,j] <- x[1:(d-j+1)+ (j-1)*(d - 1/2*(j-2))]
  
  invvechx <- invvechx + t(invvechx) - diag(diag(invvechx))
  
  return(invvechx)
}

## trace of matrix

tr <- function(A)
{   
    count <- 0
    if (is.vector(A)) return (A[1])
    if (nrow(A)!=ncol(A)) stop("Not square matrix")
    else 
        for (i in 1:nrow(A)) count <- count + A[i,i]
  
    return(count)
}

## elementary vector 
    
elem <- function(i, d)
{
    elem.vec <- rep(0, d)
    elem.vec[i] <- 1

    return(elem.vec)
}      

## commutation matrix (taken from MCMCglmmm library)

comm <- function(m,n)
{
    K <- matrix(0,m*n, m*n)
    H <- matrix(0,m,n)
    for (i in 1:m)
    {
        for(j in 1:n)
        { 
            H[i,j] <- 1
            K <- K+kronecker(H,t(H))
            H[i,j] <- 0
        }
    }

    return(K)
}

##########################################################################
## Duplication matrix
## Taken from Felipe Osorio http://www.ime.usp.br/~osorio/files/dupl.q
##########################################################################

dupl <- function(order, ret.q = FALSE)
{
    ## call
    cl <- match.call()
    time1 <- proc.time()
    if (!is.integer(order))
        order <- as.integer(order)
    n <- order - 1
    
    ## initial duplication matrix
    d1 <- matrix(0, nrow = 1, ncol = 1)
    d1[1,1] <- 1
    if (!is.integer(d1))
        storage.mode(d1) <- "integer"
    
    ## recursive formula
    if (n > 0)
    {
        for (k in 1:n)
        {
            drow <- 2*k + 1 + nrow(d1)
            dcol <- k + 1 + ncol(d1)
            d2 <- matrix(0, nrow = drow, ncol=dcol)
            storage.mode(d2) <- "integer"
            d2[1,1] <- 1
            d2[2:(k+1),2:(k+1)] <- diag(k)
            d2[(k+2):(2*k+1),2:(k+1)] <- diag(k)
            d2[(2*k+2):drow,(k+2):dcol] <- d1
            ## permutation matrix
            q <- permute.mat(k)
            ## new duplication matrix
            d2 <- q %*% d2
            storage.mode(d2) <- "integer"
            d1 <- d2
        }
    }
    else {
        d2 <- q <- d1
    }
    
    ## results
    obj <- list(call=cl, order=order, d=d2)
    if (ret.q)
        obj$q <- q
    obj$time <- proc.time() - time1
    obj
}

##########################################################################
## Pre-scaling
## Parameters
## x - data points
##
## Returns
## Pre-scaled x values
##########################################################################

pre.scale <- function(x, mean.centred=FALSE)
{
    S <- diag(diag(var(x)))
    Sinv12 <- matrix.sqrt(chol2inv(chol(S)))

    if (mean.centred) x.scaled <- sweep(x, 2, apply(x, 2, mean))
    else x.scaled <- x  
    x.scaled <- x.scaled %*% Sinv12

    return (x.scaled)
}

##########################################################################
## Pre-sphering
## Parameters
## x - data points
##
## Returns
## Pre-sphered x values
##########################################################################

pre.sphere <- function(x, mean.centred=FALSE)
{
    S <- var(x)
    Sinv12 <- matrix.sqrt(chol2inv(chol(S)))

    if (mean.centred) x.sphered <- sweep(x, 2, apply(x, 2, mean))
    else x.sphered <- x  
    x.sphered <- x.sphered %*% Sinv12

  return (x.sphered)
}

##########################################################################
## Boolean functions
##########################################################################

is.even <- function(x)
{
    y <- x[x>0] %%2
    return(identical(y, rep(0, length(y))))
}

is.diagonal <- function(x)
{
    return(identical(diag(diag(x)),x))
}

##########################################################################
## Finds row index matrix
## Parameters
## x - data points
##
## Returns
## i  - if r==mat[i,]
## NA - otherwise
##########################################################################

which.mat <- function(r, mat)
{
    ind <- numeric()

    for (i in 1:nrow(mat))
    if (identical(r, mat[i,])) ind <- c(ind,i)

    return(ind)  
}

###################################################################
## Permutation functions
###################################################################

####################################################################
## Exactly the same function as combinat:::permn
####################################################################

permn.ks <- function (x, fun = NULL, ...) 
{
    if (is.numeric(x) && length(x) == 1 && x > 0 && trunc(x) == x) x <- seq(x)
    n <- length(x)
    nofun <- is.null(fun)
    out <- vector("list", gamma(n + 1))
    p <- ip <- seqn <- 1:n
    d <- rep(-1, n)
    d[1] <- 0
    m <- n + 1
    p <- c(m, p, m)
    i <- 1
    use <- -c(1, n + 2)
    while (m != 1) {
        out[[i]] <- if (nofun) 
            x[p[use]]
        else fun(x[p[use]], ...)
        i <- i + 1
        m <- n
        chk <- (p[ip + d + 1] > seqn)
        m <- max(seqn[!chk])
        if (m < n) 
            d[(m + 1):n] <- -d[(m + 1):n]
        index1 <- ip[m] + 1
        index2 <- p[index1] <- p[index1 + d[m]]
        p[index1 + d[m]] <- m
        tmp <- ip[index2]
        ip[index2] <- ip[m]
        ip[m] <- tmp
    }
    out
}

##########################################################################
## Permutations with repetitions of the first d naturals (1:d) taking 
## k elements at a time. There are d^k of them, each having length k 
## => We arrange them into a matrix of order d^k times k
## Each row represents one permutation
## Second version: filling in the matrix comlumn-wise (slightly faster)
##########################################################################

perm.rep <- function(d,r)
{
    if (r==0) { PM <- 1 }
    if (r>0)
    {
        PM <- matrix(nrow=d^r,ncol=r)
        for (pow in 0:(r-1))
        {
            t2 <- d^pow
            p1 <- 1
            while (p1<=d^r)
            {
                for (al in 1:d)
                {
                    for (p2 in 1:t2)
                    {
                        PM[p1,r-pow] <- al
                        p1 <- p1+1
                    }
                }
            }
        }
    }

    return(PM)
}

##########################################################################
## Permute a list of values
##
## Same function as EXPAND.GRID (base package), modified to take 
## list as an argument and returns a matrix 
##########################################################################

permute <- function (args) 
{
    nargs <- length(args)
    if (!nargs) 
        return(as.data.frame(list()))
    if (nargs == 1 && is.list(a1 <- args[[1]])) 
    nargs <- length(args <- a1)
    if (nargs <= 1) 
        return(as.data.frame(if (nargs == 0 || is.null(args[[1]])) list() else args, optional = TRUE))
    
    cargs <- args
    rep.fac <- 1
    orep <- prod(sapply(args, length))

    for (i in 1:nargs) 
    {
        x <- args[[i]]
        nx <- length(x)
        orep <- orep/nx
        cargs[[i]] <- rep(rep(x, rep(rep.fac, nx)), orep)
        rep.fac <- rep.fac * nx
    }
    do.call("cbind", cargs)
} 

permute.mat <- function(order)
{
    m <- as.integer(order)
    m <- m + 1
    eye <- diag(m)
    u1 <- eye[1:m,1]
    u2 <- eye[1:m,2:m]
    q1 <- kronecker(eye, u1)
    q2 <- kronecker(eye, u2)
    q <- matrix(c(q1, q2), nrow = nrow(q2), ncol = ncol(q1) + ncol(q2))
    if (!is.integer(q)) storage.mode(q) <- "integer"
    
    return(q)
}

##########################################################################
## pinv.all generates all the permutations PR_{d,r} as described in
## Appendix B of Chacon and Duong (2014)
##########################################################################
    
pinv.all <- function(d,r)
{
    i <- 1:d^r
    n <- i-1
    dpow <- d^(0:r)
    n.mat <- matrix(rep(n,r+1),byrow=FALSE,nrow=d^r,ncol=r+1)
    dpow.mat <- matrix(rep(dpow,d^r),byrow=TRUE,nrow=d^r,ncol=r+1)
    ndf.mat <- floor(n.mat/dpow.mat)
    ans <- ndf.mat[,r:1]-d*ndf.mat[,(r+1):2]    

    return(ans+1)
} 

##########################################################################
## Block indices for double sums
##########################################################################

block.indices <- function(nx, ny, d, r=0, diff=FALSE, block.limit=1e6, npergroup)
{
    if (missing(npergroup)) 
    { 
        if (diff) npergroup <- max(c(block.limit %/% (nx*d^r), 1))
        else npergroup <- max(c(block.limit %/% nx,1))
    }
    nseq <- seq(1, ny, by=npergroup)
    if (tail(nseq,n=1) <= ny) nseq <- c(nseq, ny+1)
    if (length(nseq)==1) nseq <- c(1, ny+1)
    
    return(nseq)
}

block.indices2 <- function(nx, ny, block.limit=1e6, npergroup)
{
    if (missing(npergroup)) npergroup <- max(c(block.limit %/% nx,1))
    nseq <- seq(1, ny, by=npergroup)
    if (tail(nseq,n=1) <= ny) nseq <- c(nseq, ny+1)
    if (length(nseq)==1) nseq <- c(1, ny+1)
    
    return(nseq)
}

####################################################################
## Differences for double sums calculations
####################################################################

differences <- function(x, y, upper=FALSE, ff=FALSE, Kpow=0)
{
    if (missing(y)) y <- x
    if (is.vector(x)) x <- t(as.matrix(x))
    if (is.vector(y)) y <- t(as.matrix(y))

    nx <- nrow(x)
    ny <- nrow(y)
    d <- ncol(x)

    if (ff) difs <- ff(init=0, dim=c(nx*ny,d))
    else difs <- matrix(ncol=d,nrow=nx*ny)

    for (j in 1:d)
    {
        difs[,j] <- rep(x[,j], times=ny) - rep(y[1:ny,j], each=nx)
        ## jth column of difs contains all the differences X_{ij}-Y_{kj}
    }

    if (upper)
    {
        ind.remove <- numeric()
        for (j in 1:(nx-1))
          ind.remove <- c(ind.remove, (j*nx+1):(j*nx+j))
          
        return(difs[-ind.remove,])
    }
    else return(difs)
}

##########################################################################
## Odd factorial
##########################################################################

OF <- function(m) { factorial(m)/(2^(m/2)*factorial(m/2)) } 

##########################################################################
## Matrix square root - taken from Stephen Lake 
## http://www5.biostat.wustl.edu/s-news/s-news-archive/200109/msg00067.html
##########################################################################

matrix.sqrt <- function(A)
{
    if (length(A)==1) return(sqrt(A))
    sva <- svd(A)
    if (min(sva$d)>=0) Asqrt <- sva$u %*% diag(sqrt(sva$d)) %*% t(sva$v)
    else stop("Matrix square root is not defined")

    return(Asqrt)
}

##########################################################################
## Matrix power
##########################################################################

matrix.pow <- function(A, n)
{
    if (nrow(A)!=ncol(A)) stop("A must be a square matrix")
    if (floor(n)!=n) stop("n must be an integer")
    if (n==0) return(diag(ncol(A)))
    if (n < 0) return(matrix.pow(A=chol2inv(chol(A)), n=-n))
        
    ## trap non-integer n and return an error
    if (n == 1) return(A)
    result <- diag(1, ncol(A))
    while (n > 0) 
    {
        if (n %% 2 != 0) 
        {
          result <- result %*% A
          n <- n - 1
        }
        A <- A %*% A
        n <- n / 2
    }
    
    return(result)
}

##########################################################################
## Kmat computes the commutation matrix of orders m,n
##########################################################################
    
Kmat <- function(m, n)
{ 
    K <- matrix(0,ncol=m*n,nrow=m*n)
    i <-1:m; j <- 1:n
    rows <- rowSums(expand.grid((i-1)*n,j))
    cols <- rowSums(expand.grid(i,(j-1)*m))    
    positions <- cbind(rows,cols)
    K[positions] <- 1

    return(K)
}

##########################################################################
## mat.Kprod computes row-wise Kronecker products of matrices
## Returns a matrix with rows U[i,]%x%V[i,]
##########################################################################

mat.Kprod <- function(U,V)
{ 
    n1 <- nrow(U)
    n2 <- nrow(V)
    if (n1!=n2) stop("U and V must have the same number of vectors")
    p <- ncol(U)
    q <- ncol(V)
    onep <- rep(1,p)
    oneq <- rep(1,q)
    P <- (U%x%t(oneq))*(t(onep)%x%V)

    return(P)
}

##########################################################################
## Kpow computes the Kronecker power of a matrix A
##########################################################################

Kpow <- function(A,pow)
{    
    if (floor(pow)!=pow) stop("pow must be an integer")
    Apow <- A
    if (pow==0) { Apow <- 1 }
    if (pow>1) { for(i in 2:pow) Apow <- Apow%x%A }

    return(Apow)
} 

##########################################################################
## Kronecker sum
##########################################################################

Ksum <- function(A,B)
{
    AB <- numeric()
    for (i in 1:nrow(A))
        for (j in 1:nrow(B))
            AB <- rbind(AB, A[i,] + B[j,])

    return(AB)
}

##########################################################################
## Row-wise Kronecker product
##########################################################################

rowKpow <- function(A,B,r=1,s=1)
{
    A <- as.matrix(A)
    if (missing(B))
    {
        res <- t(apply(t(A), 2, Kpow, r))
    }
    else
    {    
        B <- as.matrix(B)
        if (nrow(A)!=nrow(B)) stop("A and B must have same number of rows") 
        res <- numeric()
        
        Ar <- t(apply(t(A), 2, Kpow, r))
        Bs <- t(apply(t(B), 2, Kpow, s))
        for (i in 1:ncol(Ar)) res <- cbind(res, apply(Bs, 2, FUN="*", Ar[,i]))
    }
    
    return(drop(res))
}

getRow <- function(object, n) { return(object[n,]) }

## Returns a matrix with the pow-th Kronecker power of A[i,] in the i-th row

mat.Kpow <- function(A,pow)
{ 
    Apow <- A
    if (pow==0) { Apow <- matrix(1,nrow=nrow(A), ncol=1) }  
    if (pow>1) 
    {
        for(i in 2:pow) Apow <- mat.Kprod(Apow,A)
    }

    return(Apow)
}

## Vector of all r-th partial derivatives of the normal density at x=0, i.e., D^{\otimes r)\phi(0)

DrL0 <- function(d,r)
{
    v <- as.vector(Kpow(A=vec(diag(d)),pow=r/2))
    DL0 <- (-1)^(r/2)*(2*pi)^(-d/2)*OF(r)*matrix(Sdrv(d=d, r=r, v=v), ncol=1) 
    
    return(DL0)
}

##########################################################################
## Symmetriser matrix
##########################################################################

## Wrapper functions for Chacon & Duong (2014) 

Sdr <- function(d, r, type="recursive")
{
    type <- match.arg(type, c("recursive", "direct"))
    Sdr.mat <- do.call(paste("Sdr", type, sep="."), list(d=d, r=r))

    return(Sdr.mat)
}

Sdrv <- function(d, r, v, type="recursive")
{
  type <- match.arg(type, c("recursive", "direct"))
  v <- as.vector(v)
  Sdrvec <- do.call(paste("Sdrv", type, sep="."), list(d=d, r=r, v=v))
  
  return(Sdrvec)
}

##########################################################################
## Sdr.direct computes the symmetrizer matrix S_{d,r} based on Equation (4)
## as described in Section 3 of Chacon and Duong (2014)
##########################################################################

Sdr.direct <- function(d,r)
{ 
    S <- matrix(0,ncol=d^r,nrow=d^r)
    per <- permn.ks(r)
    per.rep <- pinv.all(d,r)
    nper <- factorial(r)
    nper.rep <- d^r
    per <- matrix(unlist(per), byrow=TRUE, ncol=r, nrow=nper)    
    pow <- 0:(r-1)
    dpow <- d^pow
    
    if (nper.rep<=nper)
    {
        dpow.mat <- matrix(rep(dpow,nper),byrow=TRUE,ncol=r,nrow=nper)
        for(i in 1:nper.rep)
        {      
            ## Loop over no. perms with reps (d^r)
            pinvi <- per.rep[i,]
            sigpinvi <- matrix(pinvi[per],byrow=FALSE,nrow=nrow(per),ncol=ncol(per))
            psigpinvi <- drop(1+rowSums((sigpinvi-1)*dpow.mat))
            S[i,] <- tabulate(psigpinvi,nbins=d^r)
        }
    }
    
    if (nper<nper.rep)
    {
        dpow.mat <- matrix(rep(dpow,nper.rep),byrow=TRUE,ncol=r,nrow=nper.rep)
        for (s in 1:nper)
        {          
            ## Loop over no. perms (r!)
            sig <- per[s,]
            sigpinv <- per.rep[,sig]
            psigpinv <- drop(1+rowSums((sigpinv-1)*dpow.mat))
            locations <- cbind(1:d^r,psigpinv)
            S[locations] <- S[locations]+1
        }
    }    
    
    return(S/nper)
}
     
##########################################################################
## Sdr.recursive computes the symmetrizer matrix S_{d,r} based on 
## the recursive approach detailed in Algorithm 1 in Section 3 of 
## Chacon and Duong (2014)
##########################################################################
     
Sdr.recursive <- function(d,r)
{
    S <- diag(d)
    if (r==0) S <- 1
    if (r>=2) 
    {
        Id <- diag(d)
        T <- Id 
        A <- Kmat(d,d)        
        for (j in 2:r)
        {
            T <- ((j-1)/j)*(A%*%(T%x%Id)%*%A)+A/j
            S <- (S%x%Id)%*%T
            if (j<r) A <- Id%x%A
        }
    }
    
    return(S)
}

##########################################################################
## Symmetrizer matrix applied to a vector
##########################################################################
    

##########################################################################
## Sdrv.direct computes the result of multiplying the symmetrizer matrix
## S_{d,r} by a vector v of length d^r, based on Equation (4) as described
## in Section 4 of Chacon and Duong (2014)
##########################################################################

Sdrv.direct <- function(d,r,v)
{
    Sv <- rep(0,d^r)
    per <- permn.ks(r)
    per.rep <- pinv.all(d,r)
    nper <- factorial(r)
    nper.rep <- d^r
    per <- matrix(unlist(per), byrow=TRUE, ncol=r, nrow=nper)    
    pow <- 0:(r-1)
    dpow <- d^pow
    
    if (nper.rep<=nper)
    {
        dpow.mat <- matrix(rep(dpow,nper),byrow=TRUE,ncol=r,nrow=nper)
        for (i in 1:nper.rep)
        {      
            ## Loop over no. perms with reps (d^r)
            pinvi <- per.rep[i,]
            sigpinvi <- matrix(pinvi[per],byrow=FALSE,nrow=nrow(per),ncol=ncol(per))
            psigpinvi <- drop(1+rowSums((sigpinvi-1)*dpow.mat))
            Sv[i] <- sum(tabulate(psigpinvi,nbins=d^r)*v)
        }
    }
    
    if (nper<nper.rep)
    {
        dpow.mat <- matrix(rep(dpow,nper.rep),byrow=TRUE,ncol=r,nrow=nper.rep)
        for (s in 1:nper)
        {          
        ## Loop over no. perms (r!)
            sig <- per[s,]
            sigpinv <- per.rep[,sig]
            psigpinv <- drop(1+rowSums((sigpinv-1)*dpow.mat))
            Sv <- Sv+v[psigpinv]
        }
    }    
    
    return(Sv/nper)
}
   
##########################################################################
## Sdrv.recursive computes the result of multiplying the symmetrizer matrix
## S_{d,r} by a vector v of length d^r, based on the recursive Algorithm 2
## as described in Section 4 of Chacon and Duong (2014)
##########################################################################

Sdrv.recursive <- function(d,r,v)
{
    if ((!is.vector(v)) & (!is.matrix(v))) { stop("v must be a vector or a matrix") }
    if (is.vector(v)) { v <- matrix(v,nrow=1) }
    n <- nrow(v)   
    if (ncol(v)!=d^r) { stop("Length of the vector(s) must equal d^r") }
    
    if ((r==0)|(r==1)) { w <- v } 
    else
    {
        per.rep <- pinv.all(d,r)
        nper.rep <- d^r
        dpow <- d^((r-1):0)
        dpow.mat <- matrix(rep(dpow,d^r),ncol=r,nrow=d^r,byrow=TRUE) 

        w <- v
        for (p in 2:r)
        {
            Tv <- matrix(0,nrow=n,ncol=d^r)
            for (j in 1:p)
            {
                per.rep.tau <- per.rep
                per.rep.tau[,j] <- per.rep[,p]
                per.rep.tau[,p] <- per.rep[,j]
                positions <- 1+rowSums((per.rep.tau-1)*dpow.mat)
                Tv[,positions] <- Tv[,positions]+w
            }
            w <- Tv/p
        }
    }

    return(drop(w))
}

##########################################################################
## Lp norm between two functions (grid based)
##########################################################################

Lpdiff <- function(f1, f2, p=2, index=1)
{
    if (is.vector(f1$H)) d <- 1 else d <- ncol(f1$H)
    f.diff <- f1
    if (is.list(f1$estimate)) f.diff$estimate <- (abs(f1$estimate[[index]] - f2$estimate[[index]]))^p
    else f.diff$estimate <- (abs(f1$estimate - f2$estimate))^p
    if (d==1)
    {
        delta <- diff(f.diff$eval.points)
        riemann.sum <- sum(c(delta[1], delta)*f.diff$estimate)
    }
    else if (d>1)
    {
        delta <- sapply(f.diff$eval.points, diff)
        delta <- rbind(head(delta, n=1), delta)
        if (d==2) riemann.sum <- sum(outer(delta[,1], delta[,2]) * f.diff$estimate)
        else if (d==3) riemann.sum <- sum(outer(outer(delta[,1], delta[,2]), delta[,3]) * f.diff$estimate)
    }  

    return(riemann.sum)
}

## ISE of difference between two KDEs

ise.diff <- function(fhat1, fhat2, xmin, xmax)
{
    if(!isTRUE(all.equal(fhat1$eval.points, fhat2$eval.points)))
    stop("fhat1 and fhat2 need to de defined on the same grid")

    fhat.sq <- fhat1
    fhat.sq$estimate <- (fhat1$estimate - fhat2$estimate)^2

    if (missing(xmin) & missing(xmax))
    int <- integral.kde(fhat=fhat.sq, q=max(fhat.sq$eval.points)+0.1*abs(max(fhat.sq$eval.points)), density=FALSE)

    if (missing(xmin) & !missing(xmax))
    int <- integral.kde(fhat=fhat.sq, q=xmax, density=FALSE)

    if (!missing(xmin) & missing(xmax))
    int <- integral.kde(fhat=fhat.sq, q=max(fhat.sq$eval.points)+0.1*abs(max(fhat.sq$eval.points)), density=FALSE) - integral.kde(fhat=fhat.sq, q=xmin, density=FALSE)

    if (!missing(xmin) & !missing(xmax))
    int <- integral.kde(fhat=fhat.sq, q=xmax, density=FALSE) - integral.kde(fhat=fhat.sq, q=xmin, density=FALSE)

    return(int)   
}

##########################################################################
## Create transparent colour
## partial alias for plot3D::alpha.col
##########################################################################

transparency.col <- function(col, alpha)
{ 
    trans.ind <- col!="transparent"
    col[trans.ind] <- plot3D::alpha.col(col[trans.ind], alpha=alpha)
    
    return(col)
}