0

I have some issues with C++ implementation of a function that generates random integers with fixed distribution. To do this I have implemented a function urand() that generates a double uniformly distributed between 0 and 1 (included) used by the function random_integer() [which uses a standard algorithm, I think it is described in one of Knuth's bool]. The functions are implemented like this:

double urand(){
int r =rand(); ;
return ((double) r)/RAND_MAX ;
} ; // uniform random number between 0 and 1 (included);


int random_integer(std::vector<double> &p_vec) {
unsigned int n_events=p_vec.size();
double u;
double initu;
do{
    u=urand();
}while(u<=0 && u >=1);
// chose a random number uniformly distributed between 0 and 1 (excluded)
initu=u;
for (unsigned int i=0; i<n_events; i++){
    if (u<=p_vec[i]) return i;
    u-=p_vec[i];
}
cout << "Warning: get_random_event() [utilities.cpp] could not find a suitable event. Returning -1 instead.\n";
cout << "p_vec=[" ; print(p_vec); cout << "]; "; // print is just a function that displays the elements in a vector;
cout << "initu=" << initu << "; u=" << u << endl;
return -1 ;
 }

Now most of the time everything works fine, but for example I got a warning like this:

Warning: get_random_event() [utilities.cpp] could not find a suitable event. Returning -1 instead. p_vec=[0.08;0.42;0.42;0.08[; initu=1; u=2.77556e-17

Now there are two things I don't understand: initu should be strictly less than 1 (look at the condition for the while loop). But even assuming that it i

1 Answer 1

4

I think you meant

}while(u<=0 || u >=1);

not

}while(u<=0 && u >=1);

u can't ever be both less or equal to than zero and greater than or equal to one.

Sign up to request clarification or add additional context in comments.

5 Comments

You are correct of course! What about the rounding errors in that code? Would they affect the result as well? Should I use a condition like while(u<=1e-12 || u>=1-1e12) ?
Rounding will not be an issue, 0 and 1 are represented exactly.
.. but the elements in p_vec not.. and in fact in the example I posted u=1, but at the end is 2e-17 rather than 0 as expected...
So make p_vec one bigger and put a zero there. Treat a return value of n_events specially.
That's great advice Doug! Though I think I should put a 1 there.. (if I put 0 even in the case I posted it won't work, since u = 2e-17 > 0; setting it to 1 I cover all the "un-catched" cases.. Thanks a lot for the help!

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.