在本次会议中，我将演示对pnorm() 的快速而准确的近似。
在开始之前，我们需要明确：我们想通过使用近似来实现什么？效率/速度/性能，对吧？但这样的效率从何而来？
如上所述，pnorm() 计算受内存限制；即使我们进行近似计算，它仍然受内存限制（因此我们不考虑进一步并行化）。内存绑定问题有
number of floating point operations : memory access = O(1)
换句话说，这个比率是一个常数C。所以我们的目标应该是减少这个常数，即我们要减少浮点运算。
浮点运算的次数通常报告为浮点加法和乘法的次数。其他类型的浮点运算被“转换”为这种度量。现在，我们来比较几种常见的浮点运算的成本。
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
请注意，加法和乘法很便宜，根和对数更昂贵，而一些运算非常昂贵，包括幂、指数和舍入。
现在让我们回到pnorm()，甚至dnorm() 等，我们需要计算一个指数项。鉴于：
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
我们看到计算pnorm() 和dnorm() 的大部分时间都归因于计算指数。 pnorm() 比 dnorm() 需要更长的时间，因为它还有一个积分！
现在，我们的目标相当明确：我们想用非常便宜的东西替换昂贵的pnorm() 评估，理想情况下只涉及加法/乘法。我们可以吗？？
历史上有很多近似方法。 @Ben 提到了逻辑近似。 R 函数plogis() 执行此操作。但仔细阅读?plogis 会发现它仍然是基于指数的。
现在，我们可以使用非参数逼近来代替那些参数逼近吗？但我们不应该在这里做回归；相反，我们想使用一些高分辨率精确数据的插值函数来预测pnorm()。好吧，线性插值是最好的选择，因为它非常有效（由于稀疏性：线性预测矩阵是三对角矩阵）。在 R 中，approx 执行此操作。我将不熟悉的读者推荐给?approx，我将继续进行。
OP 说他只需要标准正态分布，所以我们专注于此。考虑以下近似函数（我没有使用approxfun，因为我想要自定义h）：
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
}
准确的数据是在[-4, 4] 之间的分辨率为h 的网格上获取的。 -4 以下的预测为 0，4 以上的预测为 1。这满足了 CDF 的要求。给定新值u，我们根据已知准确数据通过线性插值近似pnorm(u)。
显然，分辨率h 控制着准确性。考虑以下函数来计算 RMSE 并显示近似曲线：
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
其实h = 0.2的时候，逼近已经相当不错了。所以下面我们将使用h = 0.2。
基准测试
这应该是最激动人心的部分。在上面我们已经看到pnorm(u) 的准确计算需要 11 秒。现在
system.time(approx.pnorm(u, h = 0.2))
# user system elapsed
# 2.656 0.172 2.833
哇，我们快了将近 4 倍！！