当好的伪随机数变坏时：Multivariate Normal Sampling的坑

如果你不小心，多元正态分布抽样可能变得非常不可复现。有时比其他时候更严重。这里面有奇怪的恐怖，病态矩阵和浮点数的梦魇。教沙子做线性代数是个错误。

关键词: R, Reproducibility, Statistics

作者: Danielle Navarro

发布时间: 2025年5月18日

目录

在真实的计算机上，计算矩阵的特征分解会受到误差的影响：权威的分析是 Wilkinson (1965)。 您只能希望找到一个与 x 足够接近的问题的解决方案。 因此，即使一个真实的非对称 x 可能具有带有重复实特征值的代数解，但计算出的解可能是一个具有复共轭对特征值的相似矩阵的解。 – help("eigen")

几周前，在我还没有亲身经历感染 covid-19 带来的灵魂毁灭般的痛苦之前，1 一些同事2 找到我，谈论他们在一些 R 代码中遇到的可重复性问题。 他们一直在运行依赖于从多元正态分布生成样本的模拟，尽管做了谨慎的事情并使用 set.seed() 来控制随机数生成器 (RNG) 的状态，但结果在计算上不可重现。 相同的代码在不同的机器上执行会产生不同的随机数。 这些数字并没有像我们都疲惫地了解到当您试图强迫计算机进行数学运算时那样“只是有点不同”。 它们是痛苦地，残酷地，灾难性地，不可重现地不同。 不知何处，不知何故，出了问题。3 4 这不应该发生…… 不知何故，“Palpatine” 又回来了。 Palpatine 一直都是浮点错误。

……公平地说，在大多数情况下，它_没有_发生。 大多数计算都是可重复的，即使是随机数生成，只要您小心使用种子，也会按照它应该做的那样工作。 或者，更简单地说，用于控制 R 中计算结果的 set.seed() 方法工作得很好。 冒着让人厌烦的风险，这里是用 R 在固定种子后获得的英文字母的“随机”排列：

