模行为背后的数学答案

【问题标题】：mathematics behind modulo behavor模行为背后的数学
【发布时间】：2013-01-30 22:10:13
【问题描述】：

序言

这个问题与 (P)RNG 和 rand() 的行为无关。这是关于使用对模均匀分布的两个值的幂。

简介

我知道不应该使用模 % 将一个值从一个范围转换为另一个，例如从 rand() 函数获取一个介于 0 和 5 之间的值：会有偏差。这里解释了https://bitbucket.org/haypo/hasard/src/ebf5870a1a54/doc/common_errors.rst?at=default 和这个答案Why do people say there is modulo bias when using a random number generator?

但是今天在调查了一些看起来有问题的代码之后，我做了一个工具来演示模数的行为：https://gitorious.org/modulo-test/modulo-test/trees/master，发现还不够清楚。

骰子只有 3 位

我检查了 0..5 范围内的 6 个值。对这些值进行编码只需要 3 位。

$ ./modulo-test 10000 6 3
interations = 10000, range = 6, bits = 3 (0x00000007)
  [0..7] => [0..5]

theorical occurences    1666.67 probability 0.16666667

   [   0] occurences    2446    probability 0.24460000 ( +46.76%)
   [   1] occurences    2535    probability 0.25350000 ( +52.10%)
   [   2] occurences    1275    probability 0.12750000 ( -23.50%)
   [   3] occurences    1297    probability 0.12970000 ( -22.18%)
   [   4] occurences    1216    probability 0.12160000 ( -27.04%)
   [   5] occurences    1231    probability 0.12310000 ( -26.14%)

  minimum occurences    1216.00 probability 0.12160000 ( -27.04%)
  maximum occurences    2535.00 probability 0.25350000 ( +52.10%)
     mean occurences    1666.67 probability 0.16666667 (  +0.00%)
   stddev occurences     639.43 probability 0.06394256 (  38.37%)

使用 3 位输入，结果确实很糟糕，但表现符合预期。见答案https://stackoverflow.com/a/14614899/611560

增加输入位数

让我感到困惑的是，增加输入位数会导致结果不同。您不应该忘记增加迭代次数，例如样本数，否则结果可能是错误的（请参阅错误统计数据）。

让我们尝试 4 位：

$ ./modulo-test 20000 6 4
interations = 20000, range = 6, bits = 4 (0x0000000f)
  [0..15] => [0..5]

theorical occurences    3333.33 probability 0.16666667

   [   0] occurences    3728    probability 0.18640000 ( +11.84%)
   [   1] occurences    3763    probability 0.18815000 ( +12.89%)
   [   2] occurences    3675    probability 0.18375000 ( +10.25%)
   [   3] occurences    3721    probability 0.18605000 ( +11.63%)
   [   4] occurences    2573    probability 0.12865000 ( -22.81%)
   [   5] occurences    2540    probability 0.12700000 ( -23.80%)

  minimum occurences    2540.00 probability 0.12700000 ( -23.80%)
  maximum occurences    3763.00 probability 0.18815000 ( +12.89%)
     mean occurences    3333.33 probability 0.16666667 (  +0.00%)
   stddev occurences     602.48 probability 0.03012376 (  18.07%)

让我们尝试 5 位：

$ ./modulo-test 40000 6 5
interations = 40000, range = 6, bits = 5 (0x0000001f)
  [0..31] => [0..5]

theorical occurences    6666.67 probability 0.16666667

   [   0] occurences    7462    probability 0.18655000 ( +11.93%)
   [   1] occurences    7444    probability 0.18610000 ( +11.66%)
   [   2] occurences    6318    probability 0.15795000 (  -5.23%)
   [   3] occurences    6265    probability 0.15662500 (  -6.03%)
   [   4] occurences    6334    probability 0.15835000 (  -4.99%)
   [   5] occurences    6177    probability 0.15442500 (  -7.34%)

  minimum occurences    6177.00 probability 0.15442500 (  -7.34%)
  maximum occurences    7462.00 probability 0.18655000 ( +11.93%)
     mean occurences    6666.67 probability 0.16666667 (  +0.00%)
   stddev occurences     611.58 probability 0.01528949 (   9.17%)

让我们试试 6 位：

$ ./modulo-test 80000 6 6
interations = 80000, range = 6, bits = 6 (0x0000003f)
  [0..63] => [0..5]

