[Python-ideas] Way to check for floating point "closeness"?

[Steven D'Aprano]

On Wed, Jan 14, 2015 at 08:13:42PM -0600, Ron Adam wrote:
> 
> 
> On 01/14/2015 06:04 PM, Steven D'Aprano wrote:
[...]
> >Of course doing division by zero is an error. The whole point of this
> >discussion is that you cannot talk about errors relative to zero. At
> >zero, you can still have an absolute error, but you cannot have a
> >relative error. The definition of the error between two numeric values:
> >
> >def abserr(a, b):
> >     return abs(a - b)
> >
> >def relerr(a, b):
> >     return abs(a - b)/min(abs(a), abs(b))
> >
> >
> >If either a or b is zero, relarr() fails. This is a mathematical
> >limitation, not a bug to be fixed.
> 
> 
> I would probably write like this.
> 
>     def relerr(actual, expected):
>         return abs(actual - expected)/expected
> 
> Where expected is a non-zero reference value.
> 
> it needs a third value.

No the function doesn't require a third value. It just returns the 
relative error between two values. What you do with that error is up to 
you:

- compare it to an acceptable amount of error, as you do below
- print it
- add it to a list of error amounts
etc.

Of course, you could have other functions which do more, like decide 
whether an error is "small enough", but that is out of scope for a 
function which only calculates the error.


> In the above the expected value is the scaler.
> And then use it this way to measure a resistors tolerance.
> 
>     if relerr(218.345, 220) <= 0.05:
>         print("Good")
>     else:
>         print("Bad")

# using your implementation
if relerr(actual=200.0, expected=-5.0) < 0.05:
    print("oops")

Errors, by definition, are non-negative. You have to take care to ensure 
that the error value returned is non-negative.

The question of which to use as the denominator is more subtle. Like 
you, I used to think that you should choose ahead of time which value 
was expected and which was actual, and divide by the actual. Or should 
that be the expected? I could never decide which I wanted: error 
relative to the expected, or error relative to the actual? And then I 
could never remember which order the two arguments went.

Finally I read Bruce Dawson (I've already linked to his blog three or 
four times) and realised that he is correct and I was wrong. Error 
calculations should be symmetrical, so that

    error(a, b) == error(b, a)

regardless of whether you have absolute or relative error. Furthermore, 
for safety you normally want the larger estimate of error, not the 
smaller: given the choice between

    (abs(a - b))/abs(a)

versus

    (abs(a - b))/abs(b)


you want the *larger* error estimate, which means the *smaller* 
denominator. That's the conservative way of doing it.

A concrete example: given a=5 and b=7, we have:

absolute error = 2
relative error (calculated relative to a) = 0.4
relative error (calculated relative to b) = 0.286

That is, b is off by 40% relative to a; or a is off by 28.6% relative to 
b. Or another way to put it, given that a is the "true" value, b is 40% 
too big; or if you prefer, 28.6% of b is in error.

Whew! Percentages are hard! *wink*

The conservative, "safe" way to handle this is to just treat the error 
function as symmetrical and always report the larger of the two relative 
errors (excluding the case where the denominator is 0, in which case 
the relative error is either 100% or it doesn't exist). Worst case, you 
may reject some values which you should accept, but you will never 
accept any values that you should reject.


> Note that you would never compare to an expected value of zero.

You *cannot* compare to an expected value of zero, but you certainly can 
be in a situation where you would like to: math.sin(math.pi) should 
return 0.0, but doesn't, it returns 1.2246063538223773e-16 instead. What 
is the relative error of the sin function at x = math.pi?


>     relerr(a - b, expected_feet) < tolerance   # relative feet from b
>     relerr(a - 0, expected_feet) < tolerance   # relative feet from zero
>     relerr(a - b, ulp)    # percentage of ulp's

I don't understand what you think these three examples are showing.


> What Chris is looking for is a way to get a closeness function that works 
> most of the time. (He posted while I'm writing this.)

I think the function I have in the statistics test suite is that 
function. I would like to see ULP calculations offered as well, but 
Mark thinks that's unnecessary and I'm not going to go to the 
battlements to fight for ULPs.


> But I don't see how you can do that without specifying some scaler to give 
> it context, and a tolerance.

By the way, it is spelled "scalar", and all that means is that it is a 
number, like 23, not a vector or array or matrix.

>      def is_close(a, b, unit, good):
>          return (abs(a - b) / unit) <= good

Take a look at the statistics test suite. I'll be the first to admit 
that the error tolerances are plucked from thin air, based on what I 
think are "close enough", but they show how such a function might work:

* you provide two values, and at least one of an absolute error 
  tolerance and a relative error; 
* if the error is less than the error(s) you provided, the test 
  passes, otherwise it fails;
* NANs and INFs are handled apprpriately.


>      is_close(218.345, 220, 1, .05)   # OHMs
>      is_close(a, b, ULP, 2)     # ULPs
>      is_close(a, b, AU, .001)   # astronomical units
> 
> 
> I don't see anyway to generalise those with just a function.

Generalise in what way?


> By using objects we can do a bit more.  I seem to recall coming across 
> measurement objects some place.  They keep a bit more context with them.

A full system of <value + unit> arithmetic is a *much* bigger problem 
than just calculating error estimates correctly, and should be a 
third-party library before even considering it for the std lib.


-- 
Steve