stuck in use of monte-carlo integration method?

Question

I wanna use monte-carlo integration method, and my code is below. As u can see i determined the interval integration but the result is wrong ! Whats wrong with this code ?

any help will be appreciated .

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

What are you doing with the return value of `rand`? The first thing you have to get is a [0, 1] range via the divison you already have there. But better use `std::uniform_real_distribution`. And why are you putting the operands closer to `+` than `/`? That makes people confused - division comes always before, no matter the spaces. Lastly, why did you not use the function parameters in those places? — LogicStuff, Jun 24 '19 at 22:22
It might help to read [ask] and [mcve]. "It gives the wrong answer" doesn't tell us much, and doesn't show any effort by you to debug the code to narrow down the problem. Note: If you want to give us more details, you should edit your question instead of adding them as comments. — Dave S, Jun 24 '19 at 22:24
Take a look at this code: https://gist.github.com/NathanEpstein/f941a7f12d8630d38c72 — bgaard, Jun 24 '19 at 22:28
And you are not using your input arguments. You are using hard coded values instead of xmin, xmax, ymin, and ymax. — bgaard, Jun 24 '19 at 22:32

Alecto Irene Perez · Accepted Answer · 2019-06-24T23:00:51.150

The Monte_Carlo function is making things more complicated then they need to be. For integrating a 1-dimensional function, all we have to do is sample the value of the function a bunch of times within the region we're integrating over:

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

If we run this code with N = 1000, we get an integral value of 9.20569, which is pretty close to the exact answer (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;
}

We can also try multiple values of N, to have the program show how it converges. The following program calculates the integral with N being powers of 2 from 0 to 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
    }
}

The output shows that the monte carlo method does actually converge:

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

stuck in use of monte-carlo integration method?

1 Answers1