我在这个计算 PI 的方案程序中找不到错误

Question

我正在做蒙特卡罗实验来计算 PI 的近似值。来自 SICP：

蒙特卡罗方法包括 随机选择样本实验 从一个大集合然后制作 的基础上扣除 估计的概率 将这些结果制成表格 实验。例如，我们可以 使用以下事实进行近似 6/pi^2 是两个概率 随机选择的整数将没有 共同因素；也就是说，他们的 最大公约数为 1。 获得 的近似值，我们 进行大量实验。 在每个实验中我们选择两个 随机整数并执行测试 看看他们的 GCD 是否为 1。分数 测试通过的次数给出 我们估计 6/pi^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" = sqrt(6/(1 - 6/π²)) = sqrt(6π²/(π²- 6)) = 3.911)。

---

不过，让我们退后一步，看看蒙特卡罗方法通过此计算为我们做了什么（提示：预计收敛速度会非常慢。您运行了多少次？）...

每次试验要么给我们 0 要么1，概率 p = 6/π²。这是伯努利过程的一个示例，对于 1 的数量 m 在许多试验n中，二项式分布 。

考虑 ρ = m/n，即通过公约数测试的次数分数。该 a 的平均值为 p，方差为 p(1-p)/n，或 std dev σ_ρ = sqrt(p(1-p)/n)。对于 n = 10000，标准差应该为 0.00488。 95% 的时间您将处于平均值的 2 标准差以内，5%有时您将处于 2 标准开发之外，或介于 0.5982 和 0.6177 之间。因此，在 n=10000 的情况下，通过此方法估算的 π 95% 的时间将在 3.117 到 3.167 之间，5% 的时间在此范围之外。

如果您想将试验次数增加 100 次，那么标准差就会减少 10 倍，并且 π 的估计值会缩小到 3.1391 到 3.1441 之间，置信度为 95%。

蒙特卡罗方法非常适合粗略估计，但需要大量试验才能获得准确答案，并且通常会达到收益递减的程度。

这并不是说这是一种徒劳的近似 pi 的方法，只是要注意这个问题。

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)))