Monte Carlo 高斯函数 f(x) = exp(-x^2/2) 在 C 中的积分不正确的输出
Monte Carlo integration of the Gaussian function f(x) = exp(-x^2/2) in C incorrect output
我正在写一个小程序来逼近高斯函数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 trials = 100 的正确输出是 0.810，我的是 0.880； [-1, 0] 上的 # of trials = 50 的正确输出是 0.900，而我的是 0.940。谁能找到我哪里做错了？非常感谢。
另一个问题是，我找不到使用以下代码的参考：
double randomNumber = rand() / (double) RAND MAX;
可是导师提供的，他说会生成一个0到1的随机数，为什么他在"rand()"之后用'/'而不是'%'呢？
您的代码没有明显的错误（尽管正如@TasosPapastylianou 指出的那样，是上限计算中的错误，尽管这不是您的测试用例中的问题).在 100 次试验中，您的答案 0.880 比 0.810 更接近积分的实际值 (0.855624...)，而且这两个数字都离真实值不远，不足以表明代码中存在彻底的错误。似乎在抽样误差范围内（尽管见下文）。这是 [0,1] 上 e^(-x^2/2) 的 Monte Carlo 积分（在 R 中完成，但使用相同的算法）的 1000 次运行的直方图，有 100 次试验：
除非您的讲师非常详细地指定了算法和种子，否则您不应期望得到完全相同的答案。
至于你关于rand() / (double) RAND MAX的第二个问题：这是试图避免modulo bias。这种偏差可能会影响您的代码（特别是考虑到您四舍五入到小数点后 3 位的方式），因为它似乎确实高估了积分（基于 运行 十几次左右）。也许您可以在您的代码中使用它，看看您是否获得更好的结果。
您的代码中存在一些逻辑错误/讨论点，无论是在数学方面还是在编程方面。
首先，为了不碍事，我们这里讨论的是标准高斯，即
除了line 6上的高斯定义，省略了
规范化术语。鉴于您似乎期望的输出，这似乎是故意完成的。很公平。但是如果你想计算 actual 积分，这样 practically infinite range (e.g. [-1000, 1000]) would sum up to 1 ，那么您将需要该术语。
我的代码逻辑正确吗？
否。您的代码有两个逻辑错误：一个关于 line 29（即您的 if 语句），一个关于 line 40（即 valIntegral 的计算），这是直接后果第一个逻辑错误。
对于第一个错误，请考虑以下情节以了解原因：
您的 Monte Carlo 过程有效地考虑了一定范围内的边界框，然后说 "I will randomly place points inside this box, and then count the proportion of the total number of points that randomly fell under the curve; the integral estimate is then the area of the bounded box itself, times this proportion"。
现在，如果两者
和
均值 left（即 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); }
关于第二个逻辑错误，我们已经提到积分的值是"area of the bounding box"乘以"proportion of successes"。但是，您似乎在计算中忽略了框的 height。的确，在特殊情况下
，您的非标准化高斯函数的高度默认为“1”，因此可以省略该术语。我怀疑这就是它可能被遗漏的原因。然而，在另外两种情况下，边界框的高度必然小于小于1，因此需要计入计算。所以 line 40 的正确代码应该是：
valIntegral = yBound * (xMax-xMin) * (((double)countY)/numTrials);
为什么我没有得到正确的输出？
即使存在上述逻辑错误，正如我们上面所讨论的，您的输出应该对于特定间隔[ 0,1] 和 [-1,0]（因为它们包含均值，因此包含 1 的正确 yBound）。那你为什么仍然得到'wrong'输出？
答案是，你不是。你的输出是"correct"。除了，Monte Carlo 过程涉及随机性，100 次试验不足以产生一致的结果。如果您 运行 一次又一次地对相同范围进行 100 次试验，您会发现每次都会得到截然不同的结果（尽管总的来说，它们会围绕正确的值分布）。 运行 进行 1000000 次试验，您会发现结果变得更加精确。
那个 randomNumber 代码是怎么回事？
rand() 函数 returns 一个 整数 在 [0, RAND_MAX] 范围内，其中 RAND_MAX 是系统-特定的（看看 man 3 rand）。
modulo 方法（即 %）的工作原理如下：考虑范围 [-0.1, 0.3]。该范围跨越 0.4 个单位。 0.4 * 1000 + 1 = 401。对于从 0 到 RAND_MAX 的随机数，执行 rand() modulo 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;
你可以看到，如果你使用 normalised 高斯分布，即
，这是一个更准确的事情
double gaussian(double x) {
return exp((-pow(x,2.0))/2.0) / sqrt(2.0 * M_PI);
}
然后比较两种方法。您会看到 "effectively infinite" 范围（例如 [-1000,1000]）的 randomNumber() 方法给出了正确的结果 1，而取模方法倾向于给出大于 1 的数字。
我正在写一个小程序来逼近高斯函数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 trials = 100 的正确输出是 0.810，我的是 0.880； [-1, 0] 上的 # of trials = 50 的正确输出是 0.900，而我的是 0.940。谁能找到我哪里做错了？非常感谢。
另一个问题是，我找不到使用以下代码的参考：
double randomNumber = rand() / (double) RAND MAX;
可是导师提供的，他说会生成一个0到1的随机数，为什么他在"rand()"之后用'/'而不是'%'呢？
您的代码没有明显的错误（尽管正如@TasosPapastylianou 指出的那样，是上限计算中的错误，尽管这不是您的测试用例中的问题).在 100 次试验中，您的答案 0.880 比 0.810 更接近积分的实际值 (0.855624...)，而且这两个数字都离真实值不远，不足以表明代码中存在彻底的错误。似乎在抽样误差范围内（尽管见下文）。这是 [0,1] 上 e^(-x^2/2) 的 Monte Carlo 积分（在 R 中完成，但使用相同的算法）的 1000 次运行的直方图，有 100 次试验：
除非您的讲师非常详细地指定了算法和种子，否则您不应期望得到完全相同的答案。
至于你关于rand() / (double) RAND MAX的第二个问题：这是试图避免modulo bias。这种偏差可能会影响您的代码（特别是考虑到您四舍五入到小数点后 3 位的方式），因为它似乎确实高估了积分（基于 运行 十几次左右）。也许您可以在您的代码中使用它，看看您是否获得更好的结果。
您的代码中存在一些逻辑错误/讨论点，无论是在数学方面还是在编程方面。
首先，为了不碍事，我们这里讨论的是标准高斯，即
除了line 6上的高斯定义，省略了
我的代码逻辑正确吗？
否。您的代码有两个逻辑错误：一个关于 line 29（即您的 if 语句），一个关于 line 40（即 valIntegral 的计算），这是直接后果第一个逻辑错误。
对于第一个错误，请考虑以下情节以了解原因：
您的 Monte Carlo 过程有效地考虑了一定范围内的边界框，然后说 "I will randomly place points inside this box, and then count the proportion of the total number of points that randomly fell under the curve; the integral estimate is then the area of the bounded box itself, times this proportion"。
现在，如果两者
if 语句正确地将框的上限（即 yBound）设置为
同样，如果两者
yBound 设置为
不过,如果
yBound 设置为 两者都不
yBound 应该只是处于高斯分布的 peak，即
因此，正确的if说法如下：
if (xMax < 0.0)
{ yBound = gaussian((double)xMax); }
else if (xMin > 0.0)
{ yBound = gaussian((double)xMin); }
else
{ yBound = gaussian(0.0); }
关于第二个逻辑错误，我们已经提到积分的值是"area of the bounding box"乘以"proportion of successes"。但是，您似乎在计算中忽略了框的 height。的确，在特殊情况下
line 40 的正确代码应该是：
valIntegral = yBound * (xMax-xMin) * (((double)countY)/numTrials);
为什么我没有得到正确的输出？
即使存在上述逻辑错误，正如我们上面所讨论的，您的输出应该对于特定间隔[ 0,1] 和 [-1,0]（因为它们包含均值，因此包含 1 的正确 yBound）。那你为什么仍然得到'wrong'输出？
答案是，你不是。你的输出是"correct"。除了，Monte Carlo 过程涉及随机性，100 次试验不足以产生一致的结果。如果您 运行 一次又一次地对相同范围进行 100 次试验，您会发现每次都会得到截然不同的结果（尽管总的来说，它们会围绕正确的值分布）。 运行 进行 1000000 次试验，您会发现结果变得更加精确。
那个 randomNumber 代码是怎么回事？
rand() 函数 returns 一个 整数 在 [0, RAND_MAX] 范围内，其中 RAND_MAX 是系统-特定的（看看 man 3 rand）。
modulo 方法（即 %）的工作原理如下：考虑范围 [-0.1, 0.3]。该范围跨越 0.4 个单位。 0.4 * 1000 + 1 = 401。对于从 0 到 RAND_MAX 的随机数，执行 rand() modulo 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;
你可以看到，如果你使用 normalised 高斯分布，即
，这是一个更准确的事情double gaussian(double x) {
return exp((-pow(x,2.0))/2.0) / sqrt(2.0 * M_PI);
}
然后比较两种方法。您会看到 "effectively infinite" 范围（例如 [-1000,1000]）的 randomNumber() 方法给出了正确的结果 1，而取模方法倾向于给出大于 1 的数字。