Search code examples
c++integrationmontecarlo

C++ Monte Carlo Integration: How to run code multiple times without summing the results?

This is the code I had before I changed it to include the ability to run it multiple times:

#include <iostream>
#include <fstream>
#include <stdlib.h>
#include <math.h>
#include <random>
#include <iomanip>
#include <vector>

using namespace std; 

int main () {
    double vol;
    double hit;
    int samples;
    int i, j;
    double sum;
    double pt;
    double actual_vol;
    const double PI = 2.0*atan2(1.0,0.0);
    double abs_err;
    double rel_err;

    random_device dev;
    default_random_engine e{ dev() };
    uniform_real_distribution<double> u{0.0,1.0};

    samples = 1000000 * dim;

    actual_vol = pow(PI, double(dim/2.0)) / exp(lgamma(double(dim/2.0)+1.0));

    for (i = 0; i < samples; i++) {
            sum = 0;
            for (j = 0; j < dim; j++) {
                    pt = 2*(u(e)-0.5);
                    sum += pt*pt;
            }
            if (sqrt(sum) < 1) {
                    hit += 1;
            }
    }

    vol = ( pow(2,dim) * hit ) / samples;
    abs_err = fabs( actual_vol - vol);
    rel_err = abs_err / actual_vol;

    cout << "Average volume of your sphere: " << setprecision(7) << vol << endl;
    cout << "Actual volume: " << setprecision(7) << actual_vol << endl;
    cout << "Absolute Error: " << setprecision(7) << abs_err << endl;
    cout << "Relative Error: " << setprecision(7) << rel_err << endl;
}

And I would get the correct output which looked something like this:

Average volume of your sphere: 3.140924
Actual volume: 3.141593
Absolute Error: 0.0006686536
Relative Error: 0.000212839

Now, when I change it so that I can call that function and run it multiple times, by using this code:

#include <iostream>
#include <fstream>
#include <stdlib.h>
#include <math.h>
#include <random>
#include <iomanip>
#include <vector>

using namespace std;

double monte_carlo (int dim) {
    double vol;
    double hit;
    int samples;
    int i, j;
    double sum;
    double pt;
    double actual_vol;
    const double PI = 2.0*atan2(1.0,0.0);
    double abs_err;
    double rel_err;

    random_device dev;
    default_random_engine e{ dev() };
    uniform_real_distribution<double> u{0.0,1.0};

    samples = 1000000 * dim;

    actual_vol = pow(PI, double(dim/2.0)) / exp(lgamma(double(dim/2.0)+1.0));

    for (i = 0; i < samples; i++) {
            sum = 0;
            for (j = 0; j < dim; j++) {
                    pt = 2*(u(e)-0.5);
                    sum += pt*pt;
            }
            if (sqrt(sum) < 1) {
                    hit += 1;
            }
    }

    vol = ( pow(2,dim) * hit ) / samples;
    abs_err = fabs( actual_vol - vol);
    rel_err = abs_err / actual_vol;

    cout << "Average volume of your sphere: " << setprecision(7) << vol << endl;
    cout << "Actual volume: " << setprecision(7) << actual_vol << endl;
    cout << "Absolute Error: " << setprecision(7) << abs_err << endl;
    cout << "Relative Error: " << setprecision(7) << rel_err << endl;
}

int main (int argc, char* argv[]) {

    int dim = 0;
    int runs = 0;
    int i;

    dim =  atoi(argv[1]);
    runs = atoi(argv[2]);

    for (i = 0; i < runs; i++) {
            monte_carlo(dim);
    }

    return 0;
}

I get these results, which is now summing the previous values to the current values which is not what I want:

Average volume of your sphere: 3.141764
Actual volume: 3.141593
Absolute Error: 0.0001713464
Relative Error: 5.454126e-05
Average volume of your sphere: 6.283674
Actual volume: 3.141593
Absolute Error: 3.142081
Relative Error: 1.000156
Average volume of your sphere: 9.427502
Actual volume: 3.141593
Absolute Error: 6.285909
Relative Error: 2.000867
Average volume of your sphere: 12.56937
Actual volume: 3.141593
Absolute Error: 9.427775
Relative Error: 3.000954
Average volume of your sphere: 15.71272
Actual volume: 3.141593
Absolute Error: 12.57113
Relative Error: 4.001515
Average volume of your sphere: 18.85378
Actual volume: 3.141593
Absolute Error: 15.71219
Relative Error: 5.001345
Average volume of your sphere: 21.99504
Actual volume: 3.141593
Absolute Error: 18.85345
Relative Error: 6.001239

You'll notice the first value for the average volume of sphere is around 3.14, then the second instance of it, it is now 6.28 (or double the first one), the third instance is 9.42 (roughly three times the first one), etc.

What it should be doing is running a fresh calculation each run, and the values for each of them should all be hovering at around 3.14. How do I get it to stop summing the value from the previous run?

Thank you!!

Solution

  • This is probably because you never reinitialize your variables.

    You also have a strong "old" C bias (C headers, usage of atoi, fabs...), declare your variables when you need them, and also your paths are always going to be similar because you use the same random number generator with the same seed (default constructed).

    Still, for your problem:

    double hit = 0;
double samples = 0;

    And so on.

    Also for PI, if you have boost, use its constant instead of recomputing it with a lower precision than what it can provide.