java - 確率シミュレーション エラーが収束しない

Question

インタビューで、私は次の問題を与えられました。最初はペン/紙を使用して解決し、次にプログラムを介して結果を検証しました。

質問は次のとおりです。

A、B、C の 3 人がいます。各人は、それぞれ 6/7、4/5、3/4 の確率で的を射ることができます。それぞれが 1 発ずつ発射した場合、そのうちの 2 発が標的に命中する確率は?

答えは次のとおりです。

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 を試してみました。

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

更新 2: 問題について考えて、私は次の道をたどりました。27 と 70 で構成される比率 27/70 は互いに素であり、4x10^9 未満の 70 の因数はおよそ 57x10^6 またはすべての数の約 1.4% です。したがって、[0,4x10^9] の間でランダムに選択された 2 つの数値から 27/70 の「正確な」比率を取得する確率は、およそ 1.4% です (4x10^9 内には 27 の因数がさらにあるため)。正確な比率は非常に低く、その数は試行回数に関係なく一定です。

ここで、厚い境界について話す場合-つまり、70 + / 5の係数の範囲の数値の場合、[0,4x10 ^ 9]の範囲内で数値のペアをランダムに選択する確率が高くなります。指定された/関連する許容範囲内の比率は約 14% ですが、この手法を使用すると、正確な値と比較した場合、平均で約 5 桁の精度が得られます。この推論の仕方は正しいでしょうか？

Answer 1

モンテカルロ法はゆっくりと収束する傾向があります---n回のシミュレーション後に予想される誤差は1/sqrt（n）に比例します。実際、10^9回の試行後の5桁の精度はほぼ正しいようです。ここでは数値的なブードゥーは起こっていません。

インタビュアーがあなたが行ったようにまっすぐなモンテカルロサンプリングについて話していた場合、100万回の試行の後に彼が9桁の精度を得ることができたのは不可能です。

Answer 2

まず、いくつかの初歩的な数学では、100 万回の試行だけでは 9 桁の精度を得ることはできないことが示されています。確率がであるとすると、どれが を与えるかを27/70計算できます。正確に 385714 回の正しい試行を生成する非常に正確な一様乱数ジェネレーターを使用した場合、これは約 の誤差をもたらし、これは必要な9 桁の精度から大きく外れます。x/1000000 = 27/70x = 385714.28571abs(385714/1000000 - 0.38571428571428573) = 2.857142857304318e-07

あなたの分析は正しくないと思います。非常に正確な分布が与えられれば、必要な精度を確実に得ることができます。ただし、分布の均一性による歪みは、精度を大幅に低下させます。10 億回の試行を行う場合、期待できる最高の精度は約2.85 * 10^-10です。分布が100でも歪んでいる場合、これは約 にノックダウンされ1 * 10^-7ます。ほとんどの PRNG ディストリビューションの正確性はわかりませんが、問題はこの程度の正確さを持つことです。で簡単に遊んでみるとstd::uniform_real_distribution<double>(0.0, 1.0)、これよりも多くの分散があることは確かです。

Answer 3

インタビュアーは、シードに関係なく、100 万回の試行で結果を小数点以下 9 桁の精度に収束させるのは簡単なことだと述べました。

まあ、それは明らかにばかげています。100万回の試行で10億分の1以内の見積もりを取得することはできません. 合計が理論値と 1 つだけ異なる場合は、100 万分の 1 ずれることになります。これは、「小数点以下 9 桁」の 1,000 倍です。

ところで、c++11 には完全に優れた uniform_int_distribution 関数が付属しており、実際には丸めを正しく処理します。均一ジェネレーターの全範囲を目的の範囲の正確な倍数と剰余に分割し、生成された値を破棄します。したがって、生成される値は丸めによってバイアスされません。私はあなたのテスト プログラムにわずかな変更を加えましたが、10 億回の試行で 6 桁に収束しました。

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 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;
}

Answer 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 ロブ = 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