Расчет π с использованием ограничений моделирования Монте-Карло - PullRequest
Новая
       40

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

0 голосов
/ 12 сентября 2018

Я задал вопрос, очень похожий на этот, поэтому в конце я упомяну предыдущие решения, у меня есть веб-сайт , который вычисляет π с ЦП клиента при сохранении его на сервере, так что пока я получил:

'701.766.448.388' указывает внутри круга, и в сумме '893.547.800.000', эти цифры рассчитываются с использованием этого кода. (рабочий пример: https://jsfiddle.net/d47zwvh5/2/)

let inside = 0;
let size = 500;

for (let i = 0; i < iterations; i++) {
  var Xpos = Math.random() * size;
  var Ypos = Math.random() * size;

  var dist = Math.hypot(Xpos - size / 2, Ypos - size / 2);

  if (dist < size / 2) {
    inside++;
  }
}

Проблема

(4 * 701.766.448.388) / 893.547.800.000 = 3,141483638

Это результат, который мы получаем, и он верен до тех пор, пока четвертая цифра, 4 должна быть 5.

Предыдущие проблемы:

  1. Я испортил расчет расстояния.
  2. Я поместил кружок от 0 ... 499, который должен быть 0 ... 500
  3. Я не использовал float, что уменьшило «разрешение»

Юридическая информация

Может быть, я достиг предела, но эта демонстрация использовала 1 миллион очков и получила 3,16. учитывая, что у меня есть около 900 миллиардов, я думаю, что это может быть более точно.

Я понимаю, что если я хочу вычислить π, это неправильный способ, но я просто хочу убедиться, что все правильно, поэтому я надеялся, что кто-то может определить что-то не так, или мне просто нужно больше «точек».

РЕДАКТИРОВАТЬ: Есть довольно много упоминаний о том, насколько нереалистично, где числа, эти упоминания, где правильно, и я теперь обновил их, чтобы они были правильными.

Ответы [ 2 ]

0 голосов
/ 14 сентября 2018

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

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

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

0 голосов
/ 13 сентября 2018

любая операция 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.

...