Вычисление π с использованием ограничений моделирования Монте-Карло

Вы можете легко оценить, какую ошибку (планки погрешностей) вы должны получить, в этом вся прелесть Монте-Карло. Для этого вам нужно вычислить второй импульс и оценить дисперсию и стандартное отклонение. Хорошо, что собранное значение будет таким же, как и среднее, потому что вы просто сложили 1 после 1 после 1.

Затем вы можете получить оценку сигмы моделирования и планки погрешностей для желаемого значения. Извините, я недостаточно знаю Javascript, поэтому код здесь на C#:

using System;

namespace Pi
{
    class Program
    {
        static void Main(string[] args)
        {
            ulong N = 1_000_000_000UL; // number of samples
            var rng = new Random(312345); // RNG

            ulong v  = 0UL; // collecting mean values here
            ulong v2 = 0UL; // collecting squares, should be the same as mean
            for (ulong k = 0; k != N; ++k) {
                double x = rng.NextDouble();
                double y = rng.NextDouble();

                var r = (x * x + y * y < 1.0) ? 1UL : 0UL;

                v  += r;
                v2 += r * r;
            }

            var mean = (double)v / (double)N;
            var varc = ((double)v2 / (double)N - mean * mean ) * ((double)N/(N-1UL)); // variance
            var stdd = Math.Sqrt(varc); // std.dev, should be sqrt(Pi/4 (1-Pi/4))
            var errr = stdd / Math.Sqrt(N);

            Console.WriteLine($"Mean = {mean}, StdDev = {stdd}, Err = {errr}");

            mean *= 4.0;
            errr *= 4.0;

            Console.WriteLine($"PI (1 sigma) = {mean - 1.0 * errr}...{mean + 1.0 * errr}");
            Console.WriteLine($"PI (2 sigma) = {mean - 2.0 * errr}...{mean + 2.0 * errr}");
            Console.WriteLine($"PI (3 sigma) = {mean - 3.0 * errr}...{mean + 3.0 * errr}");
        }
    }
}

После 10⁹ образцов я получил

Mean = 0.785405665, StdDev = 0.410540627166729, Err = 1.29824345388086E-05
PI (1 sigma) = 3.14157073026184...3.14167458973816
PI (2 sigma) = 3.14151880052369...3.14172651947631
PI (3 sigma) = 3.14146687078553...3.14177844921447

что выглядит примерно правильно. Легко видеть, что в идеальном случае дисперсия была бы равна (Pi/4)*(1-Pi/4). На самом деле нет необходимости вычислять v2, просто установите его в v после симуляции.

Я, честно говоря, не знаю, почему вы получаете не то, что ожидали. Потеря точности при суммировании может быть ответом, или, как я подозреваю, ваше моделирование не создает независимых выборок из-за заполнения и перекрывающихся последовательностей (поэтому фактическое N намного ниже 900 триллионов).

Но с помощью этого метода вы контролируете ошибку и проверяете, как идут вычисления.

ОБНОВИТЬ

Я подставил ваши цифры, чтобы показать, что вы явно недооцениваете значение. Код

    N  = 893_547_800_000UL;
    v  = 701_766_448_388UL;
    v2 = v;

    var mean = (double)v / (double)N;
    var varc = ((double)v2 / (double)N - mean * mean ) * ((double)N/(N-1UL)); 
    var stdd = Math.Sqrt(varc); // should be sqrt(Pi/4 (1-Pi/4))
    var errr = stdd / Math.Sqrt(N);

    Console.WriteLine($"Mean = {mean}, StdDev = {stdd}, Err = {errr}");

    mean *= 4.0;
    errr *= 4.0;

    Console.WriteLine($"PI (1 sigma) = {mean - 1.0 * errr}...{mean + 1.0 * errr}");
    Console.WriteLine($"PI (2 sigma) = {mean - 2.0 * errr}...{mean + 2.0 * errr}");
    Console.WriteLine($"PI (3 sigma) = {mean - 3.0 * errr}...{mean + 3.0 * errr}");

И вывод

Mean = 0.785370909522692, StdDev = 0.410564786603016, Err = 4.34332975349809E-07
PI (1 sigma) = 3.14148190075886...3.14148537542267
PI (2 sigma) = 3.14148016342696...3.14148711275457
PI (3 sigma) = 3.14147842609506...3.14148885008647

Итак, у вас явно есть проблема (код? Потеря точности в представлении? Потеря точности при суммировании? Повторная/независимая выборка?)

любая операция FPU снизит вашу точность. Почему бы не сделать что-то вроде этого:

let inside = 0;
for (let i = 0; i < iterations; i++)
  {
  var X = Math.random();
  var Y = Math.random();
  if ( X*X + Y*Y <= 1.0 ) inside+=4;
  }

если мы исследуем первый квадрант единичного круга, нам не нужно изменять динамический диапазон на size, а также мы можем проверить расстояния в форме, рассчитанной на 2, которые избавляются от sqrt. Эти изменения должны повысить точность, а также скорость.

Не кодировщик JAVASCRIPT, поэтому я не знаю, какие типы данных вы используете, но вам нужно быть уверенным, что вы не нарушаете его точность. В таком случае вам нужно добавить больше переменных счетчика, чтобы уменьшить нагрузку на него. Для получения дополнительной информации см.: [edit1] точность интеграции.

Поскольку ваши числа довольно большие, держу пари, вы уже пересекли границу (не должно быть дробной части, а конечные нули также подозрительны). Например, 32-битный float может хранить только целые числа до

2^23 = 8388608

и ваш 698,565,481,000,000 намного выше этого, поэтому даже операция ++ с такой переменной приведет к потере точности, а когда показатель степени слишком велик, он даже перестанет добавлять...

Для целых чисел это не проблема, но как только вы пересекаете границу в зависимости от внутреннего формата, значение оборачивается вокруг нуля или отрицается ... Но я сомневаюсь, что это так, поскольку тогда результат будет далек от PI.