Important numbers by Monte Carlo Simulation

The code below (should) produce some decent approximations to certain mathematical numbers by Monte Carlo simulation (i.e. using random numbers). They are:
PI
ln(2) - natural log of 2
e - base of natural logarithms

Purely for recreation and general enlightenment (because I don't know!), does anyone know any other key mathematical numbers that one can produce purely by Monte Carlo methods? I'm after, in particular: sqrt(2), the golden ratio (phi) and Euler-Mascheroni constant (gamma).

You need about 100000 random numbers to get a decent approximation for pi; 1 million plus is better.
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
66
67
68
69
70
71
72
73
74
75
76
77
78
79
#include <iostream>
#include <chrono>
#include <random>
using namespace std;

unsigned seed = chrono::system_clock::now().time_since_epoch().count();
default_random_engine generator(seed);
uniform_real_distribution<double> distribution( 0.0, 1.0 );


double myRandom()           // needs to return a uniformly-distributed value on [0,1]
{
   double x;
   x = distribution( generator );
   return x;
}


double getPi( long N )
{
   long ntrue = 0;          // number meeting criterion
   double x, y;             // random points

   for ( long i = 0; i < N; i++ )
   {
      x = myRandom();
      y = myRandom();
      if ( x * x + y * y < 1 ) ntrue++;
   }

   return 4.0 * ntrue / N;
}


double getLog2( long N )
{
   long ntrue = 0;          // number meeting criterion
   double x, y;             // random points

   for ( long i = 0; i < N; i++ )
   {
      x = myRandom() + 1;
      y = myRandom();
      if ( x * y < 1 ) ntrue++;
   }

   return ( ntrue + 0.5 ) / N;
}


double getE( long N )
{
   long count;              // count in each sequence
   double sum;              // sum of random numbers
   double total;            // running count of min to get average(min) later

   for ( long i = 0; i < N; i++ )
   {
      sum = 0;
      while ( sum < 1 )
      {
         sum += myRandom();
         total++;
      }
   }

   return ( total + 0.5 ) / N;
}


int main()
{
   long N;
   cout << "Input number of random points: ";   cin >> N;
   cout << "\nMonte-Carlo estimates for: ";
   cout << "\npi:    " << getPi( N );
   cout << "\nln(2): " << getLog2( N );
   cout << "\ne:     " << getE( N );
}

Last edited on
For sqrt(2):
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
#include <random>
#include <iostream>
#include <cmath>

int main()
{
    std::random_device rd;
    std::mt19937 gen(rd());
    std::uniform_real_distribution<> dis(0, 2);
    long long ntrue{};
    long long N = {};
    for (long long n = 0; n < N; ++n)
    {
        if(pow(dis(gen), 2) <= 2)
        {
            ntrue++;
        }
        n++;
    }
    float getSqrt2 = (float)ntrue*4/N;
    std::cout<<getSqrt2;
}


N getSqrt2
10^5 1.40816
10^6 1.4155
10^8 1.41397
10^9 1.41422
Ah, thanks GunnerFunner. I missed the (nearly) obvious there. I've added your method to the code below (but translated so that random numbers are uniform on [0,1] as in the rest of the functions).

Actually, all square roots work in the same way, so this indirectly gives the golden ratio ("divine proportion") as well.

Also added another method for pi (based on ratio of volumes of sphere and cube) rather than the classic one based on the ratio of areas (circle to square). Natural logarithms for any number can also be derived from an area ratio (under xy=1) in a similar manner to log(2).

If anyone knows of any others I would be glad to hear it.
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
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
#include <iostream>
#include <iomanip>
#include <chrono>
#include <random>
using namespace std;

unsigned seed = chrono::system_clock::now().time_since_epoch().count();
default_random_engine generator(seed);
uniform_real_distribution<double> distribution( 0.0, 1.0 );


double myRandom()           // needs to return a uniformly-distributed value on [0,1]
{
   double x;
   x = distribution( generator );
   return x;
}


double getPi( long N )      // from ratio of area of a circle to the unit square
{
   long ntrue = 0; 
   double x, y;    

   for ( long i = 0; i < N; i++ )
   {
      x = myRandom();
      y = myRandom();
      if ( x * x + y * y <= 1.0 ) ntrue++;
   }

   return 4.0 * ntrue / N;
}


double getPiByVolume( long N )         // from ratio of volume of a sphere to a unit cube
{
   long ntrue = 0; 
   double x, y, z; 

   for ( long i = 0; i < N; i++ )
   {
      x = myRandom();
      y = myRandom();
      z = myRandom();
      if ( x * x + y * y + z * z <= 1.0 ) ntrue++;
   }

   return 6.0 * ntrue / N;
}


double getLog2( long N )               // from area under curve y=1/x (i.e. xy=1) from x = 1 to 2
{
   long ntrue = 0;       
   double x, y;          

   for ( long i = 0; i < N; i++ )
   {
      x = myRandom() + 1;
      y = myRandom();
      if ( x * y <= 1.0 ) ntrue++;
   }

   return ( ntrue + 0.5 ) / N;
}


double getE( long N )                  // see Wikipedia: I don't get it though!
{
   double sum;
   long long total = 0;     

   for ( long i = 0; i < N; i++ )
   {
      sum = 0;
      while ( sum <= 1 )
      {
         sum += myRandom();
         total++;
      }
   }

   return ( total + 0.5 ) / N;
}


double getRoot2( long N )              // from pdf; sqrt(A) = A * Prob( X * X < 1/A )
{
   long ntrue = 0;         
   double x;                

   for ( long i = 0; i < N; i++ )
   {
      x = myRandom();
      if ( x * x <= 0.5 ) ntrue++;
   }

   return 2.0 * ntrue / N;
}


double getPhi( long N )                // from sqrt(5) as above
{
   long ntrue = 0;         
   double x;               

   for ( long i = 0; i < N; i++ )
   {
      x = myRandom();
      if ( x * x <= 0.2 ) ntrue++;
   }

   return 0.5 + 2.5 * ntrue / N;
}


void printIt( string target, double (*f)( long ), long N, string actual )
{
   cout << setw(11) << target << setw(12) << setprecision(5) << f(N) << "    ( " << actual << " )" << endl;
}


int main()
{
   long N;
   cout << "Input number of random points: ";   cin >> N;
   cout << "\nMonte-Carlo estimates for: " << endl;
   printIt( "pi:"       , getPi        , N, "3.14159..." );
   printIt( "pi(again):", getPiByVolume, N, "3.14159..." );
   printIt( "ln(2):"    , getLog2      , N, "0.693147.." );
   printIt( "e:"        , getE         , N, "2.71828..." );
   printIt( "sqrt(2):"  , getRoot2     , N, "1.41421..." );
   printIt( "phi:"      , getPhi       , N, "1.61803..." );
}
Last edited on
There's another way I was thinking of Monte-Carlo for irrational numbers: sqrt(5) is the hypotenuse of a right triangle of sides length 1 and 2. So if we dropped N needles randomly on the perimeter of the triangle and making the reasonable assumption that the needles fell on each side in proportion to it's length (and overlooking needles falling on the vertices for sufficiently large N) we could, in the limit, get the length of sqrt(5) noting how many needles fell on the other sides, whose lengths we know. A similar approach (Buffon's needle) is used to approximate pi via Monte Carlo but I have not had a chance to think about coding it. If anybody else finds this of interest and want to take it forward they are most welcome.
Years ago at school we tried the actual physical experiment for Buffon's needle to approximate PI. Unfortunately, half the class were trying desperately hard to hit the target and the other half were mischievously trying to avoid it, so the fall of our sticks was anything but random. It would be interesting to see what a computer simulation came up with.
Buffon's needle estimate of pi (https://en.wikipedia.org/wiki/pi )

Statement: drop a needle, length L, onto paper with parallel lines, distance t apart. Provided that L <= t ("short needle") the probability that the needle crosses a line is
2L/(pi t)
and hence pi can be estimated from
pi = 2(L/t)N/n
where N is the number of drops and n is the number of times it crosses a line.

OK, the code below simulates the physical experiment, but it hits one BIG SNAG - in order to randomly orient the needle it has to use pi ... the number it is trying to estimate!
Suggestions as to how to avoid this please!
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
#include <iostream>
#include <chrono>
#include <random>
using namespace std;

unsigned seed = chrono::system_clock::now().time_since_epoch().count();
default_random_engine generator(seed);
uniform_real_distribution<double> distribution( 0.0, 1.0 );


double myRandom() { return distribution( generator ); }    // returns a uniformly-distributed value on [0,1]


void needleEnds( double L, double &y1, double &y2 )        // simulates dropping a needle, length L. Returns the end y coordinates, y1 and y2
{
   const double pi = 3.14159265358979323846;               // embarassing - this is what we are attempting to find!!!

   y1 = myRandom();                         // One end of the needle ( randomly distributed in  [0, 1] )
   double theta = myRandom() * 2.0 * pi;    // Angle of the needle   ( randomly distributed in [0, 2pi] )
   y2 = y1 + L * sin( theta );              // Other end of the needle
}


int main()
{
   long N;
   long ntrue = 0; 
   double L;
   double y1, y2;

   do
   {
      cout << "Input relative length of needle L/t (NOTE: <= 1): ";
      cin >> L;
   } while ( L > 1.0 );

   cout << "Input number of needle drops: ";
   cin >> N;

   for ( long i = 0; i < N; i++ )
   {
      needleEnds( L, y1, y2 );              // drop the needle
      if ( y2 <= 0 || y2 > 1.0 ) ntrue++;   // does the needle cross a line?
   }

   cout << "\nEstimate of pi by Buffon's needle is " << 2.0 * L * N / ntrue;
}


Last edited on
I doubt you'd see much difference for most values you could use for theta. Larger values will tend to correct any bias introduced by not using an exact multiple of 2pi.
A way to estimate pi on the principle of Buffon's needle, though not quite the exact implementation, that does not require any prior pi assumptions is to drop the needles on to a unit square containing quadrant of a circle and then counting off x^2 + y^2 <= 1; where x, y ~ UID[0,1]

http://scicomp.stackexchange.com/questions/11031/monte-carlo-estimation-of-pi-with-buffons-needle
cire - That's a good point. Certainly if you allowed theta to range over something of the order of the number of throws then it would be difficult to argue that it created any more error than using a finite number of throws already did.

gunnerfunner - thank-you for sharing that link! I'm relieved to hear that I'm not the only one coming up against this problem. I think your circle quadrant / unit square approach is essentially that which I had earlier coded, based on the ratio of areas.
Topic archived. No new replies allowed.