c - C错误输出中高斯函数f(x)= exp(-x ^ 2/2)的蒙特卡洛积分

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;
}

但是，我的代码输出与解决方案不匹配。我尝试调试并打印出所有xRand，yRand和yReal值以进行100次试验(并在Matlab中检查了带有特定xRand值的yReal值，以防万一我有错别字)，并且这些值似乎不在无论如何...我不知道我的错误在哪里。

＃次试验的正确输出= [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]，您的输出应该是正确的(因为它们包括平均值，因此正确的 yBound为1) 。那么，为什么您仍然得到“错误”的输出？

答案是，你不是。 您的输出是“正确的”。除此之外，蒙特卡洛过程涉及随机性，并且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到 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的数字。