समस्या सरल है। मेरे पास सहसंयोजक x और कुछ परिणाम y हैं और मैं x के आधार पर y का नादरिया-वाटसन अनुमान खोजना चाहता हूं। हालांकि, मैं एक ऐसा फ़ंक्शन ढूंढना चाहता हूं जो कई शर्तों को पूरा करता हो:

  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" गुठली का उपयोग करता है।

मैंने KernSmooth पैकेज से locpoly को भी आजमाया है:

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

इसके अलावा यह वज़न वापस नहीं करता है, मैं x के अपने विनिर्देश के लिए फ़ंक्शन चलाने में सक्षम नहीं था और मुझे कुछ निर्दिष्ट सीमा में समान दूरी वाले मानों का उपयोग करना पड़ा।

इसलिए मैं सोच रहा हूं कि क्या इन कार्यों में कुछ ऐसा है जो मुझे याद आ रहा है या क्या एनडब्ल्यू अनुमान के लिए आर में कोई अन्य समाधान है।

1
Adela 7 जून 2019, 15:06

2 जवाब

सबसे बढ़िया उत्तर

मैंने आपके उसी उदाहरण को Rcpp में कोडित किया है और यह R फ़ंक्शन की तुलना में तेज़ है लेकिन यह ksmooth से धीमा है। वैसे भी यह आपके इच्छित 2 तत्व लौटाता है। मैं कर्नेल को इनपुट के रूप में नहीं आने दे सकता क्योंकि Rcpp में ऐसा करना कठिन है जैसा आपने R में किया था, लेकिन आप Rcpp कोड के अंदर एक साधारण if else लिख सकते हैं जो आपके कर्नेल के आधार पर है। उपयोग करना चाहते हैं ([यहां] [1] आर में उपलब्ध वितरणों की एक सूची है)।

निम्नलिखित 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;
}

मैं 3 कार्यों की तुलना करने के लिए पैकेज microbenchmark का उपयोग करता हूं और परिणाम इस प्रकार है:

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
1
Alejandro Andrade 15 जून 2019, 19:53

यह locpol पैकेज का उपयोग करके किया जा सकता है जो सी ++ में हाथ से कार्यान्वयन की तुलना में बहुत तेजी से काम करता है:

library(locpol)
# weights
W <- locCteWeightsC(x = X, xeval = x, kernel = gaussK, bw = 0.5)$locWeig
# kernel estimate
m <- locWeightsEval(lpweig = W, y = Y)
1
Adela 8 जून 2020, 12:24