概率模拟误差不收敛答案

【问题标题】：Probability Simulation Error not Converging概率模拟误差不收敛
【发布时间】：2012-11-11 10:47:36
【问题描述】：

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

问题如下：

共有 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 位精度是微不足道的。

关于错误在我的代码中的任何想法？

更新 1： 我已经使用以下生成器尝试了上述代码，它们似乎都在大约同一时间大致试验 10 ^ 9。

std::mt19937_64
std::ranlux48_base
std::minstd_rand0

更新 2： 考虑到这个问题，我走上了下面的路。比率 27/70 由 27 和 70 组成，它们都是互质的，并且 4x10^9 下的 70 因子大约是 57x10^6 或所有数字的 1.4%。因此，从 [0,4x10^9] 之间随机选择的两个数字中获得 27/70 的“精确”比率的概率大约为 1.4%（因为在 4x10^9 内有更多的 27 因子） - 所以得到确切的比率非常低，无论试验次数如何，这个数字都是恒定的。

现在，如果在哪里谈论厚界限 - 即：在 70 +/5 的因子范围内的数字，这会增加在 [0,4x10^9] 范围内随机选择一对数字的概率将在指定/相关公差范围内给出大约 14% 的比率，但使用这种技术，与精确值相比，我们可以获得的最佳结果平均约为小数点后 5 位。这种推理方式正确吗？

【问题讨论】：

DV 标题不准确（你拒绝接受我的编辑）。它正在收敛，只是没有达到你想要的值。如果你问为什么它会收敛到 0.3857002101，这应该因为过于本地化而被关闭。如果您提出关于精度误差来自何处的问题，这是可以接受的。
mt19937 几乎是最快和最便宜的随机数引擎。您应该期望它不会产生完全均匀的分布，在应用二项式概率后，会产生一个小的偏差。你试过其他引擎吗？
@MikaelPersson：感谢您的建议，我尝试了以下方法：std::mt19937_64、std::ranlux48_base 和 minstd_rand0，它们似乎都给出了相同的结果。我会用我尝试过的生成器更新问题，希望 djechlin 的傀儡账户不会再对我投反对票了 :)
+1 有趣。就其价值而言，根据maths.surrey.ac.uk/hosted-sites/R.Knott/Fibonacci/cfCALC.html，您的结果与任何其他有理数都不是特别接近
根据经验，蒙特卡洛积分（这里的这个问题可以看作是一个非常简单的例子）的精度通常会随着1/sqrt(N) 而降低，其中N 是试验（参见例如 en.wikipedia.org/wiki/Monte_Carlo_integration ），因此您需要 O(10^18) 次试验以获得 10^-9 的精度。当然除非你使用解析公式，在这种情况下你在一次迭代后就完成了......

标签： java c++ algorithm math simulation

【解决方案1】：

首先，一些初等数学表明，仅通过一百万次试验是不可能达到 9 位精度的。给定我们的概率是27/70，我们可以计算出x/1000000 = 27/70，得到x = 385714.28571。如果我们有一个非常非常准确的均匀随机数生成器，它准确地生成了 385714 次正确的试验，这会给我们带来大约 abs(385714/1000000 - 0.38571428571428573) = 2.857142857304318e-07 的误差，这远远超出了想要的 9 位精度。

我不认为你的分析是正确的。给定一个非常准确的分布，当然有可能获得所需的精度。但是，分布均匀性的任何偏斜都会严重影响精度。如果我们进行 10 亿次试验，我们可以期望的最佳精度大约是 2.85 * 10^-10。如果分布偏斜甚至 100，这将被击倒到大约 1 * 10^-7。我不确定大多数 PRNG 分布的准确性，但问题在于有一个在这种程度上是准确的。快速玩一下std::uniform_real_distribution<double>(0.0, 1.0)，看起来肯定会有比这更多的变化。

【讨论】：

@Yuushi，我完全同意+1
@Potatoswatter：我认为随机游走是意料之中的事情，会不时发生。与均值等有很大的差异，但同时我们不应该期望会发生某种均值回归吗？

【解决方案2】：

面试官说，无论种子如何，在 100 万次试验中让结果收敛到小数点后 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;

【讨论】：

