概率模拟误差不会收敛

Question

在一次采访中,我最初使用笔/纸解决了以下问题,然后通过程序验证结果.

问题如下:

有三个人A,B和C.每个人能够分别以6/7,4/5和3/4的概率击中目标.如果他们每次开火一次,他们中的两个会击中目标的几率是多少？

答案是:

P(...) = P(A)*P(B)*(1-P(C)) +
         P(B)*P(C)*(1-P(A)) +
         P(C)*P(A)*(1-P(B))
       = 27.0/70.0
       = 38.57142857142857142857142857142857142857....%

以下是我对问题的解决方案:

#include <cstdio>
#include <cctype>
#include <ctime>
#include <random>


int main()
{
   std::mt19937 engine(time(0));

   engine.discard(10000000);

   std::uniform_real_distribution<double> uniform_real(0.0,1.0);

   double prA = (6.0 / 7.0);
   double prB = (4.0 / 5.0);
   double prC = (3.0 / 4.0);

   std::size_t trails = 4000000000;
   std::size_t total_success = 0;

   for (std::size_t i = 0; i < trails; ++i)
   {
      int current_success = 0;
      if (uniform_real(engine) < prA) ++current_success;
      if (uniform_real(engine) < prB) ++current_success;
      if (uniform_real(engine) < prC) ++current_success;

      if (current_success == 2)
         ++total_success;

      double prob = (total_success * 1.0) / (i+1);

      if ((i % 1000000) == 0)
      {
         printf("%05d Pr(...) = %12.10f  error:%15.13f\n",
                i,
                prob,
                std::abs((27.0/70.0) - prob));
      }
   }

   return 0;
}

问题如下,无论我运行多少次试验,概率平线大约在0.3857002101左右.代码中有什么问题吗？

采访者表示,无论种子如何,在100万次试验中将结果收敛到大约9位小数精度是微不足道的.

关于bug在我的代码中的位置的任何想法？

更新1: 我已经尝试使用以下生成器的上述代码,它们似乎在大约同时大致试验10 ^ 9的platau.

1. 的std :: mt19937_64
2. 的std :: ranlux48_base
3. 的std :: minstd_rand0

更新2: 考虑到这个问题,我走了下面的轨道.比率27/70由27和70组成,它们都是互质的,并且在4×10 ^ 9下的因子为约57×10 ^ 6或所有数字的约1.4%.因此,在[0,4x10 ^ 9]之间随机选择的两个数字中获得27/70的"精确"比率的概率大约为1.4%(因为在4x10 ^ 9内有更多因素为27) - 所以得到确切的比例非常低,无论试验次数多少,这个数字都是恒定的.

现在,如果要谈论厚边界 - 即:70 +/5因子范围内的数字,这会增加在[0,4x10 ^ 9]范围内随机选择一对数字的概率指定/相关容差范围内的比率大约为14%左右,但是使用这种技术,与精确值相比,我们可以获得的最佳平均值大约为5位小数.这种推理方式是否正确？

Answer 1

采访者表示,无论种子如何,在100万次试验中将结果收敛到大约9位小数精度是微不足道的.

嗯,这显然是荒谬的.一百万次试验中,你无法得到千分之一的估计值.如果总数只有一个不同于理论值,那么你将减去一百万分之一,这比"小数点后9位"大一千倍.

顺便说一下,c ++ 11带有一个非常好的uniform_int_distribution函数,它实际上正确地处理了舍入:它将统一生成器的总范围分成所需范围的精确倍数和余数,并丢弃在剩余部分,因此生成的值不会因舍入而产生偏差.我对你的测试程序进行了一些修改,它在十亿次试验中确实收敛到了六位数,这与我的预期相符:

int main() {
  std::mt19937 engine(time(0));

  std::uniform_int_distribution<int> a_distr(0,6);
  std::uniform_int_distribution<int> b_distr(0,4);
  std::uniform_int_distribution<int> c_distr(0,3);

  std::size_t trials = 4000000000;
  std::size_t total_success = 0;

  for (std::size_t i = 1; i <= trials; ++i) {
    int current_success = 0;
    if (a_distr(engine)) ++current_success;
    if (b_distr(engine)) ++current_success;
    if (c_distr(engine)) ++current_success;

    if (current_success == 2) ++total_success;

    if ((i % 1000000) == 0) {
      printf("%05d Pr(...) = %12.10f  error:%15.13f\n",
             i,
             double(total_success) / i,
             std::abs((27.0/70.0) - double(total_success) / i));
    }
  }
}

返回0;

Answer 2

首先,一些基本数学表明,只有一百万次试验才能获得9个精度.鉴于我们的概率是27/70,我们可以计算出x/1000000 = 27/70哪个给出x = 385714.28571.如果我们有一个非常非常精确的均匀随机数发生器,它可以精确地生成385714个正确的试验,那么这将给出一个大约相当于abs(385714/1000000 - 0.38571428571428573) = 2.857142857304318e-079个精度位置的误差.

我不认为你的分析是正确的.给定非常准确的分布,当然可以获得所需的精度.然而, - 分布均匀性的偏斜将严重妨碍精度.如果我们进行10亿次试验,那么我们所希望的最佳精度就在于此2.85 * 10^-10.如果分布偏差甚至是100,那么这将被降低到约1 * 10^-7.我不确定大多数PRNG发行版的准确性,但问题是如何准确到这个程度.快速玩std::uniform_real_distribution<double>(0.0, 1.0),看起来肯定会有更多的变化.

Answer 3

蒙特卡罗方法倾向于缓慢收敛 - 在n次模拟之后您期望的误差与1/sqrt(n)成比例.实际上,10 ^ 9次试验后的五位精确度似乎是正确的.这里没有数字伏都教.

如果采访者正在谈论直接蒙特卡罗采样,那么......令人难以置信的是,经过一百万次试验,他可以获得九位数的准确度.