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