坚持使用蒙特卡洛积分方法?
stuck in use of monte-carlo integration method?
我想用蒙特卡洛积分法,我的代码如下。如您所见,我确定了区间积分,但结果是错误的!这段代码有什么问题?
任何帮助将不胜感激。
#include <iostream>
#include <math.h>
#include <stdlib.h>
#define N 500
using namespace std;
double Func(double x) { return pow(x, 2) + 1; }
double Monte_Carlo(double Func(double), double xmin, double xmax, double ymin,
double ymax)
{
int acc = 0;
int tot = 0;
for (int count = 0; count < N; count++)
{
double x0 = (double)rand() / 4 + (-2);
double y0 = (double)rand() / 4 + 0;
float x = x0 / (float)RAND_MAX;
float y = y0 / (float)RAND_MAX;
cout << x << endl;
if (y <= Func(x))
acc++;
tot++;
// cout << "Dgage" << tot << '\t' << acc << endl;
}
double Coeff = acc / N;
return (xmax - xmin) * (1.2 * Func(xmax)) * Coeff;
}
int main()
{
cout << "Integral value is: " << Monte_Carlo(Func, -2, 2, 0, 4) << endl;
system("pause");
return 0;
}
Monte_Carlo 函数使事情变得比他们需要的更复杂。为了积分一维函数,我们所要做的就是在我们积分的区域内对函数的值进行多次采样:
#include <random>
double Monte_Carlo(double Func(double), double xmin, double xmax, int N)
{
// This is the distribution we're using to generate inputs
auto x_dist = std::uniform_real_distribution<>(xmin, xmax);
// This is the random number generator itself
auto rng = std::default_random_engine();
// Calculate the total of N random samples
double total = 0.0;
for(int i = 0; i < N; i++) {
double x = x_dist(rng); // Generate a value
total += Func(x);
}
// Return the size of the interval times the total,
// divided by the number of samples
return (xmax - xmin) * total / N;
}
如果我们 运行 这段代码加上 N = 1000
,我们得到一个整数值 9.20569
,这非常接近准确答案 (9.33333...
)。
// It's much more efficent to use x*x instead of pow
double func(double x) { return x * x + 1; }
int main()
{
cout << "Integral value is: " << Monte_Carlo(func, -2, 2, 1000) << endl;
getchar(); // Pause until the user presses enter
return 0;
}
我们也可以尝试 N
的多个值,让程序显示它是如何收敛的。以下程序计算积分,其中 N 是 2
的幂,从 0
到 30
#include <iostream>
#include <cmath>
#include <random>
using namespace std;
double func(double x) { return x*x + 1; }
double Monte_Carlo(double Func(double), double xmin, double xmax, int N) {
auto x_dist = std::uniform_real_distribution<>(xmin, xmax);
auto rng = std::default_random_engine();
double total = 0.0;
for(int i = 0; i < N; i++) {
double x = x_dist(rng); // Generate a value
total += Func(x);
}
return (xmax - xmin) * total / N;
}
int main() {
int N = 1;
for(int i = 0; i < 31; i++) {
std::cout << "N = " << N << "\t\tintegral = " << Monte_Carlo(func, -2, 2, N) << endl;
N *= 2; // Double N
}
}
输出显示 monte carlo 方法确实收敛:
N = 1 integral = 12.6889
N = 2 integral = 8.39917
N = 4 integral = 7.97521
N = 8 integral = 9.24233
N = 16 integral = 9.75632
N = 32 integral = 9.87064
N = 64 integral = 9.46945
N = 128 integral = 9.27281
N = 256 integral = 9.27395
N = 512 integral = 9.17546
N = 1024 integral = 9.19097
N = 2048 integral = 9.26203
N = 4096 integral = 9.37979
N = 8192 integral = 9.36167
N = 16384 integral = 9.28918
N = 32768 integral = 9.29766
N = 65536 integral = 9.31101
N = 131072 integral = 9.3227
N = 262144 integral = 9.32588
N = 524288 integral = 9.32805
N = 1048576 integral = 9.32726
N = 2097152 integral = 9.32722
N = 4194304 integral = 9.331
N = 8388608 integral = 9.33082
N = 16777216 integral = 9.33174
N = 33554432 integral = 9.33164
N = 67108864 integral = 9.33303
N = 134217728 integral = 9.33283
N = 268435456 integral = 9.33327
N = 536870912 integral = 9.33325
N = 1073741824 integral = 9.33333
我想用蒙特卡洛积分法,我的代码如下。如您所见,我确定了区间积分,但结果是错误的!这段代码有什么问题?
任何帮助将不胜感激。
#include <iostream>
#include <math.h>
#include <stdlib.h>
#define N 500
using namespace std;
double Func(double x) { return pow(x, 2) + 1; }
double Monte_Carlo(double Func(double), double xmin, double xmax, double ymin,
double ymax)
{
int acc = 0;
int tot = 0;
for (int count = 0; count < N; count++)
{
double x0 = (double)rand() / 4 + (-2);
double y0 = (double)rand() / 4 + 0;
float x = x0 / (float)RAND_MAX;
float y = y0 / (float)RAND_MAX;
cout << x << endl;
if (y <= Func(x))
acc++;
tot++;
// cout << "Dgage" << tot << '\t' << acc << endl;
}
double Coeff = acc / N;
return (xmax - xmin) * (1.2 * Func(xmax)) * Coeff;
}
int main()
{
cout << "Integral value is: " << Monte_Carlo(Func, -2, 2, 0, 4) << endl;
system("pause");
return 0;
}
Monte_Carlo 函数使事情变得比他们需要的更复杂。为了积分一维函数,我们所要做的就是在我们积分的区域内对函数的值进行多次采样:
#include <random>
double Monte_Carlo(double Func(double), double xmin, double xmax, int N)
{
// This is the distribution we're using to generate inputs
auto x_dist = std::uniform_real_distribution<>(xmin, xmax);
// This is the random number generator itself
auto rng = std::default_random_engine();
// Calculate the total of N random samples
double total = 0.0;
for(int i = 0; i < N; i++) {
double x = x_dist(rng); // Generate a value
total += Func(x);
}
// Return the size of the interval times the total,
// divided by the number of samples
return (xmax - xmin) * total / N;
}
如果我们 运行 这段代码加上 N = 1000
,我们得到一个整数值 9.20569
,这非常接近准确答案 (9.33333...
)。
// It's much more efficent to use x*x instead of pow
double func(double x) { return x * x + 1; }
int main()
{
cout << "Integral value is: " << Monte_Carlo(func, -2, 2, 1000) << endl;
getchar(); // Pause until the user presses enter
return 0;
}
我们也可以尝试 N
的多个值,让程序显示它是如何收敛的。以下程序计算积分,其中 N 是 2
的幂,从 0
到 30
#include <iostream>
#include <cmath>
#include <random>
using namespace std;
double func(double x) { return x*x + 1; }
double Monte_Carlo(double Func(double), double xmin, double xmax, int N) {
auto x_dist = std::uniform_real_distribution<>(xmin, xmax);
auto rng = std::default_random_engine();
double total = 0.0;
for(int i = 0; i < N; i++) {
double x = x_dist(rng); // Generate a value
total += Func(x);
}
return (xmax - xmin) * total / N;
}
int main() {
int N = 1;
for(int i = 0; i < 31; i++) {
std::cout << "N = " << N << "\t\tintegral = " << Monte_Carlo(func, -2, 2, N) << endl;
N *= 2; // Double N
}
}
输出显示 monte carlo 方法确实收敛:
N = 1 integral = 12.6889
N = 2 integral = 8.39917
N = 4 integral = 7.97521
N = 8 integral = 9.24233
N = 16 integral = 9.75632
N = 32 integral = 9.87064
N = 64 integral = 9.46945
N = 128 integral = 9.27281
N = 256 integral = 9.27395
N = 512 integral = 9.17546
N = 1024 integral = 9.19097
N = 2048 integral = 9.26203
N = 4096 integral = 9.37979
N = 8192 integral = 9.36167
N = 16384 integral = 9.28918
N = 32768 integral = 9.29766
N = 65536 integral = 9.31101
N = 131072 integral = 9.3227
N = 262144 integral = 9.32588
N = 524288 integral = 9.32805
N = 1048576 integral = 9.32726
N = 2097152 integral = 9.32722
N = 4194304 integral = 9.331
N = 8388608 integral = 9.33082
N = 16777216 integral = 9.33174
N = 33554432 integral = 9.33164
N = 67108864 integral = 9.33303
N = 134217728 integral = 9.33283
N = 268435456 integral = 9.33327
N = 536870912 integral = 9.33325
N = 1073741824 integral = 9.33333