Accord.Net binomial probability mass function result differs from Excel result

Question

I'm trying to implement some functionality that was formerly provided via an Excel sheet into a C# application, but the probability mass function of Accord.NET differs for some reason from the excel function.

In excel the probabilty mass function, is used this way

=BINOM.DIST(250; 3779; 0.0638; FALSE)

Result: 0.021944019794458

When I try it with Accord.NET

var binom = new BinomialDistribution(3779, 0.0638);
binom.ProbabilityMassFunction(250);

// Result: Infinity

But the cumulative distribution seems to work properly (except for the last few digits, but I assumed this is just some kind of precision error)

Excel:

=BINOM.DIST(250; 3779; 0.0638; TRUE)

Result: 0.736156366002849

Accord.NET:

var binom = new BinomialDistribution(3779, 0.0638);
binom.DistributionFunction(250);

// Result: 0.736156366002318

Why are the results so different? And is there a way to get the Excel result with Accord?

EDIT: Extreme.Numerics calculates the same result as Excel, but I don't want to use this library, as the license system of this library always led to trouble in the past.

EDIT 2: Seems like some kind of overflow error.

When I use this I get the right result:

Math.Exp(binom.LogProbabilityMassFunction(250));

Any ideas why this could be happening?

Answer 1

Inspecting the source code of Accord.NET implementation it can be seen that ProbabilityMassFunction works like this:

//...
return Special.Binomial(numberOfTrials, k) * Math.Pow(probability, k)
            * Math.Pow(1 - probability, numberOfTrials - k);

In your case numberOfTrials is 3779 and k is 250. Binomial function of 3779 over 250 is more than 10³⁹⁸, which seems too much for the double data type (approximate range is up to a bit more than 10³⁰⁸). So you get Infinity in the first multiplication term, which determines the final result here.

I'm not sure how Extreme.Numerics implements that, but your results suggest that they do it in a more clever way. Also, sadly it doesn't seem that Accord.NET provides decimal version of their interface.

---

This issue seems to have been fixed in the latest Accord.NET implementation. The implementation looks now as follows:

double log = Special.LogBinomial(numberOfTrials, k) + k * Math.Log(probability)
    + (numberOfTrials - k) * Math.Log(1 - probability);

return Math.Exp(log);