normal.r

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

      d <- ncol(x)
      difs <- differences(x, upper=FALSE)
    }
    else
    {
      n <- sqrt(nrow(x)) #(-1 + sqrt(1+8*nrow(x)))/2
      d <- ncol(x)
      difs <- x
    }
    if (is.diagonal(Sigma))
      sumval <- sum(dmvnorm.deriv.diag(difs, mu=rep(0,d), Sigma=Sigma, r=r))
    else
      sumval <- sum(dmvnorm.deriv(difs, mu=rep(0,d), Sigma=Sigma, r=r))
  }

  if (inc==0)
    sumval <- sumval - n*dmvnorm.deriv(x=rep(0,d), mu=rep(0,d), Sigma=Sigma, r=r)
  if (kfe)
    if (inc==1) sumval <- sumval/n^2
    else sumval <- sumval/(n*(n-1))
  
  return(sumval)
}









###############################################################################
## Partial derivatives of the bivariate normal (mean 0) 
## 
## Parameters
## x - points to evaluate at
## Sigma - variance matrix
## r - (r1, r2) vector of partial derivative indices 
#
## Returns
## r-th partial derivative at x
###############################################################################

dmvnorm.deriv.2d <- function(x, Sigma, r) 
{
  if (is.vector(x))
  {    
    n <- 1; x1 <- x[1]; x2 <- x[2]
  }
  else 
  {
    n <- nrow(x); x1 <- x[,1]; x2 <- x[,2]
  }

  if (sum(r) == 0)
    return(dmvnorm(x, c(0,0), Sigma))
  ## first order derivatives  
  else if (sum(r) == 1)
    derivt <- .C("dmvnormd1_2d", as.double(x1), as.double(x2), 
                 as.double(vec(Sigma)), as.integer(r), as.integer(n), 
                 as.double(rep(0, n)), PACKAGE="ks")
  ## third order derivatives
  else if (sum(r) == 2)
    derivt <- .C("dmvnormd2_2d", as.double(x1), as.double(x2), 
                 as.double(vec(Sigma)), as.integer(r), as.integer(n), 
                 as.double(rep(0, n)), PACKAGE="ks")
  ## second order derivatives
  else if (sum(r) == 3)
    derivt <- .C("dmvnormd3_2d", as.double(x1), as.double(x2), 
                 as.double(vec(Sigma)), as.integer(r), as.integer(n), 
                 as.double(rep(0, n)), PACKAGE="ks")
  ## fourth order derivatives
  else if (sum(r) == 4)
    derivt <- .C("dmvnormd4_2d", as.double(x1), as.double(x2), 
                 as.double(vec(Sigma)), as.integer(r), as.integer(n), 
                 as.double(rep(0, n)), PACKAGE="ks")
  ## fifth order derivatives
  else if (sum(r) == 5)
    derivt <- .C("dmvnormd5_2d", as.double(x1), as.double(x2), 
                 as.double(vec(Sigma)), as.integer(r), as.integer(n), 
                 as.double(rep(0, n)), PACKAGE="ks")
  ## sixth order derivatives
  else if (sum(r) == 6)
    derivt <- .C("dmvnormd6_2d", as.double(x1), as.double(x2), 
            as.double(vec(Sigma)), as.integer(r), as.integer(n), 
                 as.double(rep(0, n)), PACKAGE="ks")
  else
    stop("Only works for up to 6th order partial derivatives")
  
  return(derivt[[6]])
}         

###############################################################################
## Partial derivatives of the 3-variate normal (mean 0, diag variance matrix) 
##############################################################################

