3

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

蒙特卡洛方法包括从一个大集合中随机选择样本实验，然后根据从这些实验结果列表中估计的概率进行推断。例如，我们可以使用 6/pi^2 是随机选择的两个整数没有共同因子的概率来近似；也就是说，它们的最大公约数将是 1。为了获得 的近似值，我们进行了大量的实验。在每个实验中，我们随机选择两个整数并执行测试以查看它们的 GCD 是否为 1。通过测试的次数为我们提供了 6/pi^2 的估计值，由此我们获得了 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

谢谢你的帮助。

4

2 回答 2

11

好吧，要直接回答您的问题，您可以将 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。这是一个伯努利过程的示例，对于多次试验n中 1 的数量m ，该过程具有二项分布。

考虑 ρ = 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，这会将 std dev 减少 10 倍，并且 π 的估计值在 3.1391 和 3.1441 之间缩小，置信度为 95%。

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

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

于 2009-08-23T16:38:06.793 回答
3

我发现我的错误。谢谢大家。我在错误的地方增加了成功的案例。

更正后的代码是：

(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)))
于 2009-08-23T16:37:05.907 回答