java - 概率模拟误差不收敛

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

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

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

标准::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 个小数位。这种推理方式正确吗？

score 8 · Accepted Answer

采访者说，无论种子如何，在 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；

score 8 · Accepted Answer

首先，一些初等数学表明，仅通过一百万次试验是不可能达到 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)，看起来肯定会有比这更多的变化。

score 8 · Accepted Answer

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

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

score 3 · Accepted Answer

因为概率是以有理数形式给出的（分母中有小整数），所以您可以将可能的情况视为尺寸为 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;
}

score 1 · Accepted Answer

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

java - 概率模拟误差不收敛

5 回答 5

Related

Reference