Why are the calculated values different for the binomial upper cumulative probability?

105 views Asked by At

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

1

There are 1 answers

0
thomas On BEST ANSWER

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);