2-d monte carlo integration help

**Nominal Animal** · 09-30-2012

Originally Posted by Adak

rand() is not a great random number generator, but it's MORE than what you need for this program.

Adak, I agree it is enough for this exact program. However, in the context of the title of this thread, "2-d monte carlo integration help", rand() is not good enough.

Traditionally RAND_MAX has been very small, typically 32767. Regardless of its period (the length of the sequence of numbers it generates), you can have at most (RAND_MAX+1)² unique random 2D dart positions. On a typical desktop machine, that takes just seconds, not minutes.

In Linux, RAND_MAX == 2147483647, but the generator is still weak; the dart positions form a quite regular lattice. I also think it generates much fewer than (RAND_MAX+1)² unique 2D points, but I haven't verified that.

While using srand() and rand() is okay for the purposes of an exercise, it is never okay to use it with Monte Carlo stuff, it is simply not random enough. The implementation may work, but it is does not use true random sampling; it'll use some implementation-dependent regular or semi-regular sampling scheme. And that's not Monte Carlo, is it?

Instead of random sampling, using a linear congruential pseudorandom number generator like rand() and random(), leads to a sampling pattern that approaches a uniform lattice as the number of samples increases. This is most evident if you draw the points into a picture, and watch it fill up.

Furthermore, using (double)rand()/(double)RAND_MAX will only give you RAND_MAX+1 unique values -- just 32768 on many non-Linux systems! The smallest nonzero value (double)rand() / ((double)RAND_MAX + 1.0) returns may be as big as 0.00003, and the largest 0.99997. In other words, less than five significant digits..

The Xorshift alternative I showed in this thread does not suffer from these issues at all, and the prng_one() yields all doubles within [0, 1), about eighteen significant digits per value.

To OP:

While I understand this is just for an exercise, it is very important to realize the impact a poor random number generator has on Monte Carlo and other random-number based computations. It is very common, and very easy to just overlook that, and assume the results are good. Most often the result itself is okay, but the error bars obtained from the sampling statistics may be completely bogus, wonkified by the non-random random number generator used.

In this particular exercise you don't need to use random sampling; a regular sampling will give you the best answer anyway. The reason why should be obvious, right? Because of that, Adak is absolutely correct: rand() will work just fine for this program.

I apologize for harping about this so much; I'll try and be quiet now.

**shootingstars** · 09-30-2012

can someone show how to write the code for random, I don't understand clearly

**Adak** · 09-30-2012

Originally Posted by shootingstars

can someone show how to write the code for random, I don't understand clearly

You said it, Pal! That's one of the next issues that needs to be nailed down. Stay tuned. Run the snippet program in reply #15, and run it with both of the options (one at a time), to see what they do.

@Nominal Animal:
Please, hop off the "random" soap box. It's all true, (and great info), but not needed here (yet).

**jsuite** · 09-30-2012

Wow I appreciate all the info on rand! Unfortunately as the instructions say, I must use rand.

Adak: Your rand function helped a lot, but one minor thing i need to make it work with:

Ill right some generic code here real quick:

