Интеграция Монте-Карло для нахождения числа Пи с определенной точностью в FORTRAN

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

У меня не возникло затруднений с написанием кода с произвольным количеством попыток:

REAL(8) FUNCTION distance(xvalue, yvalue) RESULT(dist)
IMPLICIT NONE
REAL(8), INTENT(in) :: xvalue, yvalue
dist = SQRT(xvalue**2 + yvalue**2)
END FUNCTION distance 

PROGRAM ass2
  IMPLICIT NONE

  INTEGER, DIMENSION(1) :: SEED
  REAL(8) :: p, x, y

  REAL(8), EXTERNAL :: distance

  REAL(8) :: pi_last, pi
  INTEGER :: npc, npt, i

  npc = 0
  npt = 0
  pi = 1.0

  SEED(1) = 12345
  CALL RANDOM_SEED

  DO i=1, 1000000000

     CALL RANDOM_NUMBER(p)
     x = p
     CALL RANDOM_NUMBER(p)
     y = p

     npt = npt + 1

     IF (distance(x, y) < 1.0)  THEN
          npc = npc + 1
     END IF

     pi_last = pi
     pi = 4.0*(npc*1.0)/(npt*1.0)

  END DO

  PRINT*, 'Pi:', pi

END PROGRAM ass2

Я заметил, что он сходится примерно как sqrt (N шагов). Теперь мне нужно остановить алгоритм с определенной точностью, поэтому я создал бесконечный цикл DO с EXIT внутри оператора IF:

    REAL(8) FUNCTION distance(xvalue, yvalue) RESULT(dist)
IMPLICIT NONE
REAL(8), INTENT(in) :: xvalue, yvalue
dist = SQRT(xvalue**2 + yvalue**2)
END FUNCTION distance 

PROGRAM ass2
  IMPLICIT NONE

  INTEGER, DIMENSION(1) :: SEED
  REAL(8) :: p, x, y

  REAL(8), EXTERNAL :: distance

  REAL(8) :: pi_last, pi
  INTEGER :: npc, npt, i

  npc = 0
  npt = 0
  pi = 1.0

  SEED(1) = 12345
  CALL RANDOM_SEED

  DO

     CALL RANDOM_NUMBER(p)
     x = p
     CALL RANDOM_NUMBER(p)
     y = p

     npt = npt + 1

     IF (distance(x, y) < 1.0)  THEN
          npc = npc + 1
     END IF

     pi_last = pi
     pi = 4.0*(npc*1.0)/(npt*1.0)

  IF ( ABS(pi - pi_last) < 0.000001 .AND. pi - pi_last /= 0)    THEN
    EXIT
  END IF

  END DO

  PRINT*, 'Pi:', pi

END PROGRAM ass2

Проблема в том, что это возвращает значение пи, которое не имеет той точности, о которой я просил. Я понимаю логику этого: если я получаю два последовательных значения, далеко от пи, но близко друг к другу, условие будет выполнено, и программа выйдет из оператора DO. Проблема в том, что я не знаю, как изменить его, чтобы добиться определенной мной точности. Итак, вопрос:

Как мне реализовать этот алгоритм таким образом, чтобы я мог определять точность вывода числа Пи?

РЕДАКТИРОВАТЬ: Хорошо, я реализовал оба ваших решения, и они работают, но только для 10 ^ (- 1), 10 ^ (- 3) и 10 ^ (- 5). Я думаю, что это проблема псевдослучайной последовательности, если для 10 ^ (- 2) и 10 ^ (- 4) она возвращает неправильное значение числа пи.


algorithm fortran montecarlo pi
person iacolippo    schedule 11.10.2015    source источник
comment
Без прохождения кода ваша точность жестко запрограммирована на 0,000001, достигается ли это? Обратите внимание, что эта мера относится не к Pi, а к двум последовательным итерациям. Если вы хотите, чтобы это можно было изменять во время выполнения, вам нужна переменная, и вы получите значение либо в качестве аргумента вашей команды, либо прочтите его при инициализации. Ищите read для последнего и get_command_argument для первого варианта.   -  person haraldkl    schedule 11.10.2015
comment
Я думаю, я не понял, извините, я виноват. Проблема в том, что он не обеспечивает такой жестко запрограммированной точности. Я не хочу иметь возможность изменять значение во время выполнения, также нет необходимости получать его с клавиатуры с помощью команды READ * :) Я также заметил, что это связано с двумя последовательными итерациями, но обратите внимание, что я не могу определить константа пи со значением из книг или Интернета, чтобы проверить правильность цифры. Я должен получить число пи с определенной точностью, не зная его значения. Как это сделать?   -  person iacolippo    schedule 11.10.2015
comment
Примечание: если вы удалите квадратный корень, ваша программа будет работать быстрее ;-)   -  person Anthony Scemama    schedule 11.10.2015
comment
Обычно за сходимостью MC следят, вычисляя как-то статистическую ошибку. Поскольку образцы здесь статистически независимы, для этой цели можно использовать стандартное отклонение. Кроме того, лучше использовать литералы с двойной точностью, такие как 1.0d0 ...   -  person roygvib    schedule 12.10.2015
comment
Каждый раз, когда кто-то учит новичка использовать real(8), бог убивает котенка http://stackoverflow.com/questions/838310/fortran-90-kind-parameter   -  person Vladimir F    schedule 12.10.2015
comment
LOL Я буду использовать выбранный вид, начиная с этого момента   -  person iacolippo    schedule 12.10.2015
comment
Помимо всего прочего, ваш критерий сходимости не очень хорош (если посмотреть только на последние два пункта). В этом случае, если вы знаете ответ, почему бы вам не проверить значение высокой точности для Pi? Если вы не знали результата, вам потребовались бы более сложные критерии, например отслеживание нескольких последних сотен результатов.   -  person agentp    schedule 13.10.2015
comment
Да, я на самом деле написал в вопросе, что мои критерии не подходили :) Даже предлагаемое вами решение (отслеживания последних нескольких сотен) не работает, и у меня должны быть некоторые критерии того, как выбрать это количество шагов для чек. Другое решение - проверка значения высокой точности - недопустимо, я должен делать это, не зная значения числа пи. Я выполнил свою домашнюю работу, опубликую здесь, если получу решение. Всем спасибо :)   -  person iacolippo    schedule 14.10.2015

