Vektörler üzerinden Normal dağılımın cdf hesaplamak için hızlı yolu - R :: pnorm vs erfc vs?

Question

İnşallah sorularım şimdi Stackoverflow kriterlerine uyar. Lütfen aşağıdaki örneği dikkate alın. CDf'yi vektörler üzerinden hesaplamanın en çok zaman alan kısmı olduğu bir Log-Likelihood fonksiyonu yazdım. Örnek 1, R::pnorm'u kullanır, Örnek 2, erfc ile normal cdf'ye yaklaşır. Gördüğünüz gibi sonuçlar yeterince benzer, ercf versiyonu biraz daha hızlı. Pratikte (bir MLE içinde), ancak, ercf'nin kesin olmadığı, algoritmanın, kısıtlamaları doğru olarak ayarlamadığı sürece, inf alanlar içine girmesine izin verdiği anlaşılmaktadır. Sorularım:

1) Bir şey mi eksik? Bazı hata işleme (erfc için) uygulamak gerekli mi?

2) Kodu veya alternatifleri hızlandırmak için başka önerileriniz var mı? For döngüsü paralel hale getirmek için ödeme yapar mı?

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) Eh, gerçekten kullanmalıdır R en pnorm() senin 0-inci örnek olarak. Yapmıyorsunuz, buna Rcpp arayüzünü kullanıyorsunuz. R'nin pnorm() hali hazırda R-dahili olarak (yani C seviyesinde) güzel bir şekilde vektörlenmiş olup, bu yüzden Rcpp'den daha iyi veya daha hızlı olabilir. Ayrıca ,'un NA, NaN, Inf, vb. Durumlarını karşılama avantajına sahiptir.

2) MLE hakkında konuşuyorsanız ve hız ve doğruluktan endişe duyuyorsanız, hemen hemen kesinlikle çalışmak zorundasınız. logaritma ve belki de pnorm() değil, daha doğrusu dnorm()?