在 R 中寻找理想的内核 NW 估计

Question

问题很简单。我有协变量 x 和一些结果 y，我想根据 x 找到 y 的 Nadarya-Watson 估计值。但是，我想找到一个满足几个条件的函数：

除了估计它return还有权重
它不仅处理提供估计值的均匀分布点。
相当快。

我可以简单地自己实现。我天真的估计函数看起来像这样：

mNW <- function(x, X, Y, h, K = dnorm) {

  # Arguments
  # x: evaluation points
  # X: covariates
  # Y: outcome
  # h: bandwidth
  # K: kernel

  Kx <- sapply(X, function(Xi) K((x - Xi) / h))

  # Weights
  W <- Kx / rowSums(Kx) 

  # NW estimate
  m <- W %*% Y

  return(list(m = m, W = W))
}

set.seed(123)
X <- rnorm(1000)
Y <- 1 + X - 2*X^2 + rnorm(1000)
x <- c(-3, -2.1, -0.7, 0, 0.3, 0.8, 1, 1.9, 3.2)

mNW(x, X, Y, h = 0.5)

它工作正常但速度很慢。所以我试图找到已经实现的东西。第一选择是 kernsmooth:

ksmooth(X, Y, kernel = "normal", bandwidth = 0.5, x.points = x)

这个速度更快，但是它没有 return 权重。此外，它仅使用 "box" 和 "normal" 内核。

我也试过 locpoly 来自 KernSmooth 包：

locpoly(X, Y, drv = 0, kernel = "normal", bandwidth = 0.5, 
        gridsize = 9, range.x = c(-3, 3.2))

除了它没有 return 权重之外，我无法运行为我自己的 x 规范发挥作用，我不得不在某些指定范围内使用等距值。

所以我想知道这些函数中是否缺少某些东西，或者 R 中是否有另一个用于 NW 估计的解决方案。

Answer 1

我在 Rcpp 中编写了相同的示例，它比 R 函数快得多，但比 ksmooth 慢。无论如何，它 returns 您想要的 2 个元素。我不能让内核作为输入，因为它在 Rcpp 中很难像在 R 中那样做，但你可以在 Rcpp 代码中编写一个简单的 if else，具体取决于你使用的内核想要使用（[此处][1] 是 R 中可用分布的列表）。

以下是应保存在 .cpp 文件中的 C++ 代码，并使用 Rcpp::sourceCpp()

源代码到 R

#include <RcppArmadillo.h>
using namespace Rcpp;
using namespace arma;

// [[Rcpp::depends(RcppArmadillo)]]

// [[Rcpp::export]]
std::vector<arma::mat>  mNWCpp(Rcpp::NumericVector x, Rcpp::NumericVector X,Rcpp::NumericMatrix Y,
           double h){

  int number_loop = X.size();
  int number_x    = x.size();

  Rcpp::NumericMatrix Kx(number_x,number_loop);

  for(int i =0; i<number_loop;++i){
    Kx(_,i) = dnorm((x-X[i])/h);
  }

  Rcpp::NumericVector row_sums = rowSums(Kx);
  Rcpp::NumericMatrix W = Kx*0;
  for(int i =0; i<number_loop;++i){
    W(_,i) = Kx(_,i)/row_sums;
  }


  arma::mat weights = Rcpp::as<arma::mat>(W);
  arma::mat Ymat = Rcpp::as<arma::mat>(Y);

  arma::mat m = weights * Ymat;

  std::vector< arma::mat> res(2);
  res[0] = weights;
  res[1] = m;
  return res;
}

我用包microbenchmark比较了3个函数，结果如下：

Unit: microseconds
 expr    min      lq     mean median      uq    max neval
    R 1991.9 2040.25 2117.656 2070.9 2123.50 3492.5   100
 Rcpp  490.5  502.10  512.318  510.8  517.35  598.0   100
   KS  196.8  205.40  215.598  211.4  219.15  282.2   100

Answer 2

这可以使用 locpol 包来完成，它比在 C++ 中手动实现快得多：

library(locpol)
# weights
W <- locCteWeightsC(x = X, xeval = x, kernel = gaussK, bw = 0.5)$locWeig
# kernel estimate
m <- locWeightsEval(lpweig = W, y = Y)

在 R 中寻找理想的内核 NW 估计

Looking for ideal kernel NW estimate in R

r

kernel-density