# whichbcv - 1 = BCV1#          - 2 = BCV2 # Hstart - initial bandwidth matrix## Returns# H_BCV###############################################################################Hbcv <- function(x, whichbcv=1, Hstart){  n <- nrow(x)  d <- ncol(x)  D2 <- rbind(c(1,0,0), c(0,1,0), c(0,1,0), c(0,0,1))  RK <- (4*pi)^(-d/2)  # use normal reference b/w matrix for bounds  k <- (((d+8)^((d+6)/2)*pi^(d/2)*RK)/(16*n*gamma((d+8)/2)*(d+2)))^(2/(d+4))  Hmax <- k * abs(var(x))  up.bound <- Hmax  ##lo.bound <- -Hmax  xdiff <- differences(x, upper=FALSE)  if (missing(Hstart))    Hstart <- matrix.sqrt(0.9*Hmax)  bin.par <- binning(x)  bcv1.mat.temp <- function(vechH)  {    H <- invvech(vechH) %*% invvech(vechH)    # ensures that H is positive definite    return(bcv.mat(x, H, H)$bcv)  }      bcv2.mat.temp <- function(vechH)  {    H <- invvech(vechH) %*% invvech(vechH)    return(bcv.mat(x, H, 2*H)$bcv)  }  # derivatives of BCV1 function - see thesis    bcv1.mat.deriv <- function(vechH)  {    H <-  invvech(vechH) %*% invvech(vechH)    Hinv <- chol2inv(chol(H))        psi4.mat <- bcv.mat(x, H, H)$psimat    psi22 <- psi4.mat[1,3]     psi00 <- dmvnorm.sum(x, Sigma=H, inc=0)/(n*(n-1))    psi22.deriv.xxt <- dmvnorm.deriv.2d.xxt.sum(x, r=c(2,2), Sigma=H)/(n*(n-1))    psi22.deriv <- t(D2)%*% vec((Hinv %*% psi22.deriv.xxt %*% Hinv + 2* psi00 *Hinv %*% Hinv - psi22*Hinv)/2)         const <- matrix(c(0,0,1, 0,4,0, 1,0,0), nc=3, byrow=TRUE)    psi4.mat.deriv<- const %x% psi22.deriv        deriv1 <- -1/2*n^{-1}*RK*t(D2) %*% vec(chol2inv(chol(H)))    deriv2 <- 1/2 * psi4.mat %*% vech(H) + 1/4 * (t(psi4.mat.deriv) %*% (vech(H) %x% diag(c(1,1,1)))) %*% vech(H)        return(deriv1 + deriv2)        }    # derivatives of BCV2 function - see thesis     bcv2.mat.deriv <- function(vechH)  {    H <-  invvech(vechH) %*% invvech(vechH)    Hinv <- chol2inv(chol(H))        psi4.mat <- bcv.mat(x, H, 2*H)$psimat    psi22 <- psi4.mat[1,3]     psi00 <- dmvnorm.2d.sum(x, Sigma=2*H, inc=0)/(n*(n-1))    psi22.deriv.xxt <- dmvnorm.deriv.2d.xxt.sum(x,r=c(2,2),Sigma=2*H)/(n*(n-1))    psi22.deriv <- t(D2)%*% vec((Hinv %*% psi22.deriv.xxt %*% Hinv +                                 2* psi00 *Hinv %*% Hinv - psi22*Hinv)/2)         const <- matrix(c(0,0,1, 0,4,0, 1,0,0), nc=3, byrow=TRUE)    psi4.mat.deriv<- const %x% psi22.deriv        deriv1 <- -1/2*n^{-1}*RK*t(D2) %*% vec(chol2inv(chol(H)))    deriv2 <- 1/2 * psi4.mat %*% vech(H) + 1/4 *      (t(psi4.mat.deriv) %*% (vech(H) %x% diag(c(1,1,1)))) %*% vech(H)        return(deriv1 + 2*deriv2)      }  if (whichbcv==1)    result <- optim(vech(Hstart), bcv1.mat.temp, gr=bcv1.mat.deriv,                    method="L-BFGS-B", upper=vech(matrix.sqrt(up.bound)),                    lower=-vech(matrix.sqrt(up.bound)))  else if (whichbcv==2)    result <- optim(vech(Hstart), bcv2.mat.temp, gr=bcv2.mat.deriv,                    method="L-BFGS-B", upper=vech(matrix.sqrt(up.bound)),                    lower=-vech(matrix.sqrt(up.bound)))    return(invvech(result$par) %*% invvech(result$par))}################################################################################ Find the diagonal bandwidth matrix that minimises the BCV for 2-dim# # Parameters# x - data values# whichbcv - 1 = BCV1#          - 2 = BCV2# Hstart - initial bandwidth matrix## Returns# H_BCV, diag###############################################################################Hbcv.diag <- function(x, whichbcv=1, Hstart){  n <- nrow(x)  d <- ncol(x)  ##D2 <- rbind(c(1,0,0), c(0,1,0), c(0,1,0), c(0,0,1))  RK <- (4*pi)^(-d/2)    ## use maximally smoothed b/w matrix for bounds  k <- (((d+8)^((d+6)/2)*pi^(d/2)*RK)/(16*n*gamma((d+8)/2)*(d+2)))^(2/(d+4))  Hmax <- k * abs(var(x))  up.bound <- diag(Hmax)  ##lo.bound <- rep(0,d)    if (missing(Hstart))    Hstart <- 0.9*matrix.sqrt(Hmax)  bcv1.mat.temp <- function(diagH)  {    H <- diag(diagH) %*% diag(diagH)    ## ensures that H is positive definite    return(bcv.mat(x, H, H)$bcv)  }      bcv2.mat.temp <- function(diagH)  {    H <- diag(diagH) %*% diag(diagH)    return(bcv.mat(x, H, 2*H)$bcv)  }    if (whichbcv == 1)    result <- optim(diag(Hstart), bcv1.mat.temp, method="L-BFGS-B", upper=sqrt(up.bound))  else if (whichbcv == 2)    result <- optim(diag(Hstart), bcv2.mat.temp, method="L-BFGS-B", upper=sqrt(up.bound))  return(diag(result$par) %*% diag(result$par))}########################################################################### Identifying elements of Theta_6 matrix########################################################################Theta6.elem <- function(d){  Theta6.mat <- list()  Theta6.mat[[d]] <- list()  for (i in 1:d)    Theta6.mat[[i]] <- list()    for (i in 1:d)    for (j in 1:d)    {        temp <- numeric()      for (k in 1:d)             for (ell in 1:d)              temp <- rbind(temp, elem(i,d)+2*elem(k,d)+2*elem(ell,d)+elem(j,d))            Theta6.mat[[i]][[j]] <- temp    }    return(Theta6.mat)}  ################################################################################ Estimate g_AMSE pilot bandwidth for SCV for 2 to 6 dim## Parameters# Sigma.star - scaled/ sphered variance matrix# Hamise - (estimate) of H_AMISE # n - sample size## Returns# g_AMSE pilot bandwidth###############################################################################gamse.scv.nd <- function(x.star, d, Sigma.star, Hamise, n, binned=FALSE, bin.par){  psihat6 <- vector()  g6.star <- gsamse.nd(Sigma.star, n, 6)   G6.star <- g6.star^2 * diag(d)  if (!binned) x.star.diff <- differences(x.star, upper=FALSE)  ##else bin.par <- binning(x=x.star, bgridsize=bgridsize, H=sqrt(diag(G6.star)))      derivt6 <- deriv.list(d=d, r=6)  for (k in 1:nrow(derivt6))  {    r <- derivt6[k,]    if (binned)      psihat6[k] <- kfe(bin.par=bin.par, G=G6.star, r=r, binned=TRUE)    else       psihat6[k] <- kfe.scalar(x=x.star.diff, r=r, g=g6.star, diff=TRUE)  }     Theta6.mat <- matrix(0, nc=d, nr=d)  Theta6.mat.ind <- Theta6.elem(d)  for (i in 1:d)    for (j in 1:d)    {      temp <- Theta6.mat.ind[[i]][[j]]      temp.sum <- 0      for (k in 1:nrow(temp))        temp.sum <- temp.sum + psihat6[which.mat(temp[k,], derivt6)]      Theta6.mat[i,j] <- temp.sum     }      eye3 <- diag(d)  D4 <- dupl(d)$d  trHamise <- tr(Hamise) ##[1,1] + Hamise[2,2] + Hamise[3,3]   ## required constants - see thesis  Cmu1 <- 1/2*t(D4) %*% vec(Theta6.mat %*% Hamise)  Cmu2 <- 1/8*(4*pi)^(-d/2) * (2*t(D4)%*% vec(Hamise) + trHamise * t(D4) %*% vec(eye3))  num <- 2 * (d+4) * sum(Cmu2*Cmu2)  den <- -(d+2) * sum(Cmu1*Cmu2) + sqrt((d+2)^2 * sum(Cmu1*Cmu2)^2 + 8*(d+4)*sum(Cmu1*Cmu1) * sum(Cmu2*Cmu2))  gamse <- (num / (den*n))^(1/(d+6))   return(gamse)}################################################################################ Computes the smoothed cross validation function for 2 to 6 dim# # Parameters# x - data values# H - bandwidth matrix# G - pilot bandwidth matrix## Returns# SCV(H)###############################################################################scv.1d <- function(x, h, g, binned=TRUE, bin.par, inc=1){  if (!missing(x)) n <- length(x)  if (!missing(bin.par)) n <- sum(bin.par$counts)  scv1 <- dnorm.sum(x=x, bin.par=bin.par, sigma=sqrt(2*h^2+2*g^2), binned=binned, inc=inc)  scv2 <- dnorm.sum(x=x, bin.par=bin.par, sigma=sqrt(h^2+2*g^2), binned=binned, inc=inc)  scv3 <- dnorm.sum(x=x, bin.par=bin.par, sigma=sqrt(2*g^2), binned=binned, inc=inc)  bias2 <-  (scv1 - 2*scv2 + scv3)  if (bias2 < 0) bias2 <- 0  scv <- (n*h)^(-1)*(4*pi)^(-1/2) + n^(-2)*bias2  return(scv)}scv.mat <- function(x, H, G, binned=FALSE, bin.par, diff=FALSE){  n <- nrow(x)  d <- ncol(x)  scv1 <- dmvnorm.sum(x=x, Sigma=2*H + 2*G, inc=1, bin.par=bin.par, binned=binned, diff=diff)  scv2 <- dmvnorm.sum(x=x, Sigma=H + 2*G, inc=1, bin.par=bin.par, binned=binned, diff=diff)  scv3 <- dmvnorm.sum(x=x, Sigma=2*G, inc=1, bin.par=bin.par, binned=binned, diff=diff)  scvmat <- n^(-1)*det(H)^(-1/2)*(4*pi)^(-d/2) + n^(-2)*(scv1 - 2*scv2 + scv3)      return (scvmat)}################################################################################ Find the bandwidth that minimises the SCV for 1 to 6 dim# # Parameters# x - data values# pre - "scale" - pre-scaled data#     - "sphere"- pre-sphered data# Hstart - initial bandwidth matrix## Returns# H_SCV###############################################################################hscv <- function(x, nstage=2, binned=TRUE, bgridsize, plot=FALSE){  sigma <- sd(x)  n <- length(x)  d <- 1  hnorm <- sqrt((4/(n*(d + 2)))^(2/(d + 4)) * var(x))  if (missing(bgridsize)) bgridsize <- 401  ##if (missing(hmin))  hmin <- 0.1*hnorm  ##if (missing(hmax))  hmax <- 2*hnorm  ##bin.par.sub <- binning(x=x[1:min(n, 1e4)], bgridsize=bgridsize, h=hnorm)  bin.par <- binning(x=x, bgridsize=bgridsize, h=hnorm)  if (nstage==1)  {    psihat6 <- psins.1d(r=6, sigma=sigma)    psihat10 <- psins.1d(r=10, sigma=sigma)  }  else if (nstage==2)  {    ##psihat8 <- psins.1d(r=8, sigma=sigma)    ##psihat12 <- psins.1d(r=12, sigma=sigma)    g1 <- (2/(7*n))^(1/9)*2^(1/2)*sigma    g2 <- (2/(11*n))^(1/13)*2^(1/2)*sigma    psihat6 <- kfe.1d(bin.par=bin.par, binned=TRUE, r=6, g=g1, inc=1)    psihat10 <- kfe.1d(bin.par=bin.par, binned=TRUE, r=10, g=g2, inc=1)  }  g3 <- (-6/((2*pi)^(1/2)*psihat6*n))^(1/7)   g4 <- (-210/((2*pi)^(1/2)*psihat10*n))^(1/11)  psihat4 <- kfe.1d(bin.par=bin.par, binned=TRUE, r=4, g=g3, inc=1)  psihat8 <- kfe.1d(bin.par=bin.par, binned=TRUE, r=8, g=g4, inc=1)  C <- (441/(64*pi))^(1/18) * (4*pi)^(-1/5) * psihat4^(-2/5) * psihat8^(-1/9)    scv.1d.temp <- function(h)  {    return(scv.1d(x=x, bin.par=bin.par, h=h, g=C*n^(-23/45)*h^(-2), binned=binned, inc=1))  }  if (plot)  {      hseq <- seq(hmin,hmax, length=400)    hscv.seq <- rep(0, length=length(hseq))    for (i in 1:length(hseq))      hscv.seq[i] <- scv.1d.temp(hseq[i])    plot(hseq, hscv.seq, type="l", xlab="h", ylab="SCV(h)")  }    opt <- optimise(f=scv.1d.temp, interval=c(hmin, hmax))$minimum  if (n >= 1e5) warning("hscv is not always stable for large samples")    return(opt)}Hscv <- function(x, pre="sphere", Hstart, binned=FALSE, bgridsize){  d <- ncol(x)  RK <- (4*pi)^(-d/2)  if (substr(pre,1,2)=="sc")    pre <- "scale"  else if (substr(pre,1,2)=="sp")    pre <- "sphere"  if(!is.matrix(x)) x <- as.matrix(x)  ## pre-transform data  if (pre=="sphere")    x.star <- pre.sphere(x)  else if (pre=="scale")    x.star <- pre.scale(x)  S.star <- var(x.star)  n <- nrow(x.star)  if (n > 1000 & !binned)    warning("Hscv converges slowly for n > 1000 without binned estimation")    if (missing(bgridsize) & binned) bgridsize <- default.gridsize(d)    if (d > 4) binned <- FALSE   if (pre=="scale") S12 <- diag(sqrt(diag(var(x))))  else if (pre=="sphere") S12 <- matrix.sqrt(var(x))    S12inv <- chol2inv(chol(S12))  Hamise <- S12inv %*% Hpi(x=x, nstage=1, pilot="samse", pre="sphere", binned=TRUE, bgridsize=bgridsize) %*% S12inv  if (any(is.na(Hamise)))  {    warning("Pilot bandwidth matrix is NA - replaced with maximally smoothed")    Hamise <- (((d+8)^((d+6)/2)*pi^(d/2)*RK)/(16*(d+2)*n*gamma(d/2+4)))^(2/(d+4))* var(x.star)  }   if (binned)  {    bin.par <- binning(x=x.star, bgridsize=bgridsize)    gamse <- gamse.scv.nd(x.star=x.star, d=d, Sigma.star=S.star, H=Hamise, n=n, binned=TRUE, bin.par=bin.par)  }  else  {    x.star.diff <- differences(x.star, upper=FALSE)    gamse <- gamse.scv.nd(x.star=x.star, d=d, Sigma.star=S.star, H=Hamise, n=n, binned=FALSE)  }  G.amse <- gamse^2 * diag(d)    ## use normal reference bandwidth as initial condition  if (missing(Hstart))     Hstart <- (4/(n*(d + 2)))^(2/(d + 4)) * var(x.star)  else        Hstart <- S12inv %*% Hstart %*% S12inv  Hstart <- matrix.sqrt(Hstart)  scv.mat.temp <- function(vechH)  {    H <- invvech(vechH) %*% invvech(vechH)    return(scv.mat(x.star, H, G.amse))  }    ## back-transform  result <- optim(vech(Hstart), scv.mat.temp, method= "Nelder-Mead")                                        #control=list(abstol=n^(-10*d)))  H <- invvech(result$par) %*% invvech(result$par)  H <- S12 %*% H %*% S12   return(H)}Hscv.diag <- function(x, pre="scale", Hstart, binned=FALSE, bgridsize){  if(!is.matrix(x)) x <- as.matrix(x)  d <- ncol(x)  RK <- (4*pi)^(-d/2)    ## pre-transform data    if (substr(pre,1,2)=="sc")    pre <- "scale"  else if (substr(pre,1,2)=="sp")    pre <- "sphere"  if (pre=="sphere")     stop("Using pre-sphering doesn't give a diagonal bandwidth matrix\n")    if (pre=="sphere")    x.star <- pre.sphere(x)  else if (pre=="scale")    x.star <- pre.scale(x)    S.star <- var(x.star)  n <- nrow(x.star)  if (missing(bgridsize) & binned) bgridsize <- default.gridsize(d)  if (d > 4) binned <- FALSE     if (pre=="scale") S12 <- diag(sqrt(diag(var(x))))  else if (pre=="sphere") S12 <- matrix.sqrt(var(x))  S12inv <- chol2inv(chol(S12))  Hamise <- S12inv %*% Hpi.diag(x=x,nstage=1, pilot="samse", pre="scale", binned=binned, bgridsize=bgridsize) %*% S12inv  if (any(is.na(Hamise)))  {    warning("Pilot bandwidth matrix is NA - replaced with maximally smoothed")    Hamise <- (((d+8)^((d+6)/2)*pi^(d/2)*RK)/(16*(d+2)*n*gamma(d/2+4)))^(2/(d+4))* var(x.star)  }   if (binned)  {      bin.par <- binning(x=x.star, bgridsize=bgridsize, H=diag(diag(Hamise)))        gamse <- gamse.scv.nd(x.star=x.star, d=d, Sigma.star=S.star, H=Hamise, n=n, binned=binned, bin.par=bin.par)  }  else    gamse <- gamse.scv.nd(x.star=x.star, d=d, Sigma.star=S.star, H=Hamise, n=n, binned=FALSE)  G.amse <- gamse^2 * diag(d)    ## use normal reference bandwidth as initial condition  if (missing(Hstart))     Hstart <- (4/(n*(d + 2)))^(2/(d + 4)) * var(x.star)  else        Hstart <- S12inv %*% Hstart %*% S12inv  Hstart <- matrix.sqrt(Hstart)  scv.mat.temp <- function(diagH)  {    ## ensures that H is positive definite    H <- diag(diagH) %*% diag(diagH)    return(scv.mat(x.star, H, G.amse, binned=binned, bin.par=bin.par))  }    ## back-transform  result <- optim(diag(Hstart), scv.mat.temp, method= "Nelder-Mead")  H <- diag(result$par) %*% diag(result$par)  H <- S12 %*% H %*% S12    return(H)}