Thread: 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

what is srand48() and drand48()? Are they C99 functions?

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

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