Schnellste Methode zur Berechnung der CDF der Normalverteilung über Vektoren - R :: pnorm vs erfc vs?

Question

Ich hoffe, meine umformulierte Frage entspricht jetzt den Kriterien von Stackoverflow. Bitte beachten Sie das Beispiel unten. Ich schreibe eine Log-Likelihood-Funktion, in der die Berechnung der cdf über Vektoren der zeitaufwendigste Teil ist. Beispiel 1 verwendet die R::pnorm, Beispiel 2 approximiert die normale cdf mit erfc. Wie Sie sehen können, sind die Ergebnisse ausreichend ähnlich, die ercf-Version ist ein bisschen schneller.

In der Praxis (innerhalb einer MLE) jedoch stellt sich heraus, dass die Ercf nicht so präzise ist, was den Algorithmus in Inf-Bereiche laufen lässt, wenn man die Bedingungen nicht genau festlegt. Meine Fragen:

1) Fehle ich etwas? Ist es notwendig, eine Fehlerbehandlung (für die erfc) zu implementieren?

2) Haben Sie weitere Vorschläge, um den Code oder Alternativen zu beschleunigen? Lohnt es sich, die For-Schleife parallelisieren zu sehen?

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) Nun, sollten Sie wirklich R verwenden pnorm() als 0-te Beispiel. Sie nicht, verwenden Sie die Rcpp-Schnittstelle dazu. R pnorm() ist bereits gut R-intern vektorisiert (d. H. Auf C-Niveau), kann also durchaus vergleichbar oder sogar schneller als Rcpp sein. Auch es hat haben den Vorteil, Fälle von NA, NaN, Inf, etc ..

2) Wenn Sie über MLE sprechen, und Sie sind besorgt über Geschwindigkeit und Genauigkeit, sollten Sie fast sicher lieber mit arbeiten die Logarithmen, und vielleicht nicht mit pnorm() sondern dnorm()?