selector.r

r软件 另一款可以计算核估计的软件包 需安装r软件

###############################################################################
## 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