Like Tree1Likes

2-d monte carlo integration help

This is a discussion on 2-d monte carlo integration help within the C Programming forums, part of the General Programming Boards category; Originally Posted by Adak rand() is not a great random number generator, but it's MORE than what you need for ...

  1. #16
    Registered User
    Join Date
    Oct 2011
    Posts
    854
    Quote Originally Posted by Adak View Post
    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)2 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)2 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.

  2. #17
    Registered User
    Join Date
    Sep 2012
    Posts
    3
    can someone show how to write the code for random, I don't understand clearly

  3. #18
    Registered User
    Join Date
    Sep 2006
    Posts
    8,868
    Quote Originally Posted by shootingstars View Post
    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).
    Last edited by Adak; 09-30-2012 at 06:26 PM.

  4. #19
    Registered User
    Join Date
    Sep 2012
    Posts
    48
    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.

  5. #20
    Registered User
    Join Date
    Sep 2012
    Posts
    3
    how do you work on the code so far?I think we study same class

  6. #21
    Registered User
    Join Date
    Sep 2012
    Posts
    48
    I've made a lot of progress, almost finished actually.

  7. #22
    Registered User
    Join Date
    Sep 2012
    Posts
    3
    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

  8. #23
    Registered User
    Join Date
    Sep 2012
    Posts
    48
    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.

  9. #24
    Registered User
    Join Date
    Sep 2012
    Posts
    48
    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...

  10. #25
    Registered User
    Join Date
    Sep 2012
    Posts
    48
    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.

  11. #26
    TEIAM - problem solved
    Join Date
    Apr 2012
    Location
    Melbourne Australia
    Posts
    1,497
    Code:
    while (i<1)
         i++;
         rn = ((double)rand()/(double)RAND_MAX)*(w-z);
    
    
    //You are missing '{' and '}' in your while statement
    Fact - Beethoven wrote his first symphony in C

  12. #27
    TEIAM - problem solved
    Join Date
    Apr 2012
    Location
    Melbourne Australia
    Posts
    1,497
    There is no return value
    Fact - Beethoven wrote his first symphony in C

  13. #28
    Registered User
    Join Date
    Sep 2012
    Posts
    48
    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()?
    Last edited by jsuite; 10-01-2012 at 02:54 AM.

  14. #29
    Registered User
    Join Date
    Sep 2006
    Posts
    8,868
    Quote Originally Posted by shootingstars View Post
    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.

  15. #30
    Registered User
    Join Date
    Sep 2006
    Posts
    8,868
    Quote Originally Posted by jsuite View Post
    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.
    Last edited by Adak; 10-01-2012 at 03:12 AM.

Page 2 of 5 FirstFirst 12345 LastLast
Popular pages Recent additions subscribe to a feed

Similar Threads

  1. Monte Carlo Simulation Optimization
    By Phys10 in forum C++ Programming
    Replies: 5
    Last Post: 08-02-2010, 04:09 AM
  2. Help with C Monte Carlo
    By shaunmikex in forum C Programming
    Replies: 7
    Last Post: 03-21-2010, 05:11 PM
  3. Monte Carlo reliability program - SegFault :(
    By Kibble Fat in forum C Programming
    Replies: 7
    Last Post: 12-15-2009, 01:41 PM
  4. Problem with Monte Carlo(integration by rejection)
    By dionys in forum C Programming
    Replies: 9
    Last Post: 04-07-2004, 09:53 AM
  5. monte carlo methods
    By qwertiop in forum C++ Programming
    Replies: 3
    Last Post: 09-05-2001, 11:29 PM

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21