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

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

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

Кстати, c ++ 11 поставляется с очень хорошей функциейiform_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;

Методы Монте-Карло имеют тенденцию сходиться медленно — ошибка, которую вы ожидаете после n симуляций, пропорциональна 1 / sqrt (n). На самом деле пять цифр точности после 10-9 испытаний кажутся правильными. Здесь нет никакого числового вуду.

Если интервьюер говорил о прямой выборке по методу Монте-Карло, как вы сделали, это … неправдоподобно, что он мог получить девять цифр точности после миллиона испытаний.

поскольку вероятности даны в виде рациональных чисел (с маленькими целыми числами в знаменателе), вы можете рассматривать возможные ситуации как куб измерений 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;
}

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

Макс stddev 1.508242443516904E-5

Ошибка -2.2208501193610175E-6

Ошибка B 1.4871155568862982E-5

Ошибка С 1.0978161945063292E-6

Ошибка P3 -1.4134927830977695E-7

Ошибка P5 -5.363291293969397E-6

Ошибка P6 6.1072143395513034E-6

Пэтс 423900660 успех 423900660