In dieser Sitzung werde ich eine schnelle, aber genaue Annäherung an pnorm() demonstrieren.
Bevor wir beginnen, müssen wir uns klar machen: Was wollen wir erreichen, indem wir Approximation verwenden? Effizienz/Geschwindigkeit/Leistung, oder? Aber woher sollte diese Effizienz kommen?
Wie oben diskutiert, pnorm() Berechnung ist speichergebunden; Selbst wenn wir die Berechnung annähern, ist sie immer noch speichergebunden (daher berücksichtigen wir keine weitere Parallelisierung). Speichergebundene Probleme haben
number of floating point operations : memory access = O(1)
Mit anderen Worten, dieses Verhältnis eine Konstante C. Unser Ziel sollte daher sein, diese Konstante zu reduzieren, d. H. Wir wollen Gleitkommaoperationen reduzieren.
Die Anzahl der Gleitkommaoperationen wird oft als Anzahl der Gleitkommaadditionen und Multiplikationen angegeben. Andere Arten von Gleitkommaoperationen werden zu einem solchen Maß "umgewandelt". Vergleichen wir nun die Kosten mehrerer allgemeiner Gleitkommaoperationen.
u <- sample(1:10/10, 5e+7, replace = TRUE)
system.time(u + u)
# user system elapsed
# 0.468 0.168 0.639
system.time(u * u)
# user system elapsed
# 0.424 0.212 0.638
system.time(u/u)
# user system elapsed
# 0.504 0.204 0.710
system.time(u^1.1)
# user system elapsed
# 7.240 0.212 7.458
system.time(sqrt(u))
# user system elapsed
# 2.044 0.176 2.224
system.time(exp(u))
# user system elapsed
# 4.336 0.208 4.550
system.time(log(u))
# user system elapsed
# 2.748 0.172 2.925
system.time(round(u))
# user system elapsed
# 6.836 0.188 7.034
Beachten Sie, dass Addition und Multiplikation ist billig, Wurzel und Logarithmus ist teurer, während einige Operationen sehr teuer sind, einschließlich Energie, exponentielle und Rundung.
Jetzt zurück zu pnorm(), oder sogar dnorm(), usw., wo wir einen exponentiellen Term zu berechnen haben. Da:
system.time(pnorm(u))
# user system elapsed
# 11.016 0.160 11.193
system.time(dnorm(u))
# user system elapsed
# 8.844 0.164 9.022
wir sehen, dass die Mehrheit der Zeit zu berechnen pnorm() und dnorm() ist Attribut Berechnung exponentiell. pnorm() dauert länger als dnorm() weil es weiter ein Integral hat!
Jetzt ist unser Ziel ziemlich klar: Wir wollen die teure pnorm() Auswertung durch etwas wirklich billiges ersetzen, im Idealfall nur Addition/Multiplikation. Können wir??
Es gab viele Näherungsmethoden in der Geschichte. @Ben hat die logistische Approximation erwähnt. R-Funktion plogis() tut dies. Aber ein sorgfältiges Lesen auf ?plogis zeigt, dass es immer noch auf Exponentialfunktionen basiert.
Können wir nun anstelle der parametrischen Approximation eine nichtparametrische Approximation durchführen?Aber wir sollten hier keine Regression machen; stattdessen wollen wir eine Interpolationsfunktion von präzisen Daten mit feiner Auflösung verwenden, um pnorm() vorherzusagen. Nun, lineare Interpolation ist die beste Wahl, da sie super effizient ist (aufgrund von Sparsität: die lineare Prädiktormatrix ist tri-diagonal). In R macht approx dies. Ich verweise Leser, der damit nicht vertraut ist, auf ?approx, und ich werde einfach fortfahren.
Das OP sagt, er braucht nur Standard-Normalverteilung, also konzentrieren wir uns darauf. Betrachten Sie die folgende Näherungsfunktion (ich habe nicht approxfun verwenden, weil ich anpassbare h wollen):
approx.pnorm <- function(u, h = 0.2) {
x <- seq(from = -4, to = 4, by = h)
approx(x, pnorm(x), yleft = 0, yright = 1, xout = u)$y
}
Die genauen Daten auf einem Raster von Auflösung h zwischen [-4, 4] getroffen werden. Vorhersagen unter -4 ist 0, während Vorhersagen über 4 1 ist. Dies erfüllt die Anforderung einer CDF. Bei neuen Werten u nähern wir uns pnorm(u) durch lineare Interpolation basierend auf bekannten genauen Daten.
Offensichtlich steuert die Auflösung h Genauigkeit. Betrachten Sie die folgende Funktion RMSE und Anzeigenäherungskurve zu berechnen:
RMSEh <- function(h) {
x <- sort(rnorm(1000))
y <- pnorm(x)
y1 <- approx.pnorm(x, h)
plot(x, y, type = "l", lwd = 2); lines(x, y1, col = 2, lwd = 2)
mean((y - y1)^2)^0.5
}
par(mfrow = c(1, 3))
RMSEh(1) # 0.01570339
RMSEh(0.5) # 0.003968882
RMSEh(0.2) # 0.000639888
Eigentlich wenn h = 0.2, Annäherung schon ziemlich gut ist. Also werden wir im folgenden h = 0.2 verwenden.
Benchmarking
Dies sollte der spannendste Teil sein. In oben haben wir gesehen, dass die genaue Berechnung von pnorm(u) 11 Sekunden dauert. Jetzt
system.time(approx.pnorm(u, h = 0.2))
# user system elapsed
# 2.656 0.172 2.833
Wow, wir sind fast 4 mal schneller !!
Abhängig von Ihrem Anwendungsfall und Ihrer Genauigkeitsanforderung können Sie (1) eine adäquate, aber schnellere Annäherung an "pnorm" (2) "memoise" verwenden, d. H. Einige Ergebnisse in einer Nachschlagetabelle speichern. z.B. eine Old-School-Approximation zu "pnorm" ist "plogis" mit einer entsprechend angepassten Varianz ... –
z.B. [this] (http://www.jiem.org/index.php/jiem/article/viewFile/60/27) sieht nützlich aus ... –
Approximation ist eine gute Idee, denke ich, vor allem nach @ ZheyuanLi Antwort. – clemlaflemme