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.

Thank you.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; }

Regards,

Rayne