# Thread: Simulation - M(RT)^2 method

1. ## 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?

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. what is srand48() and drand48()? Are they C99 functions? 3. 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. ## 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);`
? 5. 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 