Re: tolerance in binopdf.m

[Ben Abbott]

On Sep 20, 2011, at 9:03 PM, Ben Abbott wrote:

> On Sep 20, 2011, at 8:33 PM, Rik wrote:
> 
>> On 09/20/2011 04:42 PM, Michael D Godfrey wrote:
>>> On 09/20/2011 04:28 PM, Ben Abbott wrote:
>>>> Rik,
>>>> 
>>>> I'm happy to see you complete this effort.
>>>> 
>>>> With this change, I see 5 failures. They all appear to be due to 
>>>> tolerances. All tests pass if I use a tolerance of eps(). My diff is also 
>>>> below.
>>>> 
>>>> Is adding tolerance appropriate? If so, I'm happy to push a changeset.
>>>> 
>>>> Ben
>>>> 
>>> Running on Linux, I see no failures.  So, this may need some checking.
>>> 
>>> Michael
>>> 
>> Ben,
>> 
>> I'm not too concerned with the tpdf instance so I pushed a changeset to add 
>> an eps tolerance in that case.
>> 
>> Like Michael, I'm running Linux and I don't see any issue with binopdf.  My 
>> guess is that it is the function gammaln() which is returning ever so 
>> slightly different values.  It shouldn't be too hard to track down where the 
>> first difference arises in the calculation chain.  However, after that I'm 
>> not sure what we would do since I don't think you have ready access to a 
>> different version of gammaln().  Thus, I'm tempted to just add eps to the 
>> end of all of the tests.  It's a bit of a shame as the values in question, 
>> [1/4 1/2 1/4], are simple and exact.
>> 
>> I'll let Michael comment on whether there is a better approach.
>> 
>> --Rik
> 
> I'm using gcc 4.4 (non-Applified) with "-O2".
> 
> If there are any suggestions, I'm happy to try them. If the problem turns out 
> to be with one of Apple's lib, then I'm ok with something like ...
> 
>       if (ismac ())
>         tol = eps ();
>       else
>         tol = 0;
>       endif
> 
> For reference, the fist problematic test ...
> 
> x = [-1 0 1 2 3];
> y = [0 1/4 1/2 1/4 0];
> z = binopdf (x, 2 * ones (1, 5), 0.5 * ones (1, 5));
> z - y
> ans =   0.0000e+00   0.0000e+00   1.1102e-16   0.0000e+00   0.0000e+00
> 1/2 - 0.5
> ans = 0
> 
> The result from binopdf is off by eps/2.
> 
> Ben

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

I'm not sure the simplification below is immune to roundoff, but ...

gammaln (3) + log (0.5)
ans =  1.1102e-16

(fwiw, I get the same result with ML)

Looking at these individually ...

printf ("%.17f\n", log (0.5))
-0.69314718055994529
 printf ("%.17f\n", gammaln (3))
0.69314718055994540

Does anyone have any insight regarding this?

What does Linux give for these two?

Ben