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

Question

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

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

Но когда я запускаю свою программу, я получаю значения вроде 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; ² в пределе, когда количество испытаний приближается к бесконечности. Это дает "pi" = sqrt (6 / (1 - 6 / & pi; ² )) = sqrt (6 & pi; ² / (& pi; ² - 6)) = 3,911).

---

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

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

Рассмотрим & rho; = m / n, доля проходящих тест общего делителя. Это a имеет среднее значение p и дисперсию p (1-p) / n или стандартное отклонение _{& rho;} = sqrt (p (1-p) / n). Для n = 10000 следует ожидать, что стандартное отклонение составляет 0,00488. 95% времени, когда вы будете в пределах 2 стандартных отклонений от среднего значения , и 5% времени, когда вы будете вне 2 стандартных разработок, или между 0,5982 и 0,6177. Итак, оценка & pi; из этого метода, если n = 10000, будет между 3.117 и 3.167 95% времени и вне этого диапазона 5% времени.

Если вы хотите увеличить количество испытаний на 100, это уменьшает стандартное отклонение в 10 раз и оценку & pi; сужается между 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)))