CODEWITHSUNDEEP
c++random
lucacerone
lucacerone

Reputation: 10149

Select events having different probabilities - C++

I have to run a simulation where each event E_i has a certain probability P_i to happen (i=1..4, sum_i P_i = 1).

I would like to write a function that, taken as input the four different probabilities (they depend on other things happening in my simulation), returns the integer i for the chosen event.

Which is the best way to write such a function (possibly only using standard libraries, I have to run it on a cluster where I don't own the permissions to add new libraries)?.

I was thinking to write something like:

int get_event(double p1,double p2, double p3, double p4){
    double r=((double) rand())/(RAND_MAX); // alternatively here there can be any function that
    // generates a random number uniformly distribute in (0,1);
    if (r<=p1) return 1;
    else if (r>p1 && r <=(p1+p2)) return 2;
    else if (r>(p1+p2) && r <=(p1+p2+p3)) return 3;
    else if (r>(p1+p2+p3) && r <=1) return 4;
    else return -1; // -1 is an en error code;
}

but I am not sure if this is the best approach. Any suggestion? Thanks a lot in advance.

Upvotes: 0

Views: 2224

Answers (3)

Beta
Beta

Reputation: 99154

That's pretty good, I'd suggest only a few small changes:

int get_event(double p1,double p2, double p3, double p4){
    double r=((double) rand())/(RAND_MAX);
    if (r<=p1) return 1;
    r-= p1;
    if (r<=p2) return 2;
    r-= p2;
    if (r<=p3) return 3;
    r-= p3;
    if (r<=p4) return 4;
    return -1; // -1 is an en error code;
}

Upvotes: 3

Ivaylo Strandjev
Ivaylo Strandjev

Reputation: 71009

Here is an example how you can generalize your code while making it work for larger number of probabilities. Also I have included eps to the comparison. 

const double eps = 1e-9;
int get_event(double p[], int psize){
    double r=rand()/(RAND_MAX); // alternatively here there can be any function that
    // generates a random number uniformly distribute in (0,1);
    for (int i = 1; i <= psize; ++i) {
      r -= p[i-1];
      if (r < eps) {
        return i;
      }
    }
    return -1; // -1 is an en error code;
}

Upvotes: 2

Darcy Rayner
Darcy Rayner

Reputation: 3445

There is a method that works in O(1) time.

http://www.keithschwarz.com/darts-dice-coins/

http://www.altdevblogaday.com/2012/04/04/simulating-a-loaded-dice-in-a-constant-time

Upvotes: 3

Related Questions