0

I am trying to calculate the upper cumulative probability of the binomial distribution in java. Since the probability values can go outside the double range, I am using BigDecimal. An example for total number of trials is n = 403 and number of successes k = 370. probability of success for each trial is (21/231). When I calculate the cumulative probability and subtract from 1 I get 1.1185E-14. However when I calculate the upper cumulative probability directly (from k=371 to 403) I get 5.2854E-341. I'm not able to understand why and would be grateful if someone could shed light on the discrepancy.

MathContext mtCtx = new MathContext(400, RoundingMode.HALF_UP);
BigDecimal frac1 = new BigDecimal(21.0/231.0, mtCtx);
BigDecimal frac2 = new BigDecimal(210.0/231.0, mtCtx);
Double cB = 0.0;
BigDecimal cBD = new BigDecimal(cB, mtCtx);
                
for (int num=0;num<=370;num++){
  BigDecimal pow1 = frac1.pow(num, mtCtx);
  BigDecimal pow2 = frac2.pow(403-num, mtCtx);
  
  BigInteger nminusk = bigFactorial(403-num);
  BigDecimal nminuskD = new BigDecimal(nminusk, mtCtx);
  BigInteger kFact = bigFactorial(num);
  BigDecimal kFactD = new BigDecimal(kFact, mtCtx);
  BigInteger nFact = bigFactorial(403);
  BigDecimal nFactD = new BigDecimal(nFact, mtCtx);

  BigDecimal binCoeff = (nFactD.divide((kFactD).multiply(nminuskD, mtCtx), mtCtx));
  
  
  BigDecimal prod1 = binCoeff.multiply(pow1);
  BigDecimal prod2 = prod1.multiply(pow2);
           
  cBD = cBD.add(prod2);
  
}
BigDecimal one = BigDecimal.ONE;
BigDecimal pVal = one.subtract(cBD);
System.out.println(pVal);

Thank you

Kaleido
  • 3
  • 2
  • Have edited my post to include it. – Kaleido Jun 17 '15 at 07:55
  • Waht is the return type of `bigFactoral`? I assume it is `BigInteger`. The division of a `BigInteger` returns again a `BigInteger`. Not sure but you might lose precision there. In addition: there is no need to transform the `BigInteger` to String in order to get a `BigDecimal` : `binCoeffD = new BigDecimal(binCoeff)` should be sufficent. – thomas Jun 17 '15 at 08:28
  • Yes it is BigInteger. How can there be a concept of precision if it's an integer? – Kaleido Jun 17 '15 at 08:35
  • I get what you mean about the precision now, and changed the code such that the division returns a BigDecimal. However, I still get the same discrepancy. – Kaleido Jun 17 '15 at 08:51
  • My idea was that the result of the division might be a real number and then BigInteger is not sufficent any more but I was wrong because the result is natural... But I see another problem: your definition of `frac1` and `frac2`. Before creating the BigDecimal, you do a double-division. Try `frac1 = new BigDecimal(21, mtCtx).divide(new BigDecimal(231, mtCtx), mtCtx)` – thomas Jun 17 '15 at 08:53
  • Yes, that was the problem! It works now – Kaleido Jun 17 '15 at 09:05
  • ok, posted it as an answer. Please accept (and eventually upvote ;) )the answer if it helped you out. – thomas Jun 17 '15 at 09:21

1 Answers1

0

As seen in the comments: there is a problem with the calculation of frac1 and frac2. Before creating a BigDecimal you are doing a double-division.

Try:

frac1 = new BigDecimal(21, mtCtx).divide(new BigDecimal(231, mtCtx), mtCtx);
frac2 = new BigDecimal(210, mtCtx).divide(new BigDecimal(231, mtCtx), mtCtx);
thomas
  • 5,637
  • 2
  • 24
  • 35