math - 我在这个计算 PI 的 Scheme 程序中找不到我的错误

Question

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

The Monte Carlo method consists of choosing sample experiments at random from a large set and then making deductions on the basis of the probabilities estimated from tabulating the results of those experiments. For example, we can approximate using the fact that 6/pi^2 is the probability that two integers chosen at random will have no factors in common; that is, that their greatest common divisor will be 1. To obtain the approximation to , we perform a large number of experiments. In each experiment we choose two integers at random and perform a test to see if their GCD is 1. The fraction of times that the test is passed gives us our estimate of 6/pi^2, and from this we obtain our approximation to 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/π2。这产生“pi”= sqrt(6/(1 - 6/π2)) = sqrt(6π2/(π2-6)) = 3.911)。

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

每次试验都会给我们 0 或 1，概率为 p = 6/π2。这是一个 Bernoulli process 的例子，对于试验次数 n 中 1 的数量 m，有 binomial distribution .

考虑 ρ = m/n，即通过公约数检验的次数的分数。这个 a 具有 p 的平均值和 p(1-p)/n 的方差，或 std dev σρ = sqrt(p(1-p)/n)。对于 n = 10000，您应该期望 std dev 为 0.00488。 95% of the time you will be within 2 std devs of the mean ，并且有 5% 的时间你会在 2 个标准开发者之外，或者在 0.5982 和 0.6177 之间。因此，假设 n=10000，根据此方法估计 π 将在 3.117 和 3.167 之间 95% 的时间和 5% 的时间超出此范围。

如果您想将试验次数增加 100，这会将 std dev 减少 10 倍，并且 π 的估计值在 3.1391 和 3.1441 之间缩小，置信度为 95%。

Monte Carlo 方法非常适合粗略估计，但它们需要大量和大量试验才能获得准确的答案，并且通常会达到 yield 递减的点。

并不是说这是一种没有结果的近似 pi 的方法，只是要注意这个问题。