theorical occurences   13333.33 probability 0.16666667

   [   0] occurences   13741    probability 0.17176250 (  +3.06%)
   [   1] occurences   13610    probability 0.17012500 (  +2.08%)
   [   2] occurences   13890    probability 0.17362500 (  +4.18%)
   [   3] occurences   13702    probability 0.17127500 (  +2.77%)
   [   4] occurences   12492    probability 0.15615000 (  -6.31%)
   [   5] occurences   12565    probability 0.15706250 (  -5.76%)

  minimum occurences   12492.00 probability 0.15615000 (  -6.31%)
  maximum occurences   13890.00 probability 0.17362500 (  +4.18%)
     mean occurences   13333.33 probability 0.16666667 (  +0.00%)
   stddev occurences     630.35 probability 0.00787938 (   4.73%)

问题

请解释一下为什么在更改输入位（并相应地增加样本数）时结果会有所不同？这些背后的数学推理是什么？

错误的统计数据

在之前版本的问题中，我展示了一个 32 位输入和只有 1000000 次迭代的测试，例如 10^6 个样本，并说我很惊讶得到正确的结果。太错误了，我很惭愧：必须有 N 倍以上的样本才能有信心获得生成器的所有 2^32 值。这里 10^6 与 2^32 相比要小得多。 奖励能够用数学/统计语言解释这一点的人。。

这里是错误的结果：

$ ./modulo-test 1000000 6 32
interations = 1000000, range = 6, bits = 32 (0xffffffff)
  [0..4294967295] => [0..5]

theorical occurences  166666.67 probability 0.16666667

   [   0] occurences  166881    probability 0.16688100 (  +0.13%)
   [   1] occurences  166881    probability 0.16688100 (  +0.13%)
   [   2] occurences  166487    probability 0.16648700 (  -0.11%)
   [   3] occurences  166484    probability 0.16648400 (  -0.11%)
   [   4] occurences  166750    probability 0.16675000 (  +0.05%)
   [   5] occurences  166517    probability 0.16651700 (  -0.09%)

  minimum occurences  166484.00 probability 0.16648400 (  -0.11%)
  maximum occurences  166881.00 probability 0.16688100 (  +0.13%)
     mean occurences  166666.67 probability 0.16666667 (  +0.00%)
   stddev occurences     193.32 probability 0.00019332 (   0.12%)

Zed Shaw "Programmers Need To Learn Statistics Or I Will Kill Them All"的优秀文章我还是要反复阅读。

【问题讨论】：

这一切都没有意义吗？我的意思是，你不希望他们与众不同吗？ 2/8 比 4/4294967296 大得多。
rand() % N 不建议使用，因为rand 的某些（较旧？）实现在低位内的顺序结果运行中显示模式。即使rand() 在您的实现中确实是随机的，但这并不意味着每个人都可以安全地依赖它。
@aschepler，我认为他是在询问如果生成的随机数范围不能被您试图获得的范围整除，那么对较低结果的固有偏见。这本身并不受随机数生成器实现的随机性（或缺乏随机性）的影响。
给你：stackoverflow.com/questions/2509679/…
我在en.wikipedia.org/wiki/Fisher%E2%80%93Yates_shuffle#Modulo_bias找到了一些解释

标签： math random statistics probability modulo

【解决方案1】：

本质上，你在做：

(rand() & 7) % 6

假设rand() 均匀分布在[0; RAND_MAX] 上，而RAND_MAX+1 是2 的幂。很明显，rand() & 7 可以评估为 0、1、...、7，并且结果是等概率的。

现在让我们看看将结果取模 6 时会发生什么。

0 和 6 映射到 0；
1 和 7 映射到 1；
2 映射到 2；
3 映射到 3；
4 映射到 4；
5 映射到 5。

这就解释了为什么你得到的零和一是其他数字的两倍。

第二种情况也发生了同样的事情。但是，“额外”数字的值要小得多，使得它们的贡献与噪声无法区分。

总而言之，如果你有一个整数均匀分布在 [0; M-1]，然后你取模 N，结果将偏向于零，除非 M 可以被 N 整除。

【讨论】：

+1，同样的事情也发生在 32 位的情况下。只是偏差不太明显，因为 40 亿中只有 4 个“额外”值，而不是 8 中的 2 个。
感谢您的解释。我知道基本行为，并且觉得没有必要在问题中描述它。关于 32 个输入位的示例，由于样本数减少，这是错误的。所以这基本上是我的错误。我更改了问题以修复样本计数并提出不同但相关的问题。

【解决方案2】：

rand()（或其他一些 PRNG）在区间 [0 .. RAND_MAX] 中生成值。您想使用余数运算符将这些映射到区间 [0 .. N-1]。

写

(RAND_MAX+1) = q*N + r

0 <= r < N.

那么对于区间[0 .. N-1]中的每个值都有

q+1 的值如果小于 r，则映射到该值的 rand()
如果值为>= r，则映射到rand() 的q 值。

现在，如果q 很小，则q 和q+1 之间的相对差异很大，但如果q 很大 - 例如2^32 / 6 - 差异就不容易测量。

【讨论】：