Ищем идеальную оценку ЯО в R - PullRequest
Новая
       20

Ищем идеальную оценку ЯО в R

1 голос
/ 07 июня 2019

Проблема проста.У меня есть ковариата x и некоторый результат y, и я хотел бы найти оценку Надарьи-Ватсона y на основе x.Однако я хотел бы найти функцию, которая удовлетворяет нескольким условиям:

  1. Помимо оценки, она возвращает также веса
  2. Она обрабатывает не только равномерно распределенные точки, для которых предоставляется оценка.
  3. Это достаточно быстро.

Я могу просто реализовать это самостоятельно.Моя наивная функция оценки выглядит примерно так: 

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)

Этот вариант быстрее, но он не возвращает веса.Более того, он использует только ядра "box" и "normal".

Я также пробовал locpoly из KernSmooth пакета: 

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

Кроме того, он не возвращает веса, я не смог запустить функцию для моей собственной спецификации xи мне пришлось использовать одинаково расположенные значения в некотором указанном диапазоне.

Поэтому мне интересно, есть ли что-то, чего мне не хватало в этих функциях, или есть другое решение в R для оценки NW.

1 Ответ

1 голос
/ 15 июня 2019

Я кодировал ваш тот же пример в Rcpp, и он намного быстрее, чем функция R, но медленнее, чем ksmooth.В любом случае он возвращает 2 элемента, которые вы хотели.Я не мог использовать ядро ​​в качестве входных данных, потому что это трудно сделать в Rcpp, как вы это делали в R, но вы можете написать простой if else внутри кода Rcpp в зависимости от того, какие ядра вы хотите использовать ([здесь] [1] - список доступных дистрибутивов в R).

Ниже приведен код C ++, который должен быть сохранен в файле .cpp с исходным кодом в R с Rcpp::sourceCpp()

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