蒙特卡罗积分高斯函数f(x)= exp(-x ^ 2/2)在C不正确的输出中

Question

我正在写一个简短的程序来近似高斯函数f(x)= exp(-x ^ 2/2)的定积分,我的代码如下:

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

double gaussian(double x) {
    return exp((-pow(x,2))/2);
}

int main(void) {
    srand(0);
    double valIntegral, yReal = 0, xRand, yRand, yBound;
    int xMin, xMax, numTrials, countY = 0;

    do {
        printf("Please enter the number of trials (n): ");
        scanf("%d", &numTrials);
        if (numTrials < 1) {
            printf("Exiting.\n");
            return 0;
        }  
        printf("Enter the interval of integration (a b): ");
        scanf("%d %d", &xMin, &xMax);      
        while (xMin > xMax) { //keeps looping until a valid interval is entered
            printf("Invalid interval!\n");
            printf("Enter the interval of integration (a b): ");
            scanf("%d %d", &xMin, &xMax);
        }
        //check real y upper bound
        if (gaussian((double)xMax) > gaussian((double)xMin))
            yBound = gaussian((double)xMax);
        else 
            yBound = gaussian((double)xMin);
        for (int i = 0; i < numTrials; i++) {
            xRand = (rand()% ((xMax-xMin)*1000 + 1))/1000.00 + xMin; //generate random x value between xMin and xMax to 3 decimal places             
            yRand = (rand()% (int)(yBound*1000 + 1))/1000.00; //generate random y value between 0 and yBound to 3 decimal places
            yReal = gaussian(xRand);
            if (yRand < yReal) 
                countY++;
        }
        valIntegral = (xMax-xMin)*((double)countY/numTrials);
        printf("Integral of exp(-x^2/2) on [%.3lf, %.3lf] with n = %d trials is: %.3lf\n\n", (double)xMin, (double)xMax, numTrials, valIntegral);

        countY = 0; //reset countY to 0 for the next run
    } while (numTrials >= 1);

    return 0;
}

但是,我的代码输出与解决方案不匹配.我尝试调试并打印100次试验的所有xRand,yRand和yReal值(并使用Matlab检查yReal值和特定的xRand值,以防我有任何拼写错误),这些值似乎没有超出范围无论如何......我不知道我的错误在哪里.

[0,1]中试验数= 100的正确输出为0.810,我的为0.880; 试验#的正确输出= [-1,0]上的50是0.900,而我的是0.940.谁能找到我做错的地方？非常感谢.

另一个问题是,我找不到使用以下代码的参考:

double randomNumber = rand() / (double) RAND MAX;

但它是由教师提供,他说这将产生从0随机数为1.为什么他用'/',而不是'%'之后"rand()"？

Answer 1

您的代码中存在一些逻辑错误/讨论点,包括数学和编程方面.

首先,为了解决这个问题,我们在这里讨论标准高斯,即

除了,高斯的定义line 6,省略了 规范化术语.鉴于您似乎期望的产出,这似乎是故意的.很公平.但是如果你想计算实际积分,使得实际上无限范围(例如[-1000,1000])总和为1,那么你需要那个术语.

我的代码在逻辑上是否正确？

不.您的代码有两个逻辑错误:一个打开line 29(即您的if语句),另一个打开line 40(即计算valIntegral),这是第一个逻辑错误的直接后果.

对于第一个错误,请考虑以下图表以了解原因:

你的蒙特卡罗过程有效地考虑了一定范围内的有界框,然后说"我将随机将点放在这个框内,然后计算随机落在曲线下的点总数的比例;然后积分估计有界盒子本身的面积,乘以这个比例".

现在,如果两者都有 和 在均值 的左边(即0),然后你的if语句正确设置框的上限(即yBound) 这样盒子的最上面的边界包含该曲线的最高部分.因此,例如,要估计范围[-2,-1]的积分,可以将上限设置为 .

同样,如果两者都有 和 是在正确的意思,那么你正确设置yBound到

但是,如果 ,你应该设置yBound为两个 也不 ,因为0点高于两者!所以在这种情况下,你yBound应该只是处于高斯的峰值,即 (在你的非标准化高斯的情况下,这取值为'1').

因此,正确的if陈述如下:

if (xMax < 0.0)
  { yBound = gaussian((double)xMax); }
else if (xMin > 0.0)
  { yBound = gaussian((double)xMin); }
else
  { yBound = gaussian(0.0); }

至于第二个逻辑错误,我们已经提到积分的值是"边界框的面积"乘以"成功的比例".但是,您似乎忽略了计算中框的高度.确实,在特殊情况下 ,非标准化高斯函数的高度默认为"1",因此该术语可以省略.我怀疑这就是它可能被遗漏的原因.但是,在另外两种情况下,边界框的高度必然小于 1,因此需要包括在计算中.所以正确的代码line 40应该是:

valIntegral = yBound * (xMax-xMin) * (((double)countY)/numTrials);

为什么我没有得到正确的输出？

即使存在上述逻辑错误,正如我们上面所讨论的那样,对于特定的区间[0,1]和[-1,0] ,您的输出应该是正确的(因为它们包括均值,因此它们包括1 的正确性).那么为什么你仍然得到'错误'的输出？yBound

答案是,你不是.你的输出是"正确的".除此之外,蒙特卡罗过程涉及随机性,100次试验的数量不足以导致一致的结果.如果你反复运行100次试验的相同范围,你会发现每次都会得到非常不同的结果(但总的来说,它们会分布在正确的值附近).运行1000000次试验,您将看到结果变得更加精确.

那段randomNumber代码怎么了？

该rand()函数返回[0,] 范围内的整数RAND_MAX,其中 RAND_MAX特定于系统(请查看man 3 rand).

的模的方法(即%)的工作方式如下:考虑范围[-0.1,0.3].这个范围跨越0.4个单位.0.4*1000 + 1 = 401.对于0到0之间的随机数RAND_MAX,执行rand() 模数 401将始终产生[0,400]范围内的随机数.如果再将其除以1000,则得到[0,0.4]范围内的随机数.将其添加到xmin偏移量(此处为:-0.1),您将得到[-0.1,0.3]范围内的随机数.

从理论上讲,这是有道理的.然而,不幸的是,正如在这里的另一个答案中已经指出的那样,作为一种方法,它易受模偏差的影响,因为RAND_MAX不一定完全可被401整除,因此该范围的顶部部分导致RAND_MAX与其他数字相比过多地表示某些数字. .

相比之下,老师给你的方法只是说:将rand()函数的结果除以RAND_MAX.这有效地将返回的随机数标准化为范围[0,1].这是一个更直接的事情,它避免了模数偏差.

因此,我实现这个的方法是将它变成一个函数:

double randomNumber(void) {
  return rand() / (double) RAND_MAX;
}

然后,这也简化了您的计算如下:

xRand = randomNumber() * (xMax-xMin) + xMin;
yRand = randomNumber() * yBound;

如果使用标准化高斯,即可以看到这是一个更准确的事情

double gaussian(double x) {
  return exp((-pow(x,2.0))/2.0) / sqrt(2.0 * M_PI);
}

然后比较两种方法.您将看到randomNumber()"有效无限"范围的方法(例如[-1000,1000])给出正确的结果1,而模数方法倾向于给出大于1的数字.