时间戳

首页 > 解决方案 > 卡在使用蒙特卡洛积分法？

c++ - 卡在使用蒙特卡洛积分法？

问题描述

我想使用蒙特卡罗集成方法，我的代码如下。如您所见，我确定了区间积分，但结果错误！这段代码有什么问题？

任何帮助将不胜感激 。

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

标签： c++montecarlo

解决方案

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

输出显示蒙特卡罗方法确实收敛：

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

推荐阅读