Re: tolerance in binopdf.m

[Dr. Alexander Klein]

Am 21.09.2011 um 07:23 schrieb Rik:

> On 09/20/2011 08:50 PM, Ben Abbott wrote:
> > On Sep 20, 2011, at 9:03 PM, Ben Abbott wrote:
> > More info on this particular test. After some simplification ...
> >
> > exp (gammaln (3)) * exp (2 * log (0.5)) - 0.5
> > ans = 0
> > exp (gammaln (3) + 2 * log (0.5)) - 0.5
> > ans =  1.1102e-16

Good morning,

I just tested this with 3.2.3 on OSX 10.5.8 which commes with lgamma  
in math.h, and in this case, too, the result of lgamma is off by 0.5ulp.

I did some further investigation, first on the Mac:

##Set up some factorials, Y doesn't have enough digits to warrant any  
rounding errors:
octave:37> X=1:15;Y=[1,cumprod(X)(1:end-1)];

##First let's see how well gamma fares
octave:38> Y-gamma(X)
ans =

 Columns 1 through 12:

   0.0000e+00   0.0000e+00   0.0000e+00   0.0000e+00   0.0000e+00    
0.0000e+00   0.0000e+00   0.0000e+00   0.0000e+00   0.0000e+00    
0.0000e+00  -2.9802e-08

 Columns 13 through 15:

   2.3842e-06   1.1444e-05   1.5259e-04
##That seems to be ok, even with the last few values the relative  
error is always smaller than 1e-14.


##Now let's see how lgamma behaves (assuming that log does)
octave:39> log(Y)-lgamma(X)
ans =

 Columns 1 through 12:

   0.0000e+00   0.0000e+00  -1.1102e-16  -2.2204e-16   4.4409e-16    
8.8818e-16   0.0000e+00   1.7764e-15   3.5527e-15   0.0000e+00    
0.0000e+00   0.0000e+00

 Columns 13 through 15:

  -3.5527e-15   3.5527e-15   0.0000e+00
##No real issues here, either.


##Just to be sure compare log(gamma()) to lgamma()
octave:40> log(gamma(X))-lgamma(X)
ans =

 Columns 1 through 12:

   0.0000e+00   0.0000e+00  -1.1102e-16  -2.2204e-16   4.4409e-16    
8.8818e-16   0.0000e+00   1.7764e-15   3.5527e-15   0.0000e+00    
0.0000e+00   0.0000e+00

 Columns 13 through 15:

  -7.1054e-15   0.0000e+00   0.0000e+00
##One value becomes better, one becomes worse, all others remain the  
same.



Now Xubuntu 10.04, also with 3.2.3 installed via Synaptic:

##Same setup as above:
octave:10> X=1:15;Y=[1,cumprod(X)(1:end-1)];


##Again, see how well gamma behaves
octave:11> Y-gamma(X)
ans =

 Columns 1 through 9:

   0.0000e+00   0.0000e+00   0.0000e+00   0.0000e+00  -3.5527e-15    
2.8422e-14  -1.1369e-13  -1.8190e-12  -2.1828e-11

 Columns 10 through 15:

   1.7462e-10  -2.3283e-09  -2.9802e-08   7.1526e-07  -1.0490e-05    
1.5259e-04
##That's somewhat worse than the Mac version.


##But now for the real surprises:
octave:12> log(Y)-lgamma(X)
ans =

   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0

octave:13> log(gamma(X))-lgamma(X)
ans =

   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0

Even with gamma having more roundoff error than on the mac, we get all  
zeroes here. Very interesting!

After all, I think that none of the above discrepancies warrant any  
concern. They're all in a range that is to be expected and may even be  
caused my simply moving data in and out of x87 FPU registers that  
provide some extra hidden bits of precision unless you force strict  
IEEE754 compliance for an implementation.


I think Michaels issues are much more worthwile considering.

Kind regards,

        Alex

P.s.: I was just asked to review a paper about roundoff error analysis  
where the algorithms presented differed by several orders of magnitude  
in their errors, and still none of them deserves to be labelled  
"bad". ;)

--
          Dr. Alexander Klein, Diplom-Mathematiker

Physiologisches Institut       |               TransMIT Zentrum
Raum 543                       |        für Numerische Methoden
Aulweg 129                     |          Heinrich-Buff-Ring 44
35392 Giessen                  |                  35392 Giessen