Numerical instablity of ErfInverse with very small number

Dear ROOTers,

I calculate the statistical significance for cases where the p-value is very small,

the macro being:

// b = background, d = data
void ss( int b, int d ){

Double_t p = TMath::Poisson( d, b );
Double_t Z = TMath::Sqrt(2)*TMath::ErfInverse( 1 - p );
cout << p << endl;
cout << Z << endl;

return;
}

In my case the order of magnitude is b = 300, d = 600
, the corresponding p value if 8e-53, the ErfInverse functino returns 0 sigma for such value.
If I test the ErfInverse function for gradually smaller numbers, it is fine until e-17 when it suddenly drops to zero.

It looks like some some numerical unstability or a too rough precision in the estimation of the integral.
Do you confirm ? Is there a way to bypass this problem ?

Thanks in advance,

cheers,

Xavier

Hello,

these are common numerical instabilities that, as far I am concerned, can be faced adequately with one technique:

Rewriting your formula; playing with math.

For example look what ROOT outputs for these expressions:

[quote]1-TMath::Erf(10)
(const double)0.00000000000000000e+00

TMath::Erfc(10)
(Double_t)2.08848758376254457e-45[/quote]

that according to the definition of erfc

mathworld.wolfram.com/Erfc.html

they are equal.

In a similar fashion you can write erf^-1(1-x) like erfc^-1(x) according to eq. 8 of

mathworld.wolfram.com/InverseErf.html

Now:

root [36] TMath::ErfcInverse(1.e-85)
(Double_t)1.38749667080190164e+01

Check it with Mathematica. I think it works fine.

Cheers,
Leonidas

Hej Leonidas,

Thanks for the trick,

I am left for a particularly extreme case where b = 3000 and d = 6000,
which gives me:

root [12] ss( 3000, 6000, 1 )
Error in TMath::NormQuantile: probability outside (0, 1)
p = 0
Z = -0

Is there a way to go that far in the tails ? Otherwise I may simply quote a lower bound on the significance,

cheers,

Xavier

Hi,

for getting the significance, you should not use this function but

ROOT::Math::normal_quantile_c(p, 1)

which is the inverse of the integral of the normal distribution from +inf to x.
In your case by doing 1.-p if p is less than the double eps ( ~ 2 E-16 ) the result will be always 1.

Lorenzo