布冯针 pi 估计实验接近 pi，但未获得完美值

Question

我正在做一个非常典型的大学作业。我应该编写一个 C 程序来使用 Monte Carlo 方法通过 Buffon Needle 估计圆周率。我认为我的程序运行正常，但我从来没有得到正确的 pi。它总是接近典型的 3.14，但有时是 3,148910，有时是 3,13894。其背后的原因是什么？我可以做些什么来使它 'equal' 更准确吗？（代码注释为波兰语，() 中的注释为英语。）

#include <stdio.h>
#include <math.h>
#include <time.h>  
#include <stdlib.h>

int main(int argc, char **argv)
{
    srand(time(NULL));
    int i; //wprowadzenie licznika (counter)
    double l; // dlugosc igly (lenght of needle)
    double n;// liczba prob (number of trials)
    double L; // szerokosc miedzy liniami (width between lines)
    double x; // kat miedzy igla a normalna do linii (angle)
    double pi; // szacowana wartosc liczby pi (pi)
    double P=0; //ilość prób spełniających warunek styku igły z linią (ammount of trials that worked)
    double d; //odleglosc srodka igly od linii (distance between centre of needle and line)
    double y; //minimalizacja obliczeń 
    double  stosun; //stosunek wykonanych obliczeń (ammount of trials to ammount of succesfull trials)
    printf("Program szacuje wartosc liczby pi metoda Monte Carlo na podstawie Igly Buffona. Prosze Pamietac o:\n-Podaniu liczby prób jako wartości większej od zera i całkowitej\n-Podaniu długości igły mniejszej niż szerokości między liniami\n");
    printf("Prosze podac liczbe prob:\n");
    scanf("%lf", &n);
    printf("Prosze podac dlugosc igly:\n");
    scanf( "%lf", &l);
    printf("Prosze podac szerokość miedzy liniami:\n");
    scanf("%lf", &L);
    if (n <= 0.0 || l > L) 
    {
        printf("Nastąpił błąd wynikających z podania niepoprawnych danych\n");
        return 1;
    }
    printf("%f %f %f", n,l,L); //debuging midway

    for (i=0; i<n; i++)
{
    x = ((double)rand()/RAND_MAX)*3.1415;       // losowy kat w radianach (random angle in radians)

    d = ((double)rand()/RAND_MAX)*L/2;         // losowa dlugosc igly mniejsza niz odstep miedzy liniami (random needle lenght)

    y = (l/2) * sin (x);
printf("Próba%d\n", i);
    if (d<=y)                                    
    {
        P++;                                  
    }

}

stosun=n/P;
printf("Stosun %lf", stosun);
pi=stosun*2*l/L;
printf("liczba pi=%lf", pi); 
    return 0;
}

Answer 1

Monte Carlo 模拟的方差往往非常大。（随着 n 的增加，它们会收敛到零误差，但非常缓慢。）您可以阅读有关布冯针的文章，特别是 here。我想提请您注意这个数字，其中绘制了 5 个不同的运行s，每个 1000 次投掷。检查它们中的每一个，看起来值已经收敛，但结果却非常错误（并且在每种情况下都不同）。

最好在开始之前进行分析。上述来源中严格进行的误差估计约为 5.6 / n。这就是方差，所以为了得到标准差，我们需要一个平方根，这让情况变得更糟。投掷 1000 次后，您可以预期 π 的误差约为 0.073（尽管这很容易达到两倍或更低，纯属偶然）。经过一百万次抛掷后，它应该正确到大约 ±0.002。在一亿之后，您可能会得到正确的小数点后的前三位数字和第四位或（很少）第五位的第一个错误。我不知道你的 n 但你得到的值似乎在典型范围内。 编辑： n = 1,000,000，它们是 +3 σ 和 –1 σ.

请注意，这是一种常见的误解，认为可以通过多次运行并对它们的结果进行平均来解决这个问题（毕竟，正确值应该以相同的概率上下波动）。这当然有效，但除了明显的可以分发到多台机器之外，没有任何优势。运行每 100 次模拟一千次迭代的误差与运行一次十万次迭代的误差完全相同。

使估计更好的一个明显方法是增加命中数。 Monte Carlo 方法的一大优势是各个样本通常彼此完全独立，因此为并行化（OpenMP、矢量化指令、GPU）提供了很大的机会。但是第一步确实应该尝试使您的代码运行更快，以便将 n 增加一个数量级左右的伤害较小。例如，因为您从分析中知道精度永远不会真正令人惊叹，所以您可以使用 float 而不是 double，后者在某些平台上速度高达两倍。消除对 sin 的调用（如果您不需要 π 的值作为开始，它肯定会更优雅）。仅每 1000 步左右输出一次摘要。

一个不太明显的方法是正确地进行分析并注意如何减少方差中的分子（分母将是 n 任何你做的事情）。如果我没记错的话，这个实验运行s 在针的长度等于线之间的间隔时是最好的：小于它是次优的，大于它会带来你不想要的多个交叉点。可能还有其他方法可以通过改变角度 and/or 位置的分布来增强它。

布冯针 pi 估计实验接近 pi，但未获得完美值

Buffon needle pi estimate experiment gets close to pi, but doesnt get perfect value

c

pi