Simulation - M(RT)^2 method

This is a discussion on Simulation - M(RT)^2 method within the C Programming forums, part of the General Programming Boards category; Hi all, I'm supposed to generate random variates using the M(RT)^2 method, for p(x) ~ pi - x + x^2 ...

  1. #1
    Registered User
    Join Date
    Mar 2005
    Posts
    8

    Simulation - M(RT)^2 method

    Hi all,
    I'm supposed to generate random variates using the M(RT)^2 method, for p(x) ~ pi - x + x^2 + sin(x), for 0 <= x <= pi. I normalised p(x) to get p(x) = [6/(12 + 3pi^2 + 2pi^3)]*(pi - x + x^2 + sin(x)).

    I wrote the code below, to find out which value my parameter "ar" should be in order for the acceptance ratio to be 2/3. However, I kept getting an acceptance ratio of as high as 99.5%, regardless of my "ar" value. Also, I noticed that the random variates x_old and x_new just kept on either increasing or decreasing till they are very far off the range of [0, pi]. My probabilities of x_old and x_new (p_xold and p_xnew respectively) also went into the hundreds range the more points I tried to take. This is obviously wrong, isn't it?

    What is wrong with my code?

    Plesae advise.

    Code:
    #include <stdio.h>
    #include <stdlib.h>
    #include <time.h>
    #include <math.h>
    #define pi (3.141592654)
    
    int main(void)
    {
       int accepted = 0, fails = 0, i, N = 200000;
       float y = 0.0, rate = 0.0, x_old = 0.5, x_new = 0.0, p_xnew = 0.0, p_xold = 0.0, ar = 0.3;
       srand48(time(NULL));
       for (i = 0 ; i < N ; i++)
       {
            x_new =  (x_old + ar * (drand48() - 0.5));
            p_xnew = ((6.0)/(12.0 + 3.0*pi*pi + 2.0*pi*pi*pi))*(pi - x_new + x_new*x_new + sin(x_new));
            p_xold = ((6.0)/(12.0 + 3.0*pi*pi + 2.0*pi*pi*pi))*(pi - x_old + x_old*x_old + sin(x_old));
            if (p_xnew >= p_xold)
            {
          	accepted++;
                    x_old = x_new;
            }
            else
            {
          	fails++;
                    y = drand48();
                    if (y < (p_xnew/p_xold))
                   {
             	   accepted++;
                        fails--;
                        x_old = x_new;
                   }
            }
       }
       rate = accepted / (float)N;
       printf("Rate = %f\n", rate);
       return 0;
    }
    Thank you.

    Regards,
    Rayne

  2. #2
    ---
    Join Date
    May 2004
    Posts
    1,379
    what is srand48() and drand48()? Are they C99 functions?

  3. #3
    Registered User
    Join Date
    Mar 2005
    Posts
    8
    srand48() seeds the random-number generator.
    drand48() returns a random number in the range of 0 to 1. However, I don't think drand48() is usable except in a UNIX environment, which I'm using to run and compile my programs. It should be equivalent to rand()/RAND_MAX, I think.

    Regards,
    Rayne

  4. #4
    Registered User
    Join Date
    Mar 2005
    Posts
    36

    None.

    A quick gargle didn't show what your M(RT)^2 thing is, but looking at the rest, your "normalization" constant is about .06, but your expression can be very large if once it gets too big. Without much more, it looks like there is nothing preventing
    a string of high numbers from rand48 from forcing x_new up to where it can get out of control. Thus, did you mean
    Code:
    x_new = x_old(1 + ar * drand48);
    ?
    Last edited by Karthur; 03-13-2005 at 11:43 PM. Reason: invalid content

  5. #5
    .
    Join Date
    Nov 2003
    Posts
    307
    Two suggestions:

    M_PI is defined in math.h as pi.

    Use the double datatype instead of float, unless the specification requires float. reason: floats are accurate to 6 signficiant digits minimum, doubles to 10, minimum in ANSI C.

    limits.h has two defines
    FLT_DIG -- number of digits of precision in your float
    DBL_DIG -- number of digits of precision in your double. In 32 bit Unix this is frequently closer to 15 than 10

Popular pages Recent additions subscribe to a feed

Similar Threads

  1. on method pointers and inheritance
    By BrownB in forum C++ Programming
    Replies: 2
    Last Post: 03-02-2009, 06:50 PM
  2. fork, execv, spawn, or something else?
    By DavidP in forum C++ Programming
    Replies: 8
    Last Post: 01-26-2009, 03:25 PM
  3. Best communication method to thousand childs?
    By Ironic in forum C Programming
    Replies: 8
    Last Post: 11-07-2008, 11:30 PM
  4. Replies: 2
    Last Post: 01-22-2008, 03:22 PM
  5. Simulation - Composition method
    By wu_weidong in forum Tech Board
    Replies: 1
    Last Post: 03-12-2005, 10:17 AM

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