वैक्टरों पर सामान्य वितरण के सीडीएफ की गणना करने का सबसे तेज़ तरीका - आर :: pnorm बनाम erfc बनाम?

Question

मुझे आशा है कि मेरा पुनर्निर्मित प्रश्न अब स्टैक ओवरफ्लो के मानदंडों को फिट करेगा। कृपया नीचे दिए गए उदाहरण पर विचार करें। मैं एक लॉग-लिकेलहुड फ़ंक्शन लिख रहा हूं जिसमें सीडीएफ पर वैक्टर की गणना करना सबसे अधिक समय लेने वाला हिस्सा है। उदाहरण 1 R::pnorm का उपयोग करता है, उदाहरण 2 सामान्य सीडीएफ को erfc के साथ अनुमानित करता है। जैसा कि आप देख सकते हैं कि परिणाम पर्याप्त रूप से समान हैं, ercf संस्करण थोड़ा तेज़ है।

अभ्यास में (एक एमएलई के भीतर) हालांकि यह पता चला है कि ercf सटीक नहीं है, जो एल्गोरिदम को inf क्षेत्रों में चलाता है जब तक कोई बाधाओं को सटीक रूप से सेट नहीं करता है। मेरे प्रश्न:

1) क्या मुझे कुछ याद आ रही है? क्या कुछ त्रुटि हैंडलिंग (एआरएफसी के लिए) को लागू करना आवश्यक है?

2) क्या आपके पास कोड या विकल्प को गति देने के लिए कोई अन्य सुझाव है? क्या यह फॉर-लूप समांतर करने के लिए भुगतान करता है?

require(Rcpp) 
require(RcppArmadillo) 
require(microbenchmark) 

#Example 1 : standard R::pnorm 
src1 <- ' 
NumericVector ppnorm(const arma::vec& x,const arma::vec& mu,const  arma::vec& sigma, int lt, int lg) { 
int n = x.size(); 
arma::vec res(n); 
for (int i=0; i<n; i++) { 
    res(i) = R::pnorm(x(i),mu(i),sigma(i),lt,lg); 
} 
return wrap(res); 
} 
' 

#Example 2: approximation with ercf 
src2 <- ' 
NumericVector ppnorm(const arma::vec& x,const arma::vec& mu,const arma::vec& sigma, int lt, int lg) { 
int n = x.size(); 
arma::vec res(n); 
for (int i=0; i<n; i++) { 
res(i) = 0.5 * erfc(-(x(i) - mu(i))/sigma(i) * M_SQRT1_2); 
} 
if (lt==0 & lg==0) { 
    return wrap(1 - res); 
} 
if (lt==1 & lg==0) { 
    return wrap(res); 
} 
if (lt==0 & lg==1) { 
    return wrap(log(1 - res)); 
} 
if (lt==1 & lg==1) { 
    return wrap(log(res)); 
} 
} 
' 

#some random numbers 
xex = rnorm(100,5,4) 
muex = rnorm(100,3,1) 
siex = rnorm(100,0.8,0.3) 

#compile c++ functions 
func1 = cppFunction(depends = "RcppArmadillo",code=src1) #R::pnorm 
func2 = cppFunction(depends = "RcppArmadillo",code=src2) #ercf 

#run with exemplaric data 
res1 = func1(xex,muex,siex,1,0) 
res2 = func2(xex,muex,siex,1,0) 

# sum of squared errors 
sum((res1 - res2)^2,na.rm=T) 
# 6.474419e-32 ... very small 

#benchmarking 
microbenchmark(func1(xex,muex,siex,1,0),func2(xex,muex,siex,1,0),times=10000) 
#Unit: microseconds 
#expr min  lq  mean median  uq  max neval 
#func1(xex, muex, siex, 1, 0) 11.225 11.9725 13.72518 12.460 13.617 103.654 10000 
#func2(xex, muex, siex, 1, 0) 8.360 9.1410 10.62114 9.669 10.769 205.784 10000 
#my machine: Ubuntu 14.04 LTS, i7 2640M 2.8 Ghz x 4, 8GB memory, RRO 3.2.0 based on version R 3.2.0 

Answer 1

1) ठीक है, तुम सच में उपयोग करना चाहिए आर के pnorm() अपने 0-वें उदाहरण के रूप में। आप नहीं करते हैं, आप इसे आरसीपीपी इंटरफ़ेस का उपयोग करते हैं। आर का pnorm() पहले से ही अच्छी तरह से आर-आंतरिक रूप से वेक्टरकृत है (यानी सी स्तर पर) इसलिए आरसीपीपी की तुलना में तुलनात्मक या तेज भी हो सकता है। यह इसके अलावा लाभ एनए, NaN, Inf, आदि के मामलों को कवर करने के ..

2) यदि आप MLE बारे में बात कर रहे हैं, तो करता है, और आप गति और सटीकता के बारे में चिंतित हैं, आप लगभग निश्चित रूप से नहीं बल्कि साथ काम करना चाहिए लॉगरिदम, और शायद pnorm() के साथ नहीं बल्कि dnorm()?