Я не могу найти свою ошибку в этой программе Scheme для вычисления PI

Question

Я провожу эксперимент Монте-Карло, чтобы вычислить приближение числа ПИ.От SICP:

Метод Монте-Карло состоит из случайного выбора выборочных экспериментов из большого набора и последующего выполнения выводов на основе вероятностей, оцененных по сведению результатов этих экспериментов в таблицу .Например, мы можем приблизиться, используя тот факт, что 6 / pi ^ 2 - это вероятность того, что два выбранных случайным образом целых числа не будут иметь общих факторов;то есть, что их наибольший общий делитель будет равен 1.Чтобы получить приближение к, мы проводим большое количество экспериментов.В каждом эксперименте мы выбираем два целых числа случайным образом и выполняем тест чтобы увидеть, равен ли их GCD 1.Доля раз, когда тест пройден, дает нам нашу оценку 6 / pi ^ 2, и из этого мы получаем наше приближение к pi.

Но когда я запускаю свою программу, я получаю значения типа 3.9...

Вот моя программа:

(define (calculate-pi trials)
  (define (this-time-have-common-factors?)
    (define (get-rand)
      (+ (random 9999999999999999999999999999999) 1))
    (= (gcd (get-rand) (get-rand)) 1))
  (define (execute-experiment n-times acc)
    (if (> n-times 0)
        (if (this-time-have-common-factors?)
            (execute-experiment (- n-times 1) acc)
            (execute-experiment (- n-times 1) (+ acc 1)))
        acc))
  (define n-success (execute-experiment trials 0))
  (define prob (/ n-success trials))
  (sqrt (/ 6 prob)))

Мой интерпретатор - MIT / GNU 7.7.90

Спасибо за любую помощь.

Answer 1

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

    (if (this-time-have-common-factors?)
        (execute-experiment (- n-times 1) (+ acc 1)
        (execute-experiment (- n-times 1) acc))

Итак, вы вычисляете 1 - 6 / π² в пределе, когда количество испытаний приближается к бесконечности.Это дает "pi" = sqrt(6/(1 - 6/π²)) = sqrt(6π²/(π²-6)) = 3.911).

---

Однако давайте сделаем шаг назад и посмотрим, что метод Монте-Карло делает для нас с этим вычислением (подсказка:ожидайте очень медленной конвергенции.Сколько раз вы его запускаете?)...

Каждое испытание дает нам либо 0, либо 1, с вероятностью p = 6/π².Это пример того, как Процесс Бернулли, который имеет, для числа 1-х m в ряде судебных процессов n, а биномиальное распределение.

Рассмотрим ρ = m / n, долю времени прохождения теста с общими делителями.Это a имеет среднее значение p и дисперсию p (1-p) / n, или стандартное отклонение σ_ρ = sqrt(p(1-p)/n).При n = 10000 вам следует ожидать, что std dev равен 0.00488. В 95% случаев вы будете находиться в пределах 2 стандартных значений от среднего, и 5% времени вы будете находиться за пределами 2 std-разработчиков, или между 0,5982 и 0,6177.Таким образом, оценка π из этого метода, учитывая n = 10000, будет находиться между 3.117 и 3.167 в 95% случаев и вне этого диапазона в 5% случаев.

Если вы хотите увеличить количество испытаний на 100, это уменьшает std dev в 10 раз, и оценка π сужается между 3.1391 и 3.1441 с достоверностью 95%.

Методы Монте-Карло хороши для приблизительной оценки, но для получения точного ответа требуется МНОГО-много испытаний, и обычно доходность снижается.

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

Answer 2

Я нахожу свою ошибку.Спасибо всем.Я увеличивал количество успешных случаев не в том месте.

Исправленный код является:

(define (calculate-pi trials)
  (define (this-time-have-common-factors?)
    (define (get-rand)
      (+ (random 9999999) 1))
    (= (gcd (get-rand) (get-rand)) 1))
  (define (execute-experiment n-times acc)
    (if (> n-times 0)
        (if (this-time-have-common-factors?)
            (execute-experiment (- n-times 1) (+ acc 1))
            (execute-experiment (- n-times 1) acc))
        acc))
  (define n-success (execute-experiment trials 0))
  (define prob (/ n-success trials))
  (sqrt (/ 6 prob)))