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!

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.