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 回答