mathnet / mathnet-numerics

Math.NET Numerics
http://numerics.mathdotnet.com
MIT License
3.44k stars 891 forks source link

MarcumQ Function results are not equal. #1019

Closed WuFan closed 1 year ago

WuFan commented 1 year ago

matlab: marcumq(3.1622766,1.7941,1) = 0.9432

MathNet:
SpecialFuntions.MarcumQ(1,3.1622766,1.7941) = 0.80467

The results are not equal.

WuFan commented 1 year ago

different formula

diluculo commented 1 year ago

I'm not sure about the accuracy of the implemented SpecialFunctions.MarcumQ function. If the error is significant, we can easily calculate the Marcum Q function directly using the Integrate and BesselI functions, as shown in the code below, based on the explanation from Wikipedia.

var m = 1.0;
var a = 3.1622766;
var b = 1.7941;

var expected = 0.9432355485509051327956531; // from WolframAlpha: N[MarcumQ[1, 3.1622766,1.7941], 25]
var actual = 1 - 1.0 / Math.Pow(a, m - 1) * Integrate.GaussKronrod(x => Math.Pow(x, m) * Math.Exp(-(x * x + a * a) / 2) * SpecialFunctions.BesselI(m - 1, a * x), 0, b); 
//    actual = 0.94323554855090519
//    error = |actual - expected| = 1.1102230246251565E-16
vahokie02 commented 10 months ago

FWIW, there is significant error in SpecialFunctions.MarcumQ (well outside the stated target accuracy of 1.0e-11). While it may be possible to compute the value directly per @diluculo's suggested technique, shouldn't the method return the expected value? It almost seems like SpecialFunctions.MarcumQ is implementing some other function. I think this ticket should be reopened. I've attached an example plot to show what I'm talking about.

MarcumQ