Calculating Probability C++ Bernoulli Trials

Question

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.

Answer 1

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)

Answer 2

You're exceeding the maximum value of a long. Factorial grows so quickly that you need the right type of number--what type that is will depend on what values you need.

Long is an signed integer, and as soon as you pass 2^31, the value will become negative (it's using 2's complement math).

Using an unsigned long will buy you a little time (one more bit), but for factorial, it's probably not worth it. If your compiler supports long long, then try an "unsigned long long". That will (usually, depends on compiler and CPU) double the number of bits you're using.

You can also try switching to use double. The problem you'll face there is that you'll lose accuracy as the numbers increase. A double is a floating point number, so you'll have a fixed number of significant digits. If your end result is an approximation, this may work okay, but if you need exact values, it won't work.

If none of these solutions will work for you, you may need to resort to using an "infinite precision" math package, which you should be able to search for. You didn't say if you were using C or C++; this is going to be a lot more pleasant with C++ as it will provide a class that acts like a number and that would use standard arithmetic operators.