Ответы (2)


arrow_upward
4
arrow_downward

Недостаточно указать желаемую точность - вы также должны допускать некоторую вероятность того, что цель точности не будет достигнута. Затем вы можете вычислить количество испытаний в (например) неравенстве Хёффдинга, чтобы удовлетворить желаемое точность с желаемой вероятностью (как вы заметили, n должно быть примерно квадратным корнем из единицы сверх точности, чтобы добиться успеха с постоянной вероятностью).

person David Eisenstat    schedule 11.10.2015
comment
Это имеет смысл :) Он работает для 10 ^ (- 3) и 10 ^ (- 5), но не для 10 ^ (- 4), и я думаю, что это проблема генератора псевдослучайных чисел, который дает последовательность, которая попадает в 0,1% вероятности, что цель не достигнута. Я новичок в FORTRAN, поэтому я просто принял его так, как мне его дал учитель. - person iacolippo; 11.10.2015

arrow_upward
3
arrow_downward

В идеальных условиях (идеальный генератор случайных чисел, генерирующий действительные числа в математическом смысле) ваша переменная npc является случайной величиной с биномиальное распределение B (n, / 4), где n npt из вашей программы. Его ожидаемое значение равно n * / 4, поэтому вы правильно вычислите приближение как pi=4*npc/npt. Теперь это приближение может принимать все значения от 0 до 4 независимо от того, сколько итераций цикла вы вычисляете, потому что npc может принимать все значения от 0 до npt. Для диапазона вокруг вы можете указать только вероятность (используя c как сокращение для npc; P обозначает вероятность события):

P(|pi - π| < d) = P(-d < pi - π < d) = P(-d < 4*c/n - π < d) = P(n*(π-d)/4 < c < n*(π+d)/4) =

= P(c < n*(π+d)/4) - P(c < n*(π-d)/4) ~=

~ = F N (n * (+ d) / 4) - F N (n * (- d) / 4) = 2F (d * (n / ( (4 -)))) - 1

где F N - функция вероятности нормального распределения N (n / 4; n / 4 (1- / 4)), которое аппроксимирует указанное выше биномиальное распределение, а F - функция вероятности стандартного нормального распределения. Теперь, учитывая отклонение d и вероятность p, вы можете вычислить n s.t. последний член не опускается ниже p:

n = ceil ((4 -) (F -1 ((p + 1) / 2) / d) ^ 2))

Затем с помощью n итераций цикла вы можете вычислить приближение pi с желаемой точностью с заданной вероятностью. Если мы хотим достичь вероятности p = 99%, то приведенная выше формула упрощается до

п ~ = 17,89 / д 2,

поэтому для точности d = 0,0001 необходимо примерно n = 1,789E9 итераций!

Примечание: поскольку компьютер не может соответствовать вышеуказанным идеальным настройкам, существует также (теоретический) предел точности, которого вы можете достичь с помощью этого алгоритма. На компьютере может быть представлено только конечное число чисел с плавающей запятой, поэтому ваши точки (x,y) лежат в своего рода сетке. Наилучшее приближение того, что может быть вычислено с помощью этого алгоритма, сводится к выполнению вашего цикла по всем точкам сетки в [0,1] x [0,1]. Старая добрая C-функция rand () имеет разрешение 31 бит (по крайней мере, в VS stdlib). Таким образом, нет смысла вычислять более n = 31 2 точек, что дает максимальную точность (17,89 / n) = 1,97E-9 при запросе 99% правильности.

person coproc    schedule 11.10.2015
comment
Это тоже работает, спасибо. Смотрите на редактирование, так как моя проблема теперь еще одна - person iacolippo; 12.10.2015
comment
пробовали ли вы промежуточную точность, например, 5E-3 или 3E-4 и т. д.? - person coproc; 12.10.2015