dmvnorm.deriv.3d <- function(x, Sigma, r) 
{
  ####### Diagonal variance matrix implemented ONLY at the moment
  
  ##d <- 3
  if (sum(r) > 6)
    stop("Only works for up to 6th order partial derivatives")

  if (is.vector(x))
  {    
    ##n <- 1;
    x1 <- x[1]; x2 <- x[2]; x3 <- x[3]; 
  }
  else 
  {
    ##n <- nrow(x);
    x1 <- x[,1]; x2 <- x[,2]; x3 <- x[,3]; 
  }
  
  y1 <- cbind(x1,x2)
  y2 <- x3
  r1 <- r[c(1,2)]
  r2 <- r[3]
  Sigma1 <- diag(c(Sigma[1,1], Sigma[2,2]))
  Sigma2 <- Sigma[3,3]

  ## Use existing 2-dim derivatives to compute 3-dim derivative for diag
  ## variance matrix
  derivt1 <- dmvnorm.deriv.2d(y1, Sigma1, r1)
  derivt2 <- dnorm.deriv(x=y2, mu=0, sigma=sqrt(Sigma2), r=r2)
  derivt <- derivt1*derivt2

  return(derivt)
}         

##############################################################################
## Partial derivatives of the 4-variate normal (mean 0, diag variance matrix) 
##############################################################################

dmvnorm.deriv.4d <- function(x, Sigma, r) 
{
  ####### Diagonal variance matrix implemented ONLY at the moment
  
  ##d <- 4
  if (sum(r) > 6)
    stop("Only works for up to 6th order partial derivatives")

  if (is.vector(x))
  {    
    ##n <- 1;
    x1 <- x[1]; x2 <- x[2]; x3 <- x[3]; x4 <- x[4];
  }
  else 
  {
    ##n <- nrow(x);
    x1 <- x[,1]; x2 <- x[,2]; x3 <- x[,3]; x4 <- x[,4];
  }
  
  y1 <- cbind(x1,x2)
  y2 <- cbind(x3,x4)
  r1 <- r[c(1,2)]
  r2 <- r[c(3,4)]
  Sigma1 <- diag(c(Sigma[1,1], Sigma[2,2]))
  Sigma2 <- diag(c(Sigma[3,3], Sigma[4,4]))

  ## Use existing 2-dim derivatives to compute 4-dim derivative for diag
  ## variance matrix
  derivt1 <- dmvnorm.deriv.2d(y1, Sigma1, r1)
  derivt2 <- dmvnorm.deriv.2d(y2, Sigma2, r2)
  derivt <- derivt1*derivt2

  return(derivt)
}
###############################################################################
## Partial derivatives of the 5-variate normal (mean 0, diag variance matrix) 
## 
## Parameters
## x - points to evaluate at
## Sigma - variance matrix
## r - vector of partial derivative indices 
#
## Returns
## r-th partial derivative at x
##############################################################################

dmvnorm.deriv.5d <- function(x, Sigma, r) 
{
  ####### Diagonal variance matrix implemented ONLY at the moment

  ##d <- 5
  if (sum(r) > 6)
    stop("Only works for 2nd, 4th and 6th order partial derivatives")
 
  if (is.vector(x))
  {    
    ##n <- 1;
    x1 <- x[1]; x2 <- x[2]; x3 <- x[3]; x4 <- x[4]; x5 <- x[5]; 
  }
  else 
  {
    ##n <- nrow(x);
    x1 <- x[,1]; x2 <- x[,2]; x3 <- x[,3]; x4 <- x[,4]; x5 <- x[,5]; 
  }
  
  y1 <- cbind(x1,x2)
  y2 <- cbind(x3,x4)
  y3 <- x5
  r1 <- r[c(1,2)]
  r2 <- r[c(3,4)]
  r3 <- r[5]
  Sigma1 <- diag(c(Sigma[1,1], Sigma[2,2]))
  Sigma2 <- diag(c(Sigma[3,3], Sigma[4,4]))
  Sigma3 <- Sigma[5,5]
  
  ## Use existing 2-dim derivatives to compute 5-dim derivative for diag
  ## variance matrix
  derivt1 <- dmvnorm.deriv.2d(y1, Sigma1, r1)
  derivt2 <- dmvnorm.deriv.2d(y2, Sigma2, r2)
  derivt3 <- dnorm.deriv(x=y3, mu=0, sigma=sqrt(Sigma3), r=r3)
  derivt <- derivt1*derivt2*derivt3

  return(derivt)
}         


