March 13th, 2005, 05:01 AM
|
|
Registered User
|
|
Join Date: Feb 2005
Posts: 9
Time spent in forums: 2 m 37 sec
Reputation Power: 0
|
|
|
M(RT)^2 Algorithm
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
|
Dev Articles Community Forums Sponsor:
Del.icio.us
Digg
Google
Spurl
Blink
Furl
Simpy
Y! MyWeb
Linear Mode
