Crude Monte-Carlo 积分在更多点上出错

Question

我正在使用这种粗略的蒙特卡洛积分技术来找出 $\pi$ 的值，并注意到随着我增加样本点的数量，积分值逐渐偏离实际值。代码是：

#include<iostream>
#include<cmath>
#include<cstdlib>

using namespace std;

float f(float x)//definition of the integrand
{
    return sqrt(1-x*x);
}

float rand1()//random number generator between 0 and 1
{
    float s=rand();
    return s/(RAND_MAX+1.0);
}

float calcint(float xi,float xf,float yi,float yf,float N)//integrator
{
    float n=0;
    for(int i=0;i<N;i++)
    {
        float x=(xf-xi)*rand1();float y=(yf-yi)*rand1();
        if (y<f(x))
        {
            n=n+1;
        }
    }
    return n/N*(xf-xi)*(yf-yi);
}

int main()
{
    float N=100000000;
    for (int i=1; i<N; i=i+N/10)//lists integration value for different sampling
    {
        cout<<i<<"\t"<<4*calcint(0,1,0,1,i)<<endl;
    }
    return 0;
}

输出是，

10000000 3.14188

20000000 3.14059

30000000 2.23696

40000000 1.67772

50000000 1.34218

60000000 1.11848

70000000 0.958698

80000000 0.838861

90000000 0.745654

为什么会这样？蒙特卡洛积分技术是否保证可以收敛到更多的样本点？

Answer 1

问题是 float 类型的精度有限。它有 24 位有效精度位，float 类型可以表示的最大可能整数是 16777216，而 16777217 不能表示，因为它需要 25 位有效位（二进制为 1 0000 0000 0000 0000 0000 0001）。在此处查看更多详细信息：Which is the first integer that an IEEE 754 float is incapable of representing exactly?

表示16777216.0f加上1.0f，结果是16777216.0f，而不是16777217.0f。因此，n应该使用整型而不是浮点型来统计事件的数量。

Crude Monte-Carlo 积分在更多点上出错

Crude Monte-Carlo integration goes wrong with more points

c++

montecarlo