randNumInRange(int w, int x) (where w is the max, and x is the min)
{

do stuff / arguments here

My point is (don't know the terminology) but it's a general function, so I can pass it in any other function as randNumInRange(insert max variable here, min variable here) So when I call on it, I can pass any min / max through it. Hope this makes since, you definitely have me on the right track, working on it now. Thanks!!

ps. Ill ask more about getting the function to work if I run into errors.

**shootingstars** · 09-30-2012

how do you work on the code so far?I think we study same class

**jsuite** · 09-30-2012

I've made a lot of progress, almost finished actually.

**shootingstars** · 09-30-2012

you are so great,I think I can not finish this one although I already try my best,I have no experience with C before

**jsuite** · 09-30-2012

Adak, I get errors when I try to compile your code. Here they are:

/usr/include/stdlib.h:140:1: error: unknown type name ‘size_t’
In file included from ran.c:2:0:
/usr/include/stdlib.h:337:4: error: unknown type name ‘size_t’
/usr/include/stdlib.h:367:4: error: unknown type name ‘size_t’
/usr/include/stdlib.h:471:22: error: unknown type name ‘size_t’
/usr/include/stdlib.h:473:22: error: unknown type name ‘size_t’
/usr/include/stdlib.h:473:38: error: unknown type name ‘size_t’
/usr/include/stdlib.h:485:36: error: unknown type name ‘size_t’
In file included from /usr/include/stdlib.h:497:0,
from ran.c:2:
/usr/include/alloca.h:33:22: error: unknown type name ‘size_t’
In file included from ran.c:2:0:
/usr/include/stdlib.h:503:22: error: unknown type name ‘size_t’
/usr/include/stdlib.h:508:45: error: unknown type name ‘size_t’
/usr/include/stdlib.h:508:65: error: unknown type name ‘size_t’
/usr/include/stdlib.h:756:9: error: unknown type name ‘size_t’
/usr/include/stdlib.h:756:25: error: unknown type name ‘size_t’
/usr/include/stdlib.h:761:34: error: unknown type name ‘size_t’
/usr/include/stdlib.h:761:50: error: unknown type name ‘size_t’
/usr/include/stdlib.h:840:6: error: unknown type name ‘size_t’
/usr/include/stdlib.h:843:6: error: unknown type name ‘size_t’
/usr/include/stdlib.h:847:31: error: unknown type name ‘size_t’
/usr/include/stdlib.h:851:31: error: unknown type name ‘size_t’
/usr/include/stdlib.h:860:38: error: unknown type name ‘size_t’
/usr/include/stdlib.h:864:36: error: unknown type name ‘size_t’
/usr/include/stdlib.h:871:1: error: unknown type name ‘size_t’
/usr/include/stdlib.h:872:34: error: unknown type name ‘size_t’
/usr/include/stdlib.h:874:1: error: unknown type name ‘size_t’
/usr/include/stdlib.h:875:40: error: unknown type name ‘size_t’
ran.c: In function ‘main’:
ran.c:18:7: warning: incompatible implicit declaration of built-in function ‘printf’ [enabled by default]
ran.c:20:4: warning: incompatible implicit declaration of built-in function ‘printf’ [enabled by default]

I think there's number truncation going on..

and stdlib not supported?

NEVER MIND #include ,stdio.h> got cut off when ipasted it over.

**jsuite** · 09-30-2012

Ok so the second snippet, the one between 0 and 1... That one is how they want it to be done BUT that one produces decimals. Such as .23, etc.. The first one behaves correctly, but isn't the the wat he wants us to do it. So, I need to make the second snippet produce doubles that are whole numbers (12.0, 40.0, 70,0) etc...

**jsuite** · 09-30-2012

ok so when I try this:

Code:

#include <stdio.h>
#include <stdlib.h>
double RandNumInRange (double, double);




int main (void)
{
double xmax, xmin;


xmax = 50;


xmin = 10;


double rand = RandNumInRange (xmin, xmax);


printf("\n\nThe Random Number In Range = %f \n\n\n", rand);
}


double RandNumInRange (double z, double w) //where w is any max value, and z is any min value
{
   double rn;
   int i = 0;
   srand(1);


   while (i<1)
     i++;
     rn = ((double)rand()/(double)RAND_MAX)*(w-z);




   }

I get the same number over and over again. Any ideas? D:

The number only changes when I change the seed, also it's not a whole number.

**Click_here** · 10-01-2012

Code:

while (i<1)
     i++;
     rn = ((double)rand()/(double)RAND_MAX)*(w-z);


//You are missing '{' and '}' in your while statement

**Click_here** · 10-01-2012

There is no return value

**jsuite** · 10-01-2012

Urgent question: This is needed for the loop in ThrowDarts. So Throw darts basically checks if X and Y are inside the rectangle, then at the end of the loop it takes the number in range of x, and number in range of y, and then checks if it's inside the circle. How do I know if those numbers are in the coordinates of the circle? Basically need help / explanation of DartHits function.

Part c of loop instructions:

(4) Loop over N darts thrown, and count hits
a. Generate random # between [xmin,xmax]; call this xdart (xcoord of dart)
-> call xdart = RandNumInRange(xmin, xmax);
b. Generate random # between [ymin,ymax]; call this ydart (ycoord of dart)
-> call ydart = RandNumInRange(ymin, ymax);
c. hit = DartHits(xdart, ydart);
e. if hit != 0, add one to # of hits

Found the answer:

TESTER:
(1) Write function for shape we know answer to
--> circle has area PI*r*r
--> Easy to write function for circle wt center on axis!
--> Answer 1 if distance of (x,y) from axis <= r
--> distance: sqrt(x*x + y*y) (pythagorean theorem):
** draw grid wt triangle, to show how hypotenuse is the distance
--> RETURN 1 if sqrt(x*x + y*y) <= r, ELSE 0
==> Call this function:
int DartHits(double x, double y)

So due to the mixed references to "tester" I'm guessing "tester" is more than one function, instead a couple of functions that come before ThrowDarts()?

**Adak** · 10-01-2012

Originally Posted by shootingstars

you are so great,I think I can not finish this one although I already try my best,I have no experience with C before

@shootingstars:

Post up your code, and best if you start it in a new thread so it gets the attention it deserves.

Once you get it started, we can help - just tell us what the problem/s is/are.

Do it ASAP, because it will take time.

**Adak** · 10-01-2012

Originally Posted by jsuite

Urgent question: This is needed for the loop in ThrowDarts. So Throw darts basically checks if X and Y are inside the rectangle, then at the end of the loop it takes the number in range of x, and number in range of y, and then checks if it's inside the circle. How do I know if those numbers are in the coordinates of the circle? Basically need help / explanation of DartHits function.

Part c of loop instructions:

(4) Loop over N darts thrown, and count hits
a. Generate random # between [xmin,xmax]; call this xdart (xcoord of dart)
-> call xdart = RandNumInRange(xmin, xmax);
b. Generate random # between [ymin,ymax]; call this ydart (ycoord of dart)
-> call ydart = RandNumInRange(ymin, ymax);
c. hit = DartHits(xdart, ydart);
e. if hit != 0, add one to # of hits

Found the answer:

TESTER:
(1) Write function for shape we know answer to
--> circle has area PI*r*r
--> Easy to write function for circle wt center on axis!
--> Answer 1 if distance of (x,y) from axis <= r
--> distance: sqrt(x*x + y*y) (pythagorean theorem):
** draw grid wt triangle, to show how hypotenuse is the distance
--> RETURN 1 if sqrt(x*x + y*y) <= r, ELSE 0
==> Call this function:
int DartHits(double x, double y)

So due to the mixed references to "tester" I'm guessing "tester" is more than one function, instead a couple of functions that come before ThrowDarts()?

Actually, I was wondering if "tester" was BOTH a function, AND sometimes referred to the program "testing" in the simulation. It is confusing.

This is just the Pythagorean theorem:
RETURN 1 if sqrt(x*x + y*y) <= r, ELSE 0

Where x = A, y = B, and r = C² in A² + B² = C²

So to tell if it's inside the circle, get the x and y point of the hit, and the x and y values of the center of the circle. Then:

We get a triangle from the hit, to the center of the circle:

X -x == base distance of the triangle
Y - y == height distance of the triangle

Now base² + height² == hypotenuse² (distance as the crow flies), from the hit, to the center of the circle.

If hypotenuse is greater than the radius of the circle, then it's outside the circle. If it's less than or equal to, it's inside the circle.

For instance, if the triangle had base of 3, and height of 4, then the hypotenuse would be 5, since:
3² + 4² = 5²

So if the radius of the circle was > 5, the hit would be outside the circle. If it's less than or equal to 5, then it's inside it.