I have the following c function working on implementing a rejection sampler. MersenneTwiser.h generates a random number between 0 and 1. The issue I'm having is with the rnorm function. I'm using a while loop to reject some samples. It's nonsense at the moment as it's not finished. In it's current form, the program returns 0.00000.
#include<stdio.h>
#include<math.h>
#include"MersenneTwister.h"
MTRand rng;
double p;
double prop;
int sign;
double dubexp(double theta){
p = rng.rand();
if ((p-0.5) > 0) sign= 1;
if ((p-0.5) < 0) sign= -1;
prop = -(1/theta)*sign*log(1-2*abs(p-0.5));
return(prop);
}
double u;
double theta;
double t;
double z;
double c;
double x;
double rnorm(double theta, double c){
t=rng.rand();
while (z == t)
{
x=dubexp(theta);
u=rng.rand();
z=x;
}
return z;
}
int main(){
theta=1;
c=1;
u = rnorm(theta,c);
printf("%f",u);
}
However if I remove the while loop, it returns the correct value of z. As below:
double rnorm(double theta, double c){
t=rng.rand();
x=dubexp(theta);
u=rng.rand();
z=x;
return z;
}
Aucun commentaire:
Enregistrer un commentaire