Ошибка моделирования вероятности не сходится

Ответ 1

Во-первых, какая-то элементарная математика показывает, что невозможно получить 9 точек точности только с миллионами проб. Учитывая, что наша вероятность 27/70, мы можем вычислить x/1000000 = 27/70, которая дает x = 385714.28571. Если бы у нас был очень и очень точный равномерный генератор случайных чисел, который генерировал ровно 385714 правильных испытаний, это дало бы нам погрешность приблизительно abs(385714/1000000 - 0.38571428571428573) = 2.857142857304318e-07, которая была бы значительно ниже требуемых 9 точек точности.

Я не думаю, что ваш анализ правильный. Учитывая очень точное распределение, безусловно, можно получить требуемую точность. Тем не менее, любая асимметрия от однородности в распределении серьезно затруднит точность. Если мы проведем 1 миллиард испытаний, лучшая точность, на которую мы можем надеяться, составляет около 2.85 * 10^-10. Если распределение искажено на 100, это будет сбито примерно до 1 * 10^-7. Я не уверен в точности большинства дистрибутивов PRNG, но проблема будет иметь то, что точно соответствует этой степени. Имея быструю игру с std::uniform_real_distribution<double>(0.0, 1.0), она, скорее всего, будет иметь большее отклонение от этого.

Ответ 2

Интервьюер сказал, что тривиально получить результат, чтобы сблизиться до 9 десятичных знаков в пределах 1 миллиона проб, независимо от семени.

Ну, это просто явно смешно. Вы не можете получить оценку в пределах одного из тысячи миллионов с миллионами проб. Если бы сумма была только одна, отличная от теоретического значения, вы бы отключились на один миллион, что в тысячу раз больше, чем "9 знаков после запятой".

Кстати, С++ 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));
    }
  }
}

return 0;

Ответ 3

Методы Монте-Карло имеют тенденцию сходиться медленно - ошибка, которую вы ожидаете после 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;
}

Ответ 5

FWIW следующая Java, похоже, сходится на предсказанном ответе сверху примерно на уровне, который вы ожидаете (он вычисляет стандартное отклонение ошибки наихудшего случая)

import java.util.Random;
import java.security.SecureRandom;
/** from question in Qaru */
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

Max stddev 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

Pats 423900660 успех 423900660