###############################################################################
## Partial derivatives of the 6-variate normal (mean 0, diag variance matrix) 
## 
## Parameters
## x - points to evaluate at
## Sigma - variance matrix
## r - vector of partial derivative indices 
#
## Returns
## r-th partial derivative at x
##############################################################################

dmvnorm.deriv.6d <- function(x, Sigma, r) 
{
  ####### Diagonal variance matrix implemented ONLY at the moment

  ##d <- 6
  if (sum(r) > 6)
    stop("Only works for 2nd, 4th and 6th order partial derivatives")

  if (is.vector(x))
  {    
    ##n <- 1;
    x1 <- x[1]; x2 <- x[2]; x3 <- x[3]; x4 <- x[4]; x5 <- x[5]; x6 <- x[6];
  }
  else 
  {
    ##n <- nrow(x);
    x1 <- x[,1]; x2 <- x[,2]; x3 <- x[,3]; x4 <- x[,4]; x5 <- x[,5]; x6 <- x[,6];
  }
  
  y1 <- cbind(x1,x2)
  y2 <- cbind(x3,x4)
  y3 <- cbind(x5,x6)
  r1 <- r[c(1,2)]
  r2 <- r[c(3,4)]
  r3 <- r[c(5,6)]
  Sigma1 <- diag(c(Sigma[1,1], Sigma[2,2]))
  Sigma2 <- diag(c(Sigma[3,3], Sigma[4,4]))
  Sigma3 <- diag(c(Sigma[5,5], Sigma[6,6]))
  
  ## Use existing 2-dim derivatives to compute 6-dim derivative for diag
  ## variance matrix
  derivt1 <- dmvnorm.deriv.2d(y1, Sigma1, r1)
  derivt2 <- dmvnorm.deriv.2d(y2, Sigma2, r2)
  derivt3 <- dmvnorm.deriv.2d(y3, Sigma3, r3)
  derivt <- derivt1*derivt2*derivt3

  return(derivt)
}         

###############################################################################
## Double sum  of K(X_i - X_j) used in density derivative estimation
#
## Parameters
## x - points to evaluate
## Sigma - variance matrix
## inc - 0 - exclude diagonals
##     - 1 - include diagonals
#
## Returns
## Double sum at x
###############################################################################

dmvnorm.2d.sum <- function(x, Sigma, inc=1)
{
  if (is.vector(x))
  {
    n <- 1; d <- 1; x1 <- x[1]; x2 <- x[2]
  }
  else
  {
    n <- nrow(x); d <- 2; x1 <- x[,1]; x2 <- x[,2]
  }
  
  viSigma <- vec(chol2inv(chol(Sigma)))
  result <- .C("dmvnorm_2d_sum", as.double(x1), as.double(x2),
               as.double(viSigma), as.double(det(Sigma)), as.integer(n),
               as.double(0), PACKAGE="ks")
  sumval <- result[[6]]
  
  ## above C function mvnorm_2d_sum only computes the upper triangular half
  ## so need to reflect along the diagonal and then subtract appropriate
  ## amount to compute whole sum 
  
  if (inc == 0) 
    sumval <- 2*sumval - 2*n*dmvnorm(c(0,0), mean=c(0,0), sigma=Sigma)
  else if (inc == 1)
    sumval <- 2*sumval - n*dmvnorm(c(0,0), mean=c(0,0), sigma=Sigma) 
  
  return(sumval)
}

