我的问题受到R的内置指数随机数生成器 function 的启发rexp()。在尝试生成指数分布的随机数时，许多教科书推荐使用此 Wikipedia 页面中概述的逆变换方法。我知道还有其他方法可以完成这项任务。特别是，R的源代码使用了Ahrens & Dieter (1972) 的论文中概述的算法。
我已经说服自己，Ahrens-Dieter (AD) 方法是正确的。不过，与逆变换 (IT) 方法相比，我看不到使用他们的方法的好处。AD 不仅比 IT 实施起来更复杂。似乎也没有速度优势。这是我的R代码，用于对这两种方法进行基准测试，然后是结果。
invTrans <- function(n)
-log(runif(n))
print("For the inverse transform:")
print(system.time(invTrans(1e8)))
print("For the Ahrens-Dieter algorithm:")
print(system.time(rexp(1e8)))
结果：
[1] "For the inverse transform:"
user system elapsed
4.227 0.266 4.597
[1] "For the Ahrens-Dieter algorithm:"
user system elapsed
4.919 0.265 5.213
比较两种方法的代码，AD 至少抽取两个均匀随机数（用C函数unif_rand()）得到一个指数随机数。IT 只需要一个统一的随机数。据推测，R核心团队决定不实施 IT，因为它认为取对数可能比生成更统一的随机数要慢。我知道取对数的速度可能取决于机器，但至少对我来说恰恰相反。也许 IT 的数值精度存在与 0 处的对数奇异性有关的问题？但是，R
源代码sexp.c揭示了 AD 的实现也失去了一些数值精度，因为 C 代码的以下部分从统一随机数u中删除了前导位。
double u = unif_rand();
while(u <= 0. || u >= 1.) u = unif_rand();
for (;;) {
u += u;
if (u > 1.)
break;
a += q[0];
}
u -= 1.;
u稍后在sexp.c的其余部分中作为统一随机数回收。到目前为止，似乎
- IT更容易编码，
- IT 速度更快，并且
- IT 和 AD 都可能失去数值准确性。
如果有人能解释为什么 R 仍然将 AD 作为rexp().
