################################################################################# Estimate g_AMSE pilot bandwidths for even orders - 2-dim### Parameters## r - (r1, r2) partial derivative## n - sample size## psi1 - psi_(r + (2,0))## psi2 - psi_(r + (0,2))### Returns## g_AMSE pilot bandwidths for even orders###############################################################################gamse.even.2d <- function(r, n, psi1, psi2){  d <- 2  num <- -2 * dmvnorm.deriv(x=c(0,0), r=r, Sigma=diag(2))  den <- (psi1 + psi2) * n     g.amse <- (num/den)^(1/(2 + d + sum(r)))    return(g.amse)}################################################################################# Estimate g_AMSE pilot bandwidths for odd orders - 2-dim### Parameters### r - (r1, r2) partial derivative## n - sample size## psi1 - psi_(r + (2,0))## psi2 - psi_(r + (0,2))## psi00 - psi_(0,0)## RK - R(K^(r))## ## Returns## g_AMSE pilot bandwidths for odd orders###############################################################################gamse.odd.2d <- function(r, n, psi1, psi2, psi00, RK){      d <- 2    num <- 2 * psi00 * (2 * sum(r) + d) * RK    den <- (psi1 + psi2)^2 * n^2    g.amse <- (num/den)^(1/(2*sum(r) + d + 4))    return(g.amse)}################################################################################ Estimate g_SAMSE pilot bandwidth - 2- to 6-dim ## Parameters# Sigma.star - scaled variance matrix# n - sample size## Returns# g_SAMSE pilot bandwidth###############################################################################gsamse.nd <- function(Sigma.star, n, modr, nstage=1, psihat=NULL){  d <- ncol(Sigma.star)  K <- numeric(); psi <- numeric()  derivt4 <- deriv.list(d=d, r=4)  derivt6 <- deriv.list(d=d, r=6)   ## 4th order g_SAMSE  if (modr == 4)  {    for (i in 1:nrow(derivt4))    {      r <- derivt4[i,]      if (is.even(r))      {        K <- c(K, dmvnorm.deriv.diag(x=rep(0,d), r=r, Sigma=diag(d)))        A3psi <- 0        for (j in 1:d)        {          if (nstage==1)            A3psi <- A3psi + psins(r=r+2*elem(j,d), Sigma=Sigma.star)          else if (nstage==2)            A3psi <- A3psi + psihat[which.mat(r=r+2*elem(j,d), mat=derivt6)]        }        psi <- c(psi, A3psi)          }    }  }  ## 6th order g_SAMSE  else if (modr==6)  {    for (i in 1:nrow(derivt6))    {      r <- derivt6[i,]              if (is.even(r))      {        K <- c(K, dmvnorm.deriv.diag(x=rep(0,d), r=r, Sigma=diag(d)))               A3psi <- 0        for (j in 1:d)          A3psi <- A3psi + psins(r=r+2*elem(j,d), Sigma=Sigma.star)        psi <- c(psi, A3psi)      }    }  }    ## see thesis for formula  A1 <- sum(K^2)  A2 <- sum(K * psi)    A3 <- sum(psi^2)  B1 <- (2*modr + 2*d)*A1  B2 <- (modr + d - 2)*A2  B3 <- A3  gamma <- (-B2 + sqrt(B2^2 + 4*B1*B3)) / (2*B1)  g.samse <- (gamma * n)^(-1/(modr + d + 2))    return (g.samse)      }############################################################################### Estimate psi functionals for bivariate data using 1-stage plug-in - 2-dim## Parameters# x.star - pre-transformed data points# pilot - "amse" = different AMSE pilot bandwidths#       - "samse" = optimal SAMSE pilot bandwidth## Returns# estimated psi functionals###############################################################################psifun1.2d <- function(x.star, pilot="samse", binned, bin.par){  d <- 2  derivt4 <- deriv.list(d=d, r=4) ##cbind((2*d) - 0:(2*d), 0:(2*d))   S.star <- var(x.star)  n <- nrow(x.star)    RK31 <- 15/(64*pi)  psi00 <- psins(r=c(0,0), Sigma=S.star)   psihat.star <- vector()  g.star <- vector()    ## pilots are based on 4th order derivatives  ## compute 1 pilot for SAMSE  if (pilot=="samse")    g.star <- rep(gsamse.nd(S.star, n, 4), nrow(derivt4))   ## compute 5 different pilots for AMSE  else if (pilot=="amse")    for (k in 1:nrow(derivt4))    {       r <- derivt4[k,]      psi1 <- psins(r=r + 2*elem(1, 2), Sigma=S.star)      psi2 <- psins(r=r + 2*elem(2, 2), Sigma=S.star)      ## odd order      if (prod(r) == 3)        g.star[k] <- gamse.odd.2d(r, n, psi1, psi2, psi00, RK31)            ## even order      else        g.star[k] <- gamse.even.2d(r, n, psi1, psi2)    }    if (!binned)    x.star.diff <- differences(x.star, upper=FALSE)     for (k in 1:nrow(derivt4))  {    r <- derivt4[k,]    G.star <- g.star[k]^2 * diag(2)    if (binned)      psihat.star[k] <- kfe(bin.par=bin.par, G=G.star, r=r, binned=TRUE)    else       psihat.star[k] <- kfe.scalar(x=x.star.diff, r=r, g=g.star[k], diff=TRUE)   }   return(psihat.star)}################################################################################ Estimate psi functionals for bivariate data using 2-stage plug-in - 2-dim## Parameters# x - pre-transformed data points# pilot - "amse" - different AMSE pilot#       - "samse" - SAMSE pilot# Returns# estimated psi functionals###############################################################################psifun2.2d <- function(x.star, pilot="samse", binned, bin.par){   d <- 2  derivt4 <- deriv.list(d=d, r=4)  derivt6 <- deriv.list(d=d, r=6)  S.star <- var(x.star)  n <- nrow(x.star)    RK31 <- 15/(64*pi)  RK51 <- 945/(256*pi)  RK33 <- 225/(256*pi)  psi00 <- psins(r=c(0,0), Sigma=S.star)   psihat6.star <- vector()  g6.star <- vector()  psihat.star <- vector()  g.star <- vector()  ## pilots are based on 6th order derivatives    ## compute 1 pilot for SAMSE      if (pilot=="samse")    g6.star <- rep(gsamse.nd(S.star, n, 6), nrow(derivt6))        ## compute different pilots for AMSE  else if (pilot=="amse")  {           for (k in 1:nrow(derivt6))    {      r <- derivt6[k,]      psi1 <- psins(r=r + 2*elem(1, 2), Sigma=S.star)      psi2 <- psins(r=r + 2*elem(2, 2), Sigma=S.star)      if (prod(r) == 5)        g6.star[k] <- gamse.odd.2d(r, n, psi1, psi2, psi00, RK51)      else if (prod(r) == 9)        g6.star[k] <- gamse.odd.2d(r, n, psi1, psi2, psi00, RK33)       else          g6.star[k] <- gamse.even.2d(r, n, psi1, psi2)    }  }  if (!binned) x.star.diff <- differences(x.star, upper=FALSE)    for (k in 1:nrow(derivt6))  {    r <- derivt6[k,]    G6.star <- g6.star[k]^2 * diag(d)    if (binned)      psihat6.star[k] <- kfe(bin.par=bin.par, G=G6.star, r=r, binned=TRUE)    else      psihat6.star[k] <- kfe.scalar(x=x.star.diff, r=r, g=g6.star[k], diff=TRUE)  }   ## pilots are based on 4th order derivatives using 6th order psi functionals  ## computed above 'psihat6.star'      if (pilot=="samse")    g.star <- rep(gsamse.nd(S.star, n, 4, nstage=2, psihat=psihat6.star), nrow(derivt4))    else if (pilot=="amse")    for (k in 1:nrow(derivt4))    {      r <- derivt4[k,]      psi1 <- psihat6.star[7 - (r + 2*elem(1,2))[1]]      psi2 <- psihat6.star[7 - (r + 2*elem(2,2))[1]]            if (prod(r) == 3)        g.star[k] <- gamse.odd.2d(r, n, psi1, psi2, psi00, RK31)      else        g.star[k] <- gamse.even.2d(r, n, psi1, psi2)    }    for (k in 1:nrow(derivt4))  {    r <- derivt4[k,]    G.star <- g.star[k]^2 * diag(2)    if (binned)      psihat.star[k] <- kfe(bin.par=bin.par, G=G.star, r=r, binned=TRUE)    else       psihat.star[k] <- kfe.scalar(x=x.star.diff, r=r, g=g.star[k], diff=TRUE)  }  return(psihat.star)}################################################################################ Estimate psi functionals for 3-variate data using 1-stage plug-in - 3-dim## Parameters# x.star - pre-transformed data points# pilot - "samse" = optimal SAMSE pilot bandwidth# Returns# estimated psi functionals###############################################################################psifun1.nd <- function(x.star, d, pilot="samse", binned, bin.par){   derivt <- Psi4.list(d)$psi  derivt4 <- deriv.list(d, r=4)  S.star <- var(x.star)  n <- nrow(x.star)  psihat.star <- rep(0, length=nrow(derivt))  g.star <- vector()  if (!binned) x.star.diff <- differences(x.star, upper=FALSE)     ## compute 1 pilot for SAMSE  g.star <- gsamse.nd(S.star, n, 4, nstage=1)  G.star <- g.star^2 * diag(d)   for (k in 1:nrow(derivt4))  {    r <- derivt4[k,]    kind <- which.mat(r, derivt)    if (binned)      psihat.star[kind] <- kfe(bin.par=bin.par, G=G.star, r=r, binned=TRUE)    else       psihat.star[kind] <- kfe.scalar(x=x.star.diff, r=r, g=g.star, diff=TRUE)  }  return(psihat.star)}################################################################################ Estimate psi functionals for 3-variate data using 2-stage plug-in - 3-dim## Parameters# x.star - pre-transformed data points# pilot - "amse" = different AMSE pilot bandwidths#       - "samse" = optimal SAMSE pilot bandwidth## Returns# estimated psi functionals###############################################################################psifun2.nd <- function(x.star, d, pilot="samse", binned, bin.par){  derivt <- Psi4.list(d)$psi  derivt4 <- deriv.list(d, r=4)  derivt6 <- deriv.list(d, r=6)   S.star <- var(x.star)  n <- nrow(x.star)  if (!binned)    x.star.diff <- differences(x.star, upper=FALSE)  psihat6.star <- vector()  g6.star <- vector()  psihat.list.star <- vector()  psihat.star <- vector()  g.star <- vector()  ## pilots are based on 6th order derivatives     ## compute 1 pilot for SAMSE      if (pilot=="samse")    g6.star <- gsamse.nd(Sigma.star=S.star, n=n, modr=6)  G6.star <- g6.star^2 * diag(d)    for (k in 1:nrow(derivt6))  {    r <- derivt6[k,]           if (binned)      psihat6.star[k] <- kfe(bin.par=bin.par, G=G6.star, r=r, binned=TRUE)    else       psihat6.star[k] <- kfe.scalar(x=x.star.diff, r=r, g=g6.star, diff=TRUE)  }    ## pilots are based on 4th order derivatives using 6th order psi functionals  ## computed above 'psihat6.star'      if (pilot=="samse")    g.star <- gsamse.nd(S.star, n, 4, nstage=2, psihat=psihat6.star)   G.star <- g.star^2 * diag(d)    for (k in 1:nrow(derivt4))  {    r <- derivt4[k,]    kind <- which.mat(r, derivt)    if (binned)      psihat.star[kind] <- kfe(bin.par=bin.par, G=G.star, r=r, binned=TRUE)    else       psihat.star[kind] <- kfe.scalar(x=x.star.diff, r=r, g=g.star, diff=TRUE)  }    return(psihat.star)}############################################################################### Estimate psi functionals for 6-variate data using 1-stage plug-in ## with unconstrained pilot#### Parameters## x - data points## Sd4, Sd6 - symmetrizer matrices of order 4 and 6#### Returns## estimated psi functionals#############################################################################psifun1.unconstr.nd <- function(x, Sd4, Sd6, rel.tol=10^-10){  n <- nrow(x)  d <- ncol(x)  S <- var(x)    nlim <- 1e4  upper <- TRUE   ## stage 1 of plug-in  G4 <-(2^(d/2+3)/((d+4)*n))^(2/(d+8))*S  vecPsi4 <- vecPsir(x=x, Sdr=Sd4, Gr=G4, r=4, upper=upper, nlim=nlim)    return (vecPsi4)}############################################################################### Estimate psi functionals for 6-variate data using 2-stage plug-in ## with unconstrained pilot#### Parameters## x - data points## Sd4, Sd6 - symmetrizer matrices of order 4 and 6#### Returns## estimated psi functionals############################################################################psifun2.unconstr.nd <- function(x, Sd4, Sd6, rel.tol=10^-10){  d <- ncol(x)  n <- nrow(x)  S <- var(x)  Hstart <- (4/(d+2))^(2/(d+4))*n^(-2/(d+4))*S  Hstart <- matrix.sqrt(Hstart)  nlim <- 1e4   ## matrix of pairwise differences  upper <- TRUE  difs <- differences(x, upper=upper)   ## constants for normal reference  D4phi0 <- D4L0(d=d, Sd4=Sd4)  Id1 <- diag(d)  vId <- vec(Id1)  ## stage 1 of plug-in  G6 <- (2^(d/2+5)/((d+6)*n))^(2/(d+8))*S  G612 <- matrix.sqrt(G6)   vecPsi6 <- vecPsir(x=x, Sdr=Sd6, Gr=G6, r=6, upper=upper, nlim=nlim)     Id4 <- diag(d^4)  Id2 <- diag(d^2)  Kdd2 <- K.mat(m=d,n=d^2)  Psi6 <- (Id1%x%Kdd2%x%Id2)%*%vecPsi6    ## asymptotic squared bias for r = 4  AB2r4<-function(vechG){    r <- 4    G <- invvech(vechG)%*%invvech(vechG)    G12 <- matrix.sqrt(G)    Ginv12 <- chol2inv(chol(G12))    AB <- n^(-1)*det(Ginv12)*(K.pow(A=Ginv12,pow=r)%*%D4phi0)+      (1/2)*(t(vec(G))%x%Id4)%*%Psi6    return (sum(AB^2))  }    res <- optim(vech(Hstart),AB2r4, control=list(reltol=rel.tol))  V4 <- res$value  G4 <- res$par  G4 <- invvech(G4)%*%invvech(G4)   ## stage 2 of plug-in  vecPsi4 <- vecPsir(x=x, Sdr=Sd4, Gr=G4, r=4, upper=upper, nlim=nlim)    return (vecPsi4)}############################################################################## Psi_4 matrix of 4th order psi functionals used in AMISE - 2 to 6 dim#### Parameters## x - data points