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

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

[0, 1] 上 # of trial = 100 的正确输出是 0.810，而我的是 0.880； [-1, 0] 上 # of trial = 50 的正确输出为 0.900，而我的为 0.940。谁能找到我做错的地方？非常感谢。

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

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

但是是导师提供的，他说会生成一个从0到1的随机数。为什么"rand()"后面用'/'而不是'%'？

Answer 1

您的代码没有明显的错误（尽管在上限计算中存在错误，正如@TasosPapastylianou 指出的那样，尽管这不是您的测试用例中的问题）。在 100 次试验中，您的答案 0.880 比 0.810 更接近积分的实际值 (0.855624...)，而且这些数字都与真实值相差甚远，无法表明代码中存在彻底的错误。似乎在抽样误差范围内（尽管见下文）。这是在[0,1] 上运行 1000 次蒙特卡洛积分（在 R 中完成，但使用相同算法）的直方图，其中 [0,1] 进行了 100 次试验：

除非您的讲师详细指定了算法和种子，否则您不应期望得到完全相同的答案。

至于您关于rand() / (double) RAND MAX 的第二个问题：这是一种避免modulo bias 的尝试。这种偏差可能会影响您的代码（尤其是考虑到您四舍五入到小数点后 3 位的方式），因为它似乎确实高估了积分（基于运行它十几次左右）。也许你可以在你的代码中使用它，看看你是否能得到更好的结果。

Answer 2

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

首先，为了不碍事，我们在这里讨论的是标准高斯，即

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

---

我的代码逻辑正确吗？

否。您的代码有两个逻辑错误：一个在line 29 上（即您的if 语句），一个在line 40 上（即valIntegral 的计算），这是第一个逻辑错误的直接后果。

对于第一个错误，请考虑以下情节以了解原因：

您的蒙特卡洛过程有效地考虑了某个范围内的有界框，然后说“我将在该框内随机放置点，然后计算随机落在 以下的点的总数的比例 em> 曲线；积分估计就是有界框本身的面积乘以这个比例”。

现在，如果两者都 和 位于均值的左侧（即 0），那么您的 if 语句正确地将框的上限（即 yBound）设置为 使得盒子的最高边界包含该曲线的最高部分。因此，例如，要估计范围 [-2,-1] 的积分，您可以将上限设置为 .

同样，如果两者都 和 位于均值的右侧，那么您正确地将yBound 设置为

但是，如果 ，您应该将yBound 设置为都不 也不 ，因为 0 点高于两者！。所以在这种情况下，您的 yBound 应该只是在高斯的 peak 处，即 （在您的未归一化高斯的情况下，它的值为'1'）。

因此，正确的if语句如下：

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

至于第二个逻辑错误，我们已经提到积分的值是“边界框的面积”乘以“成功率”。但是，您似乎在计算中忽略了框的 height。确实，在特殊情况下 ，非归一化高斯函数的高度默认为“1”，因此可以省略此项。我怀疑这就是它可能被错过的原因。但是，在其他两种情况下，边界框的高度必然小于 1，因此需要包含在计算中。所以line 40的正确代码应该是：

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

---

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

尽管存在上述逻辑错误，正如我们上面所讨论的，您的输出应该对于特定区间 [0,1] 和 [-1 ,0] （因为它们包括平均值，因此是正确的yBound of 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);
}

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