Search code examples
c++combinationsprobabilityfactorial

Calculating Probability C++ Bernoulli Trials


The program asks the user for the number of times to flip a coin (n; the number of trials).

A success is considered a heads.

Flawlessly, the program creates a random number between 0 and 1. 0's are considered heads and success.

Then, the program is supposed to output the expected values of getting x amount of heads. For example if the coin was flipped 4 times, what are the following probabilities using the formula

nCk * p^k * (1-p)^(n-k)
Expected 0 heads with n flips: xxx
Expected 1 heads with n flips: xxx
...
Expected n heads with n flips: xxx

When doing this with "larger" numbers, the numbers come out to weird values. It happens if 15 or twenty are put into the input. I have been getting 0's and negative values for the value that should be xxx.

Debugging, I have noticed that the nCk has come out to be negative and not correct towards the upper values and beleive this is the issue. I use this formula for my combination:

double combo = fact(n)/fact(r)/fact(n-r);

here is the psuedocode for my fact function:

long fact(int x)
{
     int e; // local counter
     factor = 1;
     for (e = x; e != 0; e--)
     {
         factor = factor * e;
     }
     return factor;
}

Any thoughts? My guess is my factorial or combo functions are exceeding the max values or something.


Solution

  • You haven't mentioned how is factor declared. I think you are getting integer overflows. I suggest you use double. That is because since you are calculating expected values and probabilities, you shouldn't be concerned much about precision.

    Try changing your fact function to.

    double fact(double x)
    {
         int e; // local counter
         double factor = 1;
         for (e = x; e != 0; e--)
         {
             factor = factor * e;
         }
         return factor;
    }
    

    EDIT: Also to calculate nCk, you need not calculate factorials 3 times. You can simply calculate this value in the following way.

    if k > n/2, k = n-k.
    
           n(n-1)(n-2)...(n-k+1)
    nCk = -----------------------
              factorial(k)