dmvnorm.deriv.2d.sum <- function(x, Sigma, r, inc=1)
{
  if (is.vector(x))
  {
    n <- 1; x1 <- x[1]; x2 <- x[2]
  }
  else
  {
    n <- nrow(x); x1 <- x[,1]; x2 <- x[,2]
  }
     
  if (sum(r)==0)
    return(dmvnorm.2d.sum(x, Sigma, inc=inc))
  ## 1st order derivatives
  else if (sum(r)==1)
    derivt <- .C("dmvnormd1_2d_sum", as.double(x1), as.double(x2), 
                 as.double(vec(Sigma)), as.integer(r), as.integer(n), as.double(0), PACKAGE="ks")
  ## 2nd order derivativeselse if (sum(r)==2)
  else if (sum(r)==2)
    derivt <- .C("dmvnormd2_2d_sum", as.double(x1), as.double(x2), 
                 as.double(vec(Sigma)), as.integer(r), as.integer(n), as.double(0), PACKAGE="ks")
  ## 3rd order derivativeselse if (sum(r)==3)
  else if (sum(r)==3)
    derivt <- .C("dmvnormd3_2d_sum", as.double(x1), as.double(x2), 
                 as.double(vec(Sigma)), as.integer(r), as.integer(n), as.double(0), PACKAGE="ks")
  ## fourth order derivatives
  else if (sum(r) == 4)
    derivt <- .C("dmvnormd4_2d_sum", as.double(x1), as.double(x2), 
                   as.double(vec(Sigma)), as.integer(r), as.integer(n),
                 as.double(0), PACKAGE="ks")
  ## fifth order derivatives
  else if (sum(r) == 5)
    derivt <- .C("dmvnormd5_2d_sum", as.double(x1), as.double(x2), 
                 as.double(vec(Sigma)), as.integer(r), as.integer(n),
                 as.double(0), PACKAGE="ks")
  ## sixth order derivatives
  else if (sum(r) == 6)
    derivt <- .C("dmvnormd6_2d_sum", as.double(x1), as.double(x2), 
                 as.double(vec(Sigma)), as.integer(r), as.integer(n),
                   as.double(0), PACKAGE="ks")
  else
    stop("Only works for up to 6th order partial derivatives")
  
  if (sum(r)==0)
    sumval <- derivt
    else
      sumval <- derivt[[6]]
  
  ## above C functions mvnorm?_2d_sum only computes the upper triangular half
  ## so need to reflect along the diagonal and then subtract appropriate
  ## amount to compute whole sum 
  if (inc == 0) 
    sumval <- 2*sumval - 2*n*dmvnorm.deriv.2d(x=c(0,0), r=r, Sigma=Sigma)
  else if (inc == 1)
    sumval <- 2*sumval - n*dmvnorm.deriv.2d(x=c(0,0), r=r, Sigma=Sigma) 

  return(sumval)
}         


dmvnorm.deriv.2d.xxt.sum <- function(x, Sigma, r)
{
  if (is.vector(x))
  {
    n <- 1; x1 <- x[1]; x2 <- x[2]
  }
  else
  {
    n <- nrow(x); x1 <- x[,1]; x2 <- x[,2]
  }
  viSigma <- vec(chol2inv(chol(Sigma)))
  result <- .C("dmvnormd4_2d_xxt_sum", as.double(x1), as.double(x2),
               as.double(viSigma), as.integer(r), as.integer(n), 
               as.double(rep(0,4)), PACKAGE="ks")
  
  ## above C functions dmvnorm4_2d_xxt_sum only computes the upper triangular
  ## half so need to reflect along the diagonal to compute whole sum
  sumval <- 2*result[[6]]
  
  return(invvec(sumval))
} 


###############################################################################
## Double sum  of K(X_i - X_j) used in density derivative estimation - 3-dim
#
## Parameters
## x - points to evaluate
## Sigma - variance matrix
## inc - 0 - exclude diagonals
##     - 1 - include diagonals
#
## Returns
## Double sum at x
###############################################################################

