Monte Carlo methods

Hey, I would like to know how to implement Monte Carlo three methods in c++? Also where could I read about three different Monte Carlo sampling methods? Currently, I have written a pretty bad code, so I don't want to post anything until I write something that makes a bit of sense. Unfortunately, I don't understand Monte Carlo methods used in integration as well as I would. Need to read more. Anyhow, the problem I have is this:

Evaluate integral: I = Int(x(1 + 0.3x(1-x)))dx between 1 and 0, and find its divergence from 0.5.

Any help is very welcome!
Last edited on
The idea is that an integral is an area under a curve. That area can be estimated as a fraction of a known rectangular area with the same base by generating a random set of points in that rectangle and finding the fraction of them that lie below the curve.

Suppose you used a random-number generator to generate a set of points (x,y) with x uniformly distributed between 0 and 1 and y uniformly distributed between 0 and 2 (say - anything that spans the range of your integrand 1 + 0.3x(1-x) between those values of x will do). Find the fraction of points that lie below your integrand (f) and multiply it by the area of the whole rectangle (2 in this instance). If you have a large number of points, that will be the Monte-Carlo randomly-approximated area under your curve ... and hence an approximation to your integral.

It's similar in principle to a well-known Monte-Carlo method of approximating pi from the area of a quadrant compared with the area of a square.


BTW, this doesn't seem like"lounge" material - you could move the topic to "General C++ programming".




Thank you lastchance.
Last edited on
I need to find the integral - 0.5 for function x(1 + 0.3x(1-x)). So far I have successfully managed to find integral using first Monte Carlo method. Now I have two Monte Carlo methods left. I have written code but instead of getting the approximate integral value, I get -0.5 (II and III column in the result below). I believe there might be a problem somewhere for Monte_Carlo_2 and Monte_Carlo_3 functions. Basically, I would like to know what I need to change to get approximate values for integral - 0.5 for second and third option. The code and result are shown below:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
#include <iostream>
#include <stdio.h>
#include <cstdlib>
#include <cmath>

double rnd();

double function (double x);

void Monte_Carlo_1 (long* s);
void Monte_Carlo_2 (long* d);
void Monte_Carlo_3 (double x, long* f);

int main()
{
      double x;
      long max = 5, n = 0, s = 0, d = 0, f = 0;
      printf("\n*** Integralas ( I - 0.5 ): Monte Karlo metodas ***\n");
      printf("\nTikslus rezultatas = 0.0250 \n \n");
      printf("Iter. \t I \t II \t III \n \n");

      for (int j=0; j < 13; j++)
      {
              for (int i = 0; i < max; i++)
              {
                Monte_Carlo_1(&s);
                Monte_Carlo_2(&d);
                Monte_Carlo_3(x,&f);
                n++;
              }
              printf("%6ld %7.4lf %7.4lf %7.4lf \n", max, (double)s / n - 0.5, (double)d / n - 0.5, f / n - 0.5);
              max *= 2;
      }
}

double rnd()
{
      return (double)rand() / ( RAND_MAX );
}

void Monte_Carlo_1 (long* s)
{
      double x = rnd();
      double y = rnd();
      if(y < function(x))
      {
              (*s)++;
      }
}

void Monte_Carlo_2 (long* d)
{
      (*d) += rnd();
}

void Monte_Carlo_3 (double x, long* f)
{
      x = sqrt(rnd());
      *f += (x / (1 + 0.3 * x * (1 - x)) / x);
}

double function (double x)
{
      return (x + 0.3 * x * x - 0.3 * x * x * x);
}


*** Integral ( I - 0.5 ): Monte Carlo method ***

Exact  result = 0.0250 
 
Iter. 	 I 	 II 	 III 
 
     5  0.1000 -0.5000 -0.5000 
    10 -0.1000 -0.5000 -0.5000 
    20 -0.0429 -0.5000 -0.5000 
    40  0.0067 -0.5000 -0.5000 
    80 -0.0419 -0.5000 -0.5000 
   160 -0.0111 -0.5000 -0.5000 
   320  0.0402 -0.5000 -0.5000 
   640  0.0161 -0.5000 -0.5000 
  1280  0.0104 -0.5000 -0.5000 
  2560  0.0130 -0.5000 -0.5000 
  5120  0.0226 -0.5000 -0.5000 
 10240  0.0224 -0.5000 -0.5000 
 20480  0.0232 -0.5000 -0.5000  
Last edited on
Two points:
(1) Monte_Carlo_1 is the method that I would expect; I can't see what you are actually trying to do (mathematically) with your other two methods; please explain.

EDIT: one other method is not to randomise x - just take equally-spaced values; it is whether y lies below or above the function that matters.

(2) Why don't you just return the counts as the value of the functions? You also have no need to use pointers. Something like (untested):
1
2
3
4
5
6
7
8
long Monte_Carlo_1 ()
{
      long s = 0;
      double x = rnd();
      double y = rnd();
      if (y < function(x) ) s++;
      return s;
}

Actually, given the number of simulations you will need I should use unsigned long long rather than long.


Other things:
- you are presuming that your integrand lies between 0 and 1 when you set y - it happens to be true, but you can't guarantee it;
- there are better random number generators, including one which will give you a continuous uniform distribution (which is what you want).
Last edited on
From the book I read apparently there are three Monte Carlo sampling methods. First is calculating the area behind the curve where you divide area in small squares of sizes i and j. The second method (depicted Monte_Carlo_2) to get more accurate results is where you divide the area behind the curve into rectangles of size i. The third method involves non homogenous distributed numbers ONLY under the curve whereas we randomly generate numbers ALSO above the curve in the first and second methods, which isn't the case in the third MC method. In third MC we generate numbers according to function (p(x)) = 2x) in my case. However, I solved the problem anyhow. After some hard work I managed to solve it.

I could return, but there is little difference whether I use void or long function and return a value. Why would you suggest to return a value?
little difference whether I use void or long function and return a value. Why would you suggest to return a value?


It doesn't really matter, I guess; however, if a function produces a single output it feels slightly more natural (to me) to use that as the return value of the function than doing so through a function parameter ...
... which sort of leads to another suggestion: if you do want to return a value as a function parameter then make that a reference; then you won't need the unnecessary pointers.
Topic archived. No new replies allowed.