机器算法验证 - Box-Muller 相对于逆 CDF 方法模拟正态分布的优势？ - 吾爱随笔录

Box-Muller 相对于逆 CDF 方法模拟正态分布的优势？

机器算法验证正态分布模拟均匀分布

2022-02-07 22:22:44

为了从一组均匀变量中模拟正态分布，有几种技术：

Box-Muller 算法，其中一个对两个独立的均匀变量进行采样 $(0,1)$ 并通过以下方式将它们转换为两个独立的标准正态分布：
$Z_{0} = \sqrt{- 2 ln U_{1}} cos (2 π U_{0}) Z_{1} = \sqrt{- 2 ln U_{1}} sin (2 π U_{0})$ $Z_0 = \sqrt{-2\text{ln}U_1}\text{cos}(2\pi U_0)\\ Z_1 = \sqrt{-2\text{ln}U_1}\text{sin}(2\pi U_0)$
CDF方法，可以等同于普通的cdf $(F(Z))$ 到统一变量：
$F (Z) = U$ $F(Z) = U$ 并得出 $Z = F^{- 1} (U)$ $Z = F^{-1}(U)$

我的问题是：哪个计算效率更高？我认为是后一种方法——但我读过的大多数论文都使用 Box-Muller——为什么？

附加信息：

正常 CDF 的倒数已知并由下式给出：

F^{- 1} (Z) = \sqrt{2} \erf^{- 1} (2 Z - 1), Z \in (0, 1) .

$F^{-1}(Z)\; =\; \sqrt2\;\operatorname{erf}^{-1}(2Z - 1), \quad Z\in(0,1).$

因此：

Z = F^{- 1} (U) = \sqrt{2} \erf^{- 1} (2 U - 1), U \in (0, 1) .

$Z = F^{-1}(U)\; =\; \sqrt2\;\operatorname{erf}^{-1}(2U - 1), \quad U\in(0,1).$

1个回答

从纯概率的角度来看，这两种方法都是正确的，因此是等效的。从算法的角度来看，比较必须同时考虑精度和计算成本。

Box-Muller 依赖于一个统一的生成器，成本与这个统一的生成器大致相同。正如我在评论中提到的，如果不是没有对数，你可以不用正弦或余弦调用就可以逃脱：

产生 $U_{1}, U_{2} \overset{iid}{\sim} U (- 1, 1)$ $U_1,U_2\stackrel{\text{iid}}{\sim}\mathcal{U}(-1,1)$ 直到 $S=U_1^2+U_2^2\le 1$
拿 $Z=\sqrt{-2\log(S)/S}$ 并定义 $X_{1} = Z U_{1}, X_{2} = Z U_{2}$ $X_1=ZU_1\,,\ X_2=Z U_2$

通用反演算法需要调用逆正态 cdf，例如qnorm(runif(N))在 R 中，这可能比上述更昂贵，更重要的是可能在精度方面失败，除非分位数函数编码良好。

按照whuber的评论，比较rnorm(N)和qnorm(runif(N))在执行时间上都具有逆 cdf 的优势：

> system.time(qnorm(runif(10^8)))
sutilisateur     système      écoulé
 10.137           0.120      10.251 
> system.time(rnorm(10^8))
utilisateur     système      écoulé
 13.417           0.060      13.472` `

使用更准确的 R 基准

        test replications elapsed relative user.self sys.self
3 box-muller          100   0.103    1.839     0.103        0
2    inverse          100   0.056    1.000     0.056        0
1    R rnorm          100   0.064    1.143     0.064        0

并在适合尾巴方面：在此处输入图像描述

在我的博客上对 Radford Neal 发表评论后，我想指出rnormR 中的默认值使用了反转方法，因此上述比较反映的是界面，而不是模拟方法本身！引用 RNG 上的 R 文档：

‘normal.kind’ can be ‘"Kinderman-Ramage"’, ‘"Buggy
 Kinderman-Ramage"’ (not for ‘set.seed’), ‘"Ahrens-Dieter"’,
 ‘"Box-Muller"’, ‘"Inversion"’ (the default), or ‘"user-supplied"’.
 (For inversion, see the reference in ‘qnorm’.)  The
 Kinderman-Ramage generator used in versions prior to 1.7.1 (now
 called ‘"Buggy"’) had several approximation errors and should only
 be used for reproduction of old results.  The ‘"Box-Muller"’
 generator is stateful as pairs of normals are generated and
 returned sequentially.  The state is reset whenever it is selected
 (even if it is the current normal generator) and when ‘kind’ is
 changed.

其它你可能感兴趣的问题

上一篇在随机森林中包含交互项下一篇相关与因果混淆的真实例子