dmvnorm.3d.sum <- function(x, Sigma, inc=1)
{
  if (is.vector(x))
  {
    n <- 1; d <- 4; x1 <- x[1]; x2 <- x[2]; x3 <- x[3]; 
  }
  else
  {
    n <- nrow(x); d <- ncol(x); x1 <- x[,1]; x2 <- x[,2]; x3 <- x[,3];
  }
  
  viSigma <- vec(chol2inv(chol(Sigma)))
  result <- .C("dmvnorm_3d_sum", as.double(x1),as.double(x2), as.double(x3), 
               as.double(viSigma), as.double(det(Sigma)), as.integer(n),
               as.double(0), PACKAGE="ks")
  sumval <- result[[7]]
  
  ## above C function mvnorm_3d_sum only computes the upper triangular half
  ## so need to reflect along the diagonal and then subtract appropriate
  ## amount to compute whole sum 
  
  if (inc == 0) 
    sumval <- 2*sumval - 2*n*dmvnorm(rep(0,d), rep(0,d), sigma=Sigma)
  else if (inc == 1)
    sumval <- 2*sumval - n*dmvnorm(rep(0,d), rep(0,d), sigma=Sigma) 
  
  return(sumval)
}

dmvnorm.deriv.3d.sum <- function(x, Sigma, r, inc=1, binned=FALSE, bin.par)
{
  if (is.vector(x))
  {
    n <- 1; x1 <- x[1]; x2 <- x[2] ; x3 <- x[3];  
  }
  else
  {
    n <- nrow(x); x1 <- x[,1]; x2 <- x[,2]; x3 <- x[,3];   
  }

  d <- 3
  sumval <- 0

  if (binned)
  {
    ## drvkde computes include-diagonals estimate
    fhatr <- drvkde(x=bin.par$counts, drv=r, bandwidth=sqrt(diag(Sigma)),
                       binned=TRUE, range.x=bin.par$range.x, se=FALSE)$est
    sumval <- sum(bin.par$counts * n * fhatr)
    if (inc == 0) 
      sumval <- sumval - n*dmvnorm.deriv.3d(x=rep(0,d), r=r, Sigma=Sigma)
  }
  else
  {  
    for (j in 1:n)
    {  
      y1 <- x1 - x1[j]
      y2 <- x2 - x2[j]
      y3 <- x3 - x3[j]
      sumval <- sumval + sum(dmvnorm.deriv.3d(cbind(y1, y2, y3), Sigma, r))
    }
    
    if (inc==0)
      sumval <- sumval - n*dmvnorm.deriv.3d(rep(0,d), Sigma, r)
  }

  
  return(sumval)
}

###############################################################################
# Double sum  of K(X_i - X_j) used in density derivative estimation - 4-dim
#
## Parameters
## x - points to evaluate
## Sigma - variance matrix
## inc - 0 - exclude diagonals
##     - 1 - include diagonals
#
## Returns
## Double sum at x
###############################################################################

dmvnorm.4d.sum <- function(x, Sigma, inc=1)
{
  if (is.vector(x))
  {
    n <- 1; d <- 4; x1 <- x[1]; x2 <- x[2]; x3 <- x[3]; x4 <- x[4]; 
  }
  else
  {
    n <- nrow(x); d <- ncol(x); x1 <- x[,1]; x2 <- x[,2]; x3 <- x[,3]; x4 <- x[,4];
  }

  viSigma <- vec(chol2inv(chol(Sigma)))
  result <- .C("dmvnorm_4d_sum", as.double(x1), as.double(x2),
               as.double(x3), as.double(x4),
               as.double(viSigma), as.double(det(Sigma)), as.integer(n),
               as.double(0), PACKAGE="ks")
  sumval <- result[[8]]
    
  ## above C function mvnorm_2d_sum only computes the upper triangular half
  ## so need to reflect along the diagonal and then subtract appropriate
  ## amount to compute whole sum 
  
  if (inc == 0) 
    sumval <- 2*sumval - 2*n*dmvnorm(rep(0,d), rep(0,d), sigma=Sigma)
  else if (inc == 1)
    sumval <- 2*sumval - n*dmvnorm(rep(0,d), rep(0,d), sigma=Sigma) 

  return(sumval)
}

dmvnorm.deriv.4d.sum <- function(x, Sigma, r, inc=1, binned=FALSE, bin.par)
{