粗略的蒙特卡洛积分在更多点上出错答案

【问题标题】：Crude Monte-Carlo integration goes wrong with more points粗略的蒙特卡洛积分在更多点上出错
【发布时间】：2017-10-01 19:18:55
【问题描述】：

我正在使用这种粗略的蒙特卡洛积分技术来找出 $\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

为什么会这样？蒙特卡洛积分技术能否保证在大量样本点上收敛？

【问题讨论】：

首先，当 N 为 0 时，您在 calcInt 中除以 0。其次，您正在混合 float 和 int ，预计不会损失精度。
我解决了第一个问题。但我看不到我在哪里混合int和float。我只是在 for 语句的条件下比较它们。我并没有真正转换它们。您能否详细说明这可能会如何影响结果？为什么它们会稳步下降？如果在比较过程中可能发生转换问题，那么结果应该突然减小到一个值并保持在那里。
离题：if (y < f(x))可以表示为if (y * y < 1 - x * x)。
您在 1) for 循环限制 i=i+N/10 2) 您对 calcInt 的调用（查看最后一个参数）。您正在转换它们，即使不是明确的。
@PaulMcKenzie 我看到了这些转换，但它们仍然没有解释为什么集成值会下降。请参阅接受的答案，我应该将n 声明为整数，或与N 具有相同类型。否则，n 会停止增加，但 N 会增加，这会导致错误（更小）的 n/N 值。无论如何，感谢您的努力。我很感激。

标签： c++ montecarlo

【解决方案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?

这意味着当您将 1.0f 添加到 16777216.0f 时，结果将是 16777216.0f 而不是 16777217.0f。因此，n 应该使用整数类型而不是浮点类型来统计事件数。

【讨论】：

哇！你太棒了。我现在看到了问题。那太愚蠢了。是的，它解决了问题。