Ben Jann <[email protected]> asks,
> in a Mata program I am using the expression
>
> abs((x-y)/z) < 1
>
> where x, y and z are real scalars. Note that in some cases y has been
> set to x + z earlier in the program. Although the expression should
> evaluate to "false" in these cases, it sometimes does evaluate to "true"
> due to roundoff error.
>
> My question now is: Would it be reasonable to code
>
> abs((x-y)/z) < 1-epsilon(1)
>
> Or how would one approach such a problem?
The first step in approaching this problem is to understand the source
of the error. Ben reports that when y = x + z, he can have problems.
In the y=x+z case,
abs((x-y)/z) = abs((x-(x+z))/z)
= abs(-z/z) (inifinite precision)
= 1 (inifinite precision)
Ben is hoping that x-(x+z) == z, and that is not necessarily true when
we perform finite-precision arithmetic.
So now let's consider (x-(x+z)), the source of the problem.
I am going to use the following notation
+ infinite-precision addition
- infinite-precision subtraction
PLUS finite-precision addition
SUB finite-precision subtraction
In this case
x PLUS z = x+z + error
The error can be so large so that
x PLUS z = x (1)
or
x PLUS z = z (2)
All that is required is that (1) x be very large relative to z, or
(2) vice versa.
The importance of that point is that we cannot bound the error by
epsilon(1), as Ben was hoping.
I cannot go further with this problem, in fact, without knowing something
about the relative magnitudes of x and z.
However, Ben has just revealed that it is important to him to identify
abs((x-y)/z) == 1 when y = x+z. In that case, I suggest Ben remove y from his
code and substitue (d+z) everywhere for it. Then abs((x-y)/z) == 1 is
equivalent to abs(d)==abs(x).
Perhaps there are other possible code-specific solutions, or perhaps,
Ben does know something about magnitudes and we can bound the error, and
Ben will be willing to accept abs((x-y)/z)<1 will evaluate to false in
some cases where, in infinite precision, it would evaluate to true.
But one thing I can guarantee: There is no reasons to think about how
one might code
abs((x-y)/z) < 1
so that it is guaranteed to be false when y = x PLUS z, for all values
of x and z.
-- Bill
[email protected]
*
* For searches and help try:
* http://www.stata.com/support/faqs/res/findit.html
* http://www.stata.com/support/statalist/faq
* http://www.ats.ucla.edu/stat/stata/