我喜欢你的解决方案，但我想知道，如果概率不同，它就不起作用，例如，如果 Pr(A) = 3/7 你就必须使用"
有一个问题，请您详细说明“不偏向于四舍五入”这句话，我不确定我是否理解，举个例子就好了。
@Switzy：有一些小于 1 的浮点数，而且只有有限的一组可用于 uniform_real_distribution，因为它不能偏向 0。无论 6.0/7.0 (例如）计算出不太可能实际选择这些值的 6/7；首先，因为可能值的数量可能不能被 7 整除，即 2 的幂。另一方面，uniform_int_distribution 在给定范围内均匀分布。

【解决方案3】：

Monte Carlo 方法趋向于缓慢收敛 --- 您在 n 次模拟后预期的误差与 1/sqrt(n) 成正比。实际上，经过 10^9 次试验后，五位数的准确度似乎是正确的。这里没有数字巫术。

如果面试官像你一样谈论直接的蒙特卡洛抽样，那么……经过一百万次试验后，他能获得九位数的准确度，这是令人难以置信的。

【讨论】：

【解决方案4】：

因为概率以有理数的形式给出（分母中有小整数），所以您可以将可能的情况视为尺寸为 7x5x4 的立方体（形成 140（分母的乘积）子立方体）。您可以按如下方式显式访问每个子立方体，并在 140 次迭代中得到确切的数字，而不是随机跳来跳去：

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

int main()
{
  std::size_t total_success = 0, num_trials = 0;

  for (unsigned a = 1; a <= 7; ++a)
  {
    unsigned success_a = 0;

    if (a <= 6)
      // a hits 6 out of 7 times
      success_a = 1;

    for (unsigned b = 1; b <= 5; ++b)
    {
      unsigned success_b = 0;

      if (b <= 4)
        // b hits 4 out of 5 times
        success_b = 1;

        for (unsigned c = 1; c <= 4; ++c)
        {
          unsigned success_c = 0;

          // c hits 3 out of 4 times
          if (c <= 3)
            success_c = 1;

          // count cases where exactly two of them hit
          if (success_a + success_b + success_c == 2)
            ++total_success;

          ++num_trials;

        } // loop over c
    } // loop over b
  } // loop over a

  double prob = (total_success * 1.0) / num_trials;

  printf("Pr(...) = %12.10f  error:%15.13f\n",
         prob,
         std::abs((27.0/70.0) - prob));

   return 0;
}

【讨论】：

你的回答让我很感兴趣。您是否有可能获得有关此信息的链接？
在搜索了一下之后，this link 说明了这实际上是在做什么。然而，如果概率在分母中有大整数（使得它们的乘积变得太大），或者如果概率是无理数等，这将是不可行的。
@AndreHolzner 看这个问题的好方法！ +1！我相信在使用这种技术的 SO 上有一个问题：stackoverflow.com/a/842183/1739918

【解决方案5】：

FWIW 以下 Java 似乎以您预期的速度收敛于上面的预测答案（它计算最坏情况错误的标准偏差）

import java.util.Random;
import java.security.SecureRandom;
/** from question in Stack Overflow */
public class SoProb
{
  public static void main(String[] s)
  {
    long seed = 42;


/*
In an interview, I was given the following problem to solve initially using pen/paper, then via a program to verify the result.

The question is as follows:

There are three people A,B and C. Each person is capable of hitting a target with a probability of 6/7, 4/5 and 3/4 respectively. What is the probability that if they were to each fire one shot that exactly two of them will hit the target?

The answer is:

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....%

Below is my solution to the problem:
*/

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

   Random r = new Random(seed);
   // Random r = new SecureRandom(new byte[] {(byte)seed});
   // 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);
   // double prB = (6.0 / 7.0);
   // double prC = (4.0 / 5.0);
   // double prA = (3.0 / 4.0);

   double pp = prA*prB*(1-prC) +
         prB*prC*(1-prA) +
         prC*prA*(1-prB);
   System.out.println("Pp " + pp);
   System.out.println("2870 " + (27.0 / 70.0));

   // std::size_t trails = 4000000000;
   int trails = Integer.MAX_VALUE;
   // std::size_t total_success = 0;
   int total_success = 0;

   int aCount = 0;
   int bCount = 0;
   int cCount = 0;

   int pat3 = 0; // A, B
   int pat5 = 0; // A, C
   int pat6 = 0; // B, C
   double pat3Prob = prA * prB * (1.0 - prC);
   double pat5Prob = prA * prC * (1.0 - prB);
   double pat6Prob = prC * prB * (1.0 - prA);
   System.out.println("Total pats " + 
     (pat3Prob + pat5Prob + pat6Prob));

   for (int i = 0; i < trails; ++i)
   {
      int current_success = 0;
      // if (uniform_real(engine) < prA) ++current_success;
      int pat = 0;
      if (r.nextDouble() < prA) 
      {
        ++current_success;
        aCount++;
        pat += 1;
      }
      // if (uniform_real(engine) < prB) ++current_success;
      if (r.nextDouble() < prB) 
      {
        ++current_success;
        bCount++;
        pat += 2;
      }
      // if (uniform_real(engine) < prC) ++current_success;
      if (r.nextDouble() < prC) 
      {
        ++current_success;
        cCount++;
        pat += 4;
      }
      switch (pat)
      {
        case 3:
          pat3++;
          break;
        case 5:
          pat5++;
          break;
        case 6:
          pat6++;
          break;
      }

      if (current_success == 2)
         ++total_success;

      double prob = (total_success + 1.0) / (i+2);

      if ((i % 1000000) == 0)
      {
         /*
         printf("%05d Pr(...) = %12.10f  error:%15.13f\n",
                i,
                prob,
                std::abs((27.0/70.0) - prob));
         */
         System.out.println(i + "P rob = " + prob +
           " error " +  Math.abs((27.0 / 70.0) - prob));
         Double maxVar = 0.25 / i;
         System.out.println("Max stddev " + Math.sqrt(maxVar));
         double ap = (aCount + 1.0) / (i + 2.0);
         double bp = (bCount + 1.0) / (i + 2.0);
         double cp = (cCount + 1.0) / (i + 2.0);
         System.out.println("A error " + (ap - prA));
         System.out.println("B error " + (bp - prB));
         System.out.println("C error " + (cp - prC));
         double p3Prob = (pat3 + 1.0) / (i + 2.0);
         double p5Prob = (pat5 + 1.0) / (i + 2.0);
         double p6Prob = (pat6 + 1.0) / (i + 2.0);
         System.out.println("P3 error " + (p3Prob - pat3Prob));
         System.out.println("P5 error " + (p5Prob - pat5Prob));
         System.out.println("P6 error " + (p6Prob - pat6Prob));
         System.out.println("Pats " + (pat3 + pat5 + pat6) +
           " success " + total_success);
      }
   }

  }

}

当前输出：

1099000000P rob = 0.3857148864682168 错误 6.00753931045972E-7

最大标准差 1.508242443516904E-5

一个错误-2.2208501193610175E-6

B 错误 1.4871155568862982E-5

C 错误 1.0978161945063292E-6

P3 错误 -1.4134927830977695E-7

P5 错误 -5.363291293969397E-6

P6 错误 6.1072143395513034E-6

拍拍 423900660 成功 423900660

【讨论】：