MonteCarlo Integral Program bad results

[jgialis]

First shot at coding MCI. I can't figure out how to generate random floating point integers which would probably increase the accuracy of my results, so for now I am just generating random integer ordered pairs.


Still, my results are way off even with a large sample size. Where am I going wrong? Is it a conceptual issue? Or a programming logic error?

  #include <iostream>
#include <ctime>
#include <cmath>
using namespace std;



int main()
{
    int width = 25;
    int height = 5;

    double tP = 0;
    double iP = 0;
    
    int sampleSize = 10000;
    
    int randX = 0;
    int randY = 0;
    
    double functionY = 0;
    double value = 0;
  

    srand(time(0));
    for (int i = 0; i < sampleSize; i++) {
        
        randX = (rand() % width) + 1;
        randY = (rand() % height + 1);
        functionY = pow(randX, 0.5);
        
        if(randY <= functionY){
            tP++;
            iP++;
        }
        else
            tP++;
    }
    
    cout << "tP: " << tP << endl;
    cout << "iP: " << iP << endl;
    
    
    value = (width * height) * ((iP)/(iP + tP));
    cout << "Sample size: " << sampleSize << endl;
    cout << "Actual integral: " << (0.667 * pow(25, 1.5) - (0.667 * pow(0, 1.5))); //Please ignore bad programmers practice of hardcoding 
    cout << "\nMonteCarlo Integral Value: " << value << endl;
    
}

[lastchance]

(1) You are trying to approximate an integral ... but all your test points are integers (lines 28 and 29). This will seriously bias whether your y values lie under the curve or not, which is what you are trying to find the fraction of. Consider using the uniform real-number distribution in <random>.

(2) You don't need tP (it is the same as the sample size).

(3) Use sqrt rather than pow(.,0.5). Don't bother to find 0 to any power. Use 2.0/3 rather than the inaccurate 0.667.

(4) If you want to make your program more useful, put functionY in a separate function, rather than hard-coding it in int main().

[jgialis]

Thanks for your response. Over all, do I have the right idea/logic? I’ll go ahead and use the uniform real number distribution.

[lastchance]

The logic for the Monte-Carlo method is OK (find the fraction of the points lying under the curve and then multiply by the outer rectangular area). The random-number generation is wrong: you are making it discrete when it should be continuous.

See the uniform real distribution example at
http://www.cplusplus.com/reference/random/uniform_real_distribution/

[jgialis]

I went ahead and implemented a continuous random number generator using the uniform real-number distribution. I'm still getting results way off however.

Thank you for all your help! I really appreciate it.

int main()
{
    int width = 25;
    int height = 5;
    default_random_engine generator;
    uniform_real_distribution<double> distributionX(0.0,width);
    uniform_real_distribution<double> distributionY(0.0,height);
    
    double x = 0.0;
    double y = 0.0;
    double iP = 0.0;

    int sampleSize = 10000;
    
    double functionY = 0;
 
    for (int i = 0; i < sampleSize; i++) {

        x = distributionX(generator);
        y = distributionY(generator);
        functionY = sqrt(x);
        
        cout << "(" << x << ", " << y << ")\n"; //Included this to see that I was indeed getting random values from 0 < x < 25, 0 < y < 5
        if(y <= functionY)
            iP++;

    }

  
    cout << "iP: " << iP << endl;
    cout << "Sample size: " << sampleSize << endl;
    cout << "Actual integral: " << (2.0/3 * pow(25, 1.5) - (2.0/3 * pow(0, 1.5)));
    cout << "\nMonteCarlo Integral Value: " << (width * height) * (iP)/(iP + sampleSize) << endl;
    
}

[Ganado]

When you post a code example, can you please also include any "#include" and "using" information before main?

Wouldn't the monte carlo integral value just be
(width * height) * (iP)/(sampleSize)
?

It's the ratio of "success points" (iP) to all points (sampleSize).

[lastchance]

Missed that! Yes, with @Ganado's correction it works perfectly.

[jgialis]

Ahh, how could I miss that!? I had it in my pseudocode too. Thanks for pointing it out! I appreciate all the help.

@Ganado I will post all the information next time. Sorry about that.

Have a great day.