[](https://blog.djnavarro.net/posts/2025-05-18_multivariate-normal-sampling-floating-point/<https:/blog.djnavarro.net/posts/2025-05-18_multivariate-normal-sampling-floating-point/#cb1-1>)set.seed(1L)
[](https://blog.djnavarro.net/posts/2025-05-18_multivariate-normal-sampling-floating-point/<https:/blog.djnavarro.net/posts/2025-05-18_multivariate-normal-sampling-floating-point/#cb1-2>)sample(letters)

 [1] "y" "d" "g" "a" "b" "k" "n" "r" "w" "j" "f" "t" "q" "x" "i" "e" "u" "l"
[19] "s" "p" "o" "m" "v" "z" "c" "h"

这是这些字母的_完全相同的排列_，第二次产生是因为作者感到无聊，并且她希望您理解这是一个 quarto 文档，因此您在此处看到的一切都是在实际的 R 会话中完成的，并且她没有虚构：

[](https://blog.djnavarro.net/posts/2025-05-18_multivariate-normal-sampling-floating-point/<https:/blog.djnavarro.net/posts/2025-05-18_multivariate-normal-sampling-floating-point/#cb3-1>)set.seed(1L)
[](https://blog.djnavarro.net/posts/2025-05-18_multivariate-normal-sampling-floating-point/<https:/blog.djnavarro.net/posts/2025-05-18_multivariate-normal-sampling-floating-point/#cb3-2>)sample(letters)

 [1] "y" "d" "g" "a" "b" "k" "n" "r" "w" "j" "f" "t" "q" "x" "i" "e" "u" "l"
[19] "s" "p" "o" "m" "v" "z" "c" "h"

因为原始代码和复制代码都在调用 sample() 将字母表中的字母洗牌成随机顺序之前使用了 set.seed(1)，所以在两种情况下我们都得到了_相同的_随机排列。 尽管我在同一台机器上执行此代码两次，但没有理由期望如果我在运行 R 4.4.3 的我的 ubuntu 笔记本电脑或运行 R 4.3.1 的 window 机器上运行原始代码会有所不同。 无论哪种方式，它都“应该”是相同的结果。 我的意思是……这就是 set.seed() 的全部目的。 它“应该”使您的代码可重现，即使您使用该代码执行“随机”操作。 对吧？ 好吧。 它可能会。 如果你幸运的话。

当生成多元随机正态样本时，这并不是我的同事用他们的代码所经历的。 在他们的情况下，即使代码使用 set.seed() 来（表面上）防止这种情况发生，不同的机器也会产生不同的随机样本。 相当合理的是，他们的第一个想法是“这他妈的是什么???”。

遇到问题后，他们的第二个想法是这是他们的错：_他们_一定在他们的代码中做错了什么，导致模拟中断。 事实并非如此。 他们的代码没问题。 他们的第三个想法是 R 本身有问题，或者更准确地说是 MASS::mvrnorm() 有问题，因为这是造成所有困难的函数。 但是他们看不到该函数有任何问题，所以他们请我调查一下。

当我如此天真地在我的办公桌上扔了一颗核武器时，我的第一个想法是 MASS 不可能是问题所在。 我的意思是，不是真的。 这是一个经过非常仔细测试的包，它被非常广泛地使用，它很好。 对吧？ 虽然……实际上，现在我想起来了，MASS::mvrnorm() 是一个非常旧的函数，并且它对较低级别发生的事情做出了一些假设，这些假设可能并不正确。 所以……也许这就是问题所在？？？ 我真的不知道…… 听着，这很复杂，好吗？

所以让我回头再跟你说一遍，因为当你追溯到足够远的时候，我的同事们遇到的问题的根本原因实际上并不是 MASS 包。 真正的问题原来是 浮点运算 的不便事实不像实数运算那样运行，这反过来又产生了计算机通常不做人们期望的事情的 tiny 极小的副作用。

唉。

我应该知道的。 事后看来，很明显会是浮点问题5，我不应该花几个小时/几天/几周的时间来调查。 正如我当时在 Mastodon 上所说的那样：

今后，每当有人问我“为什么[某事]即使我设置了 RNG 种子也无法重现？” 我甚至懒得去研究具体情况：我只会回复“浮点运算”，然后等到我最终被证明是正确的。

人们会对我的预知能力和我显然百科全书般的计算知识感到惊讶。 他们会赞扬我的天才，我的信徒会将我抬上轿子，我将像女王和女神一样生活在凡人之中……

矩阵是奇怪的物体。 对它们执行的在理论上讲得通的变换，当它们在物理世界中被执行时，有时会产生奇怪的效果。 （显然，这是我创作的艺术作品）

……好吧，也许不是。

难题

所以，好吧。 我的同事遇到的问题是“跨不同机器”的问题，因此我无法使用像渲染此帖子的单台机器轻松复制这个精确问题。 尽管如此，我可以使用一个“技巧”来近似这个问题，如果您以前从未遇到过它，这个技巧似乎非常聪明，但实际上非常平淡。6 考虑这两个协方差矩阵，cov1 和 cov2。 这是生成它们的代码：

[](https://blog.djnavarro.net/posts/2025-05-18_multivariate-normal-sampling-floating-point/<https:/blog.djnavarro.net/posts/2025-05-18_multivariate-normal-sampling-floating-point/#cb5-1>)cov1 <- matrix(
[](https://blog.djnavarro.net/posts/2025-05-18_multivariate-normal-sampling-floating-point/<https:/blog.djnavarro.net/posts/2025-05-18_multivariate-normal-sampling-floating-point/#cb5-2>) data = c(
[](https://blog.djnavarro.net/posts/2025-05-18_multivariate-normal-sampling-floating-point/<https:/blog.djnavarro.net/posts/2025-05-18_multivariate-normal-sampling-floating-point/#cb5-3>)  4.58, -1.07, 2.53, 0.14, 
[](https://blog.djnavarro.net/posts/2025-05-18_multivariate-normal-sampling-floating-point/<https:/blog.djnavarro.net/posts/2025-05-18_multivariate-normal-sampling-floating-point/#cb5-4>)  -1.07, 5.83, 1.15, -1.45, 
[](https://blog.djnavarro.net/posts/2025-05-18_multivariate-normal-sampling-floating-point/<https:/blog.djnavarro.net/posts/2025-05-18_multivariate-normal-sampling-floating-point/#cb5-5>)  2.53, 1.15, 2.26, -0.79, 
[](https://blog.djnavarro.net/posts/2025-05-18_multivariate-normal-sampling-floating-point/<https:/blog.djnavarro.net/posts/2025-05-18_multivariate-normal-sampling-floating-point/#cb5-6>)  0.14, -1.45, -0.79, 4.93
[](https://blog.djnavarro.net/posts/2025-05-18_multivariate-normal-sampling-floating-point/<https:/blog.djnavarro.net/posts/2025-05-18_multivariate-normal-sampling-floating-point/#cb5-7>) ), 
[](https://blog.djnavarro.net/posts/2025-05-18_multivariate-normal-sampling-floating-point/<https:/blog.djnavarro.net/posts/2025-05-18_multivariate-normal-sampling-floating-point/#cb5-8>) nrow = 4L, 
[](https://blog.djnavarro.net/posts/2025-05-18_multivariate-normal-sampling-floating-point/<https:/blog.djnavarro.net/posts/2025-05-18_multivariate-normal-sampling-floating-point/#cb5-9>) ncol = 4L
[](https://blog.djnavarro.net/posts/2025-05-18_multivariate-normal-sampling-floating-point/<https:/blog.djnavarro.net/posts/2025-05-18_multivariate-normal-sampling-floating-point/#cb5-10>))
[](https://blog.djnavarro.net/posts/2025-05-18_multivariate-normal-sampling-floating-point/<https:/blog.djnavarro.net/posts/2025-05-18_multivariate-normal-sampling-floating-point/#cb5-11>)
[](https://blog.djnavarro.net/posts/2025-05-18_multivariate-normal-sampling-floating-point/<https:/blog.djnavarro.net/posts/2025-05-18_multivariate-normal-sampling-floating-point/#cb5-12>)set.seed(1L)
[](https://blog.djnavarro.net/posts/2025-05-18_multivariate-normal-sampling-floating-point/<https:/blog.djnavarro.net/posts/2025-05-18_multivariate-normal-sampling-floating-point/#cb5-13>)tol <- 10e-13
[](https://blog.djnavarro.net/posts/2025-05-18_multivariate-normal-sampling-floating-point/<https:/blog.djnavarro.net/posts/2025-05-18_multivariate-normal-sampling-floating-point/#cb5-14>)eps <- matrix(rnorm(16L) * tol, 4L, 4L)
[](https://blog.djnavarro.net/posts/2025-05-18_multivariate-normal-sampling-floating-point/<https:/blog.djnavarro.net/posts/2025-05-18_multivariate-normal-sampling-floating-point/#cb5-15>)eps <- eps + t(eps)
[](https://blog.djnavarro.net/posts/2025-05-18_multivariate-normal-sampling-floating-point/<https:/blog.djnavarro.net/posts/2025-05-18_multivariate-normal-sampling-floating-point/#cb5-16>)
[](https://blog.djnavarro.net/posts/2025-05-18_multivariate-normal-sampling-floating-point/<https:/blog.djnavarro.net/posts/2025-05-18_multivariate-normal-sampling-floating-point/#cb5-17>)cov2 <- cov1 + eps

打印出来时，它们看起来相同：

[](https://blog.djnavarro.net/posts/2025-05-18_multivariate-normal-sampling-floating-point/<https:/blog.djnavarro.net/posts/2025-05-18_multivariate-normal-sampling-floating-point/#cb6-1>)cov1

   [,1] [,2] [,3] [,4]
[1,] 4.58 -1.07 2.53 0.14
[2,] -1.07 5.83 1.15 -1.45
[3,] 2.53 1.15 2.26 -0.79
[4,] 0.14 -1.45 -0.79 4.93

[](https://blog.djnavarro.net/posts/2025-05-18_multivariate-normal-sampling-floating-point/<https:/blog.djnavarro.net/posts/2025-05-18_multivariate-normal-sampling-floating-point/#cb8-1>)cov2

   [,1] [,2] [,3] [,4]
[1,] 4.58 -1.07 2.53 0.14
[2,] -1.07 5.83 1.15 -1.45
[3,] 2.53 1.15 2.26 -0.79
[4,] 0.14 -1.45 -0.79 4.93

但是，当然，既然您已经看到了代码，您会毫不惊讶地发现 cov2 实际上是 cov1 的一个非常轻微的扰动版本：

[](https://blog.djnavarro.net/posts/2025-05-18_multivariate-normal-sampling-floating-point/<https:/blog.djnavarro.net/posts/2025-05-18_multivariate-normal-sampling-floating-point/#cb10-1>)cov1 - cov2

       [,1]     [,2]     [,3]     [,4]
[1,] 1.253220e-12 -5.131451e-13 2.597922e-13 -9.740264e-13
[2,] -5.131451e-13 1.641354e-12 -1.820766e-13 1.476375e-12
[3,] 2.597922e-13 -1.820766e-13 -3.023359e-12 -1.514788e-12
[4,] -9.740264e-13 1.476375e-12 -1.514788e-12 8.970602e-14

当发生浮点截断误差时，我们会遇到像这样的微小差异。7 为了使用经典的例子，当我们尝试计算一个简单的和，比如 0.1 + 0.2 - 0.3，结果应该是零（duh），但是如果我们在计算机上这样做，通常不是。 我们_实际_得到的答案是非常轻微的错误，因为 0.1 在浮点表示中的二进制表示8 是无限长的，不能用双精度浮点数完全表示，所以……

[](https://blog.djnavarro.net/posts/2025-05-18_multivariate-normal-sampling-floating-point/<https:/blog.djnavarro.net/posts/2025-05-18_multivariate-normal-sampling-floating-point/#cb12-1>)0.1 + 0.2 - 0.3

[1] 5.551115e-17

……答案并不_完全_正确。

更糟糕的是，应用浮点截断误差的计算的_精确_结果不一定是跨系统不变的。 操作系统差异、编译器设置以及许多其他因素都会影响结果。 这些细节对于这篇文章并不重要，坦率地说，我自己也不完全理解它们。 就我所知，笔记本电脑外壳的颜色，或者程序员的猫的名字可能有关。 一旦您的计算开始遇到浮点精度的限制，就会出现各种怪异现象。

那么，仅为了辩论起见，让我们假设在一些花哨的模拟过程中，您和我计算了不同机器上的协方差矩阵。 它应该是一个相同的协方差矩阵，但由于浮点数的怪异，您的机器计算出 cov1，而我的机器计算出 cov2。 这些差异非常小，但它们足够大，以至于——不知何故，违反了上帝和人类的所有法则——发生了这种情况：

[](https://blog.djnavarro.net/posts/2025-05-18_multivariate-normal-sampling-floating-point/<https:/blog.djnavarro.net/posts/2025-05-18_multivariate-normal-sampling-floating-point/#cb14-1>)set.seed(1L)
[](https://blog.djnavarro.net/posts/2025-05-18_multivariate-normal-sampling-floating-point/<https:/blog.djnavarro.net/posts/2025-05-18_multivariate-normal-sampling-floating-point/#cb14-2>)mvr1 <- MASS::mvrnorm(n = 1L, mu = rep(0L, 4L), Sigma = cov1)
[](https://blog.djnavarro.net/posts/2025-05-18_multivariate-normal-sampling-floating-point/<https:/blog.djnavarro.net/posts/2025-05-18_multivariate-normal-sampling-floating-point/#cb14-3>)mvr1

[1] -0.4391833 0.2560893 0.8542052 -2.2883238

[](https://blog.djnavarro.net/posts/2025-05-18_multivariate-normal-sampling-floating-point/<https:/blog.djnavarro.net/posts/2025-05-18_multivariate-normal-sampling-floating-point/#cb16-1>)set.seed(1L)
[](https://blog.djnavarro.net/posts/2025-05-18_multivariate-normal-sampling-floating-point/<https:/blog.djnavarro.net/posts/2025-05-18_multivariate-normal-sampling-floating-point/#cb16-2>)mvr2 <- MASS::mvrnorm(n = 1L, mu = rep(0L, 4L), Sigma = cov2)
[](https://blog.djnavarro.net/posts/2025-05-18_multivariate-normal-sampling-floating-point/<https:/blog.djnavarro.net/posts/2025-05-18_multivariate-normal-sampling-floating-point/#cb16-3>)mvr2

[1] 0.2878007 0.5942556 -0.2596185 -2.3919000

是的……呃……这些甚至根本不是一回事。

在文章的这一点上，您可能会根据您的背景做出以下两种反应之一。 如果您曾经有过阅读数值线性代数教科书并以某种方式幸存下来的创伤经历，那么您会疲惫地叹息并说“是的 Danielle，这就是您试图用沙子做数学运算时会发生的事情”。 但是，如果您生活在日常应用科学中阳光明媚的土地上，并且您的神灵没有被 IEEE-754 残酷谋杀，您可能会想一些更像是 “这他妈的是什么疯狂 Danielle?????????” 的话。

所以。 是的。 这是关于计算机的令人尴尬的事情之一。 当我们尝试生成具有多元正态分布的随机向量时，cov1 和 cov2 之间非常小的差异（有时，取决于用于进行采样的精确方法）会导致生成的数字产生很大的差异，即使 RNG 种子相同。 更奇怪的是，即使相关的随机数_都_具有正确的分布属性，这种情况也可能（并且确实）发生。 也就是说：结果是正确的，但它们不可重现。

这有点令人困惑。 依赖于随机数生成的其他类型的计算不会受到相同的影响。 直观上，我们期望_输入_参数中的微小差异会导致_输出_值中的微小差异。 事实上，如果我们从具有略微不同的标准差的单变量正态分布生成随机样本，就会发生这种情况：

[](https://blog.djnavarro.net/posts/2025-05-18_multivariate-normal-sampling-floating-point/<https:/blog.djnavarro.net/posts/2025-05-18_multivariate-normal-sampling-floating-point/#cb18-1>)s1 <- sqrt(cov1[1L, 1L])
[](https://blog.djnavarro.net/posts/2025-05-18_multivariate-normal-sampling-floating-point/<https:/blog.djnavarro.net/posts/2025-05-18_multivariate-normal-sampling-floating-point/#cb18-2>)s2 <- sqrt(cov2[1L, 1L])
[](https://blog.djnavarro.net/posts/2025-05-18_multivariate-normal-sampling-floating-point/<https:/blog.djnavarro.net/posts/2025-05-18_multivariate-normal-sampling-floating-point/#cb18-3>)
[](https://blog.djnavarro.net/posts/2025-05-18_multivariate-normal-sampling-floating-point/<https:/blog.djnavarro.net/posts/2025-05-18_multivariate-normal-sampling-floating-point/#cb18-4>)set.seed(1L)
[](https://blog.djnavarro.net/posts/2025-05-18_multivariate-normal-sampling-floating-point/<https:/blog.djnavarro.net/posts/2025-05-18_multivariate-normal-sampling-floating-point/#cb18-5>)r1 <- rnorm(n = 1L, mean = 0L, sd = s1)
[](https://blog.djnavarro.net/posts/2025-05-18_multivariate-normal-sampling-floating-point/<https:/blog.djnavarro.net/posts/2025-05-18_multivariate-normal-sampling-floating-point/#cb18-6>)r1

[1] -1.34067

[](https://blog.djnavarro.net/posts/2025-05-18_