Re: eps function and numerical accuracy

[Daniel J Sebald]

On 07/08/2013 07:21 PM, Nir Krakauer wrote:
> Yes, it looks like assert is giving the right result but the wrong
> error message.

I've revamped the assert.m script to allow non-scalar tolerance values. 
 The changeset is here:

https://savannah.gnu.org/patch/index.php?8119

Rik, could you please look through the changeset?  I think it has some 
improvements and the variable tolerance is quite nice.  Only the 
comparisons which fall outside of tolerance are displayed.  The error 
messages are formatted better and make sense in terms of what is out of 
tolerance and by how much.

For that reason I've removed the FIXME about formatting at the top of 
the script file.  Let's tweak the output format as desired and be done 
with that.  I think it is pretty close as changed.

Instead of keeping track of error state via variable "iserror" all the 
way to the end of the routine, it uses a subfunction for printing the 
error.  In addition to a little more flexibility of output error 
message, this has the added advantage of displaying the line number 
about where the tests fails.  With the former approach, the line number 
was always somewhere at the end of the file.

The varying tolerance comes in handy for testing function results for 
which there is a wide numerical range--when a large tolerance value (say 
eps) doesn't really test most of the routine and a small tolerance value 
(say eps(min(x)) is so small that the vast majority of the other values 
will fail.  To illustrate, here is a test for Rik's changes to binocdf.m:

> > p = 1-1e-3;
> > pval = binopdf ([0:50]', 50, p);
> > assert (cumsum (pval), binocdf ([0:50]', 50, p), 1000 * eps (pval))

Notice that the comparisons are within tolerance for all values.  I've 
had to relax the tolerance a bit because as we discovered with the 
discussions about binocdf there are different numerical techniques at 
the heart of the two routines binopdf and binocdf.

The following is what results on binocdf.m prior to Rik's change:

> > assert (cumsum (pval), binocdf ([0:50]', 50, p), 1000 * eps (pval))
error: assert (cumsum (pval),binocdf ([0:50]', 50, p),1000 * eps (pval)) 
expected
   1.0000e-150
   4.9951e-146
[snip]
    2.0408e-09
    2.2198e-07
but got
   0.0000e+00
   0.0000e+00
[snip]
   2.0408e-09
   2.2198e-07
maximum absolute error 1e-150 exceeds tolerance 1.35666e-163
maximum absolute error 4.9951e-146 exceeds tolerance 8.89103e-159
[snip]
maximum absolute error 5.36716e-17 exceeds tolerance 4.1359e-22
maximum absolute error 1.89715e-17 exceeds tolerance 2.64698e-20

Notice that only the comparisons that failed are shown so that one 
doesn't have to weed through all the values.  Also notice that the 
comments about the error and tolerance now make sense, and furthermore 
the tolerance nicely scales relative to the magnitude of the expected 
item.  In this case, Rik knew about the string of zeros at the front. 
However, I think this sort of setup would be a good test for the 
plethora of other routines using 1 - betainc().

Dan