Stata The Stata listserver
[Date Prev][Date Next][Thread Prev][Thread Next][Date index][Thread index]

Re: st: Another precision problem with Mata


From   wgould@stata.com (William Gould, Stata)
To   statalist@hsphsun2.harvard.edu
Subject   Re: st: Another precision problem with Mata
Date   Mon, 22 Aug 2005 09:55:42 -0500

Ben Jann <ben.jann@soz.gess.ethz.ch> 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
wgould@stata.com
*
*   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/



© Copyright 1996–2014 StataCorp LP   |   Terms of use   |   Privacy   |   Contact us   |   What's new   |   Site index