Monte Carlo 集成以在 FORTRAN 中以一定的精度查找 pi
Monte Carlo integration to find pi with a certain precision in FORTRAN
我正在学习数值方法课程，我被要求实施著名的 Monte Carlo 算法来找到 pi，您可以找到 here.
我通过任意次数的试验编写代码没有困难：
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 of steps)。现在我必须以一定的精度停止算法，所以我在 IF 语句中创建了一个带有 EXIT 的无限 DO 循环：
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
问题是这个 returns 的 pi 值不符合我要求的精度。我明白了它背后的逻辑：如果我得到两个远离 pi 但彼此接近的连续值，则条件将得到满足，程序将退出 DO 语句。问题是我不知道如何修改它以获得我决定的精度。所以问题是：
如何以可以决定输出中 pi 精度的方式实现此算法？
编辑：好的，我实施了您的两个解决方案并且它们有效，但仅适用于 10^(-1)、10^(-3) 和 10^(-5)。如果对于 10^(-2) 和 10^(-4) 它 returns 是一个不正确的 pi 值，我认为这是伪随机序列的问题。
仅指定所需的精度是不够的——您还需要允许不满足精度目标的一些可能性。然后，您可以解决（例如）Hoeffding's inequality 中的试验次数，以达到所需概率的所需精度（正如您所观察到的，n 需要大约为 1 的平方根才能成功）常数概率）。
在理想情况下（生成数学意义上的实数的完美随机数生成器）您的变量 npc 是一个随机变量 binomial distribution B(n,π/4) 其中 n 是npt 来自您的程序。它的 期望值 是 n*π/4，所以你正确地计算出 π 的近似值 pi=4*npc/npt。现在这个近似值可以取从 0 到 4 的所有值，不管你计算了多少次循环迭代，因为 npc 可以取从 0 到 npt 的所有值。对于 π 附近的范围你只能给出一个概率（用 c 作为 shorthand for 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) ~=
~= FN(n*(π+d)/4) -
FN(n*(π-d)/4) = 2F(d*√(n/(π(4-π))))-1
其中 FN 是正态分布 N(nπ/4;nπ/4(1-π/4)) 的概率函数，其中 approximates the above binomial distribution 和F 是标准正态分布的概率函数。现在给定偏差 d 和概率 p，您可以计算 n s.t。最后一项不低于 p:
n = ceil(π(4-π)(F-1((p+1)/2)/d)^2))
然后通过 n 次循环迭代，您可以以给定的概率将 π 的近似值 pi 计算到所需的精度。 如果我们要达到p=99%的概率，那么上面的公式就简化为
n ~= 17.89/d2 ,
所以对于精度 d=0.0001 大约需要 n=1.789E9 次迭代！
注意：由于计算机无法满足上述理想设置，因此使用此算法可以达到的精度也有（理论上的）限制。在计算机中只能表示有限多个浮点数，因此您的点 (x,y) 位于一种网格上。可以使用此算法计算的 π 的最佳近似值归结为对 [0,1]x[0,1] 中的所有网格点执行循环。好的旧 C 函数 rand() 具有 31 位的分辨率（至少在 VS stdlib 中）。因此，计算超过 n=312 个点是没有意义的，这在要求 99% 的正确性时给出了 √(17.89/n) = 1.97E-9 的最大精度。
我正在学习数值方法课程，我被要求实施著名的 Monte Carlo 算法来找到 pi，您可以找到 here.
我通过任意次数的试验编写代码没有困难：
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 of steps)。现在我必须以一定的精度停止算法，所以我在 IF 语句中创建了一个带有 EXIT 的无限 DO 循环：
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
问题是这个 returns 的 pi 值不符合我要求的精度。我明白了它背后的逻辑：如果我得到两个远离 pi 但彼此接近的连续值，则条件将得到满足，程序将退出 DO 语句。问题是我不知道如何修改它以获得我决定的精度。所以问题是：
如何以可以决定输出中 pi 精度的方式实现此算法？
编辑：好的，我实施了您的两个解决方案并且它们有效，但仅适用于 10^(-1)、10^(-3) 和 10^(-5)。如果对于 10^(-2) 和 10^(-4) 它 returns 是一个不正确的 pi 值，我认为这是伪随机序列的问题。
仅指定所需的精度是不够的——您还需要允许不满足精度目标的一些可能性。然后，您可以解决（例如）Hoeffding's inequality 中的试验次数，以达到所需概率的所需精度（正如您所观察到的，n 需要大约为 1 的平方根才能成功）常数概率）。
在理想情况下（生成数学意义上的实数的完美随机数生成器）您的变量 npc 是一个随机变量 binomial distribution B(n,π/4) 其中 n 是npt 来自您的程序。它的 期望值 是 n*π/4，所以你正确地计算出 π 的近似值 pi=4*npc/npt。现在这个近似值可以取从 0 到 4 的所有值，不管你计算了多少次循环迭代，因为 npc 可以取从 0 到 npt 的所有值。对于 π 附近的范围你只能给出一个概率（用 c 作为 shorthand for 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) ~=
~= FN(n*(π+d)/4) - FN(n*(π-d)/4) = 2F(d*√(n/(π(4-π))))-1
其中 FN 是正态分布 N(nπ/4;nπ/4(1-π/4)) 的概率函数，其中 approximates the above binomial distribution 和F 是标准正态分布的概率函数。现在给定偏差 d 和概率 p，您可以计算 n s.t。最后一项不低于 p:
n = ceil(π(4-π)(F-1((p+1)/2)/d)^2))
然后通过 n 次循环迭代，您可以以给定的概率将 π 的近似值 pi 计算到所需的精度。 如果我们要达到p=99%的概率，那么上面的公式就简化为
n ~= 17.89/d2 ,
所以对于精度 d=0.0001 大约需要 n=1.789E9 次迭代！
注意：由于计算机无法满足上述理想设置，因此使用此算法可以达到的精度也有（理论上的）限制。在计算机中只能表示有限多个浮点数，因此您的点 (x,y) 位于一种网格上。可以使用此算法计算的 π 的最佳近似值归结为对 [0,1]x[0,1] 中的所有网格点执行循环。好的旧 C 函数 rand() 具有 31 位的分辨率（至少在 VS stdlib 中）。因此，计算超过 n=312 个点是没有意义的，这在要求 99% 的正确性时给出了 √(17.89/n) = 1.97E-9 的最大精度。