Ben Jann <[email protected]> asked about the precision of abs((x-y)/z)
and in particular, abs((x-y)/z)<1 when x=y+z. Actually, yesterday he asked
the question when y=x+z, but one way or the other, it does not make much
difference.
I replied, rather negatively and certainly in purist mode. I said that there
was no reason to think about how one might code abs((x-y/z)<1 so that it is
guaranteed to be false in all the right circumstances because, for all
possible values, there is no way to do that.
I could have been more helpful, however, by eliminating cases when
calculating y PLUS z results in y or results in z because the other is too
small.
I have now thought about that problem, and I believe the solution to
Ben's problem is to code
if ( abs(abs(x-y)-z) < epsilon(x) ) {
...
}
The above is not guaranteed to work, but it will work in lots of
reasonable cases.
Below my signature is a do-file that Ben should run. It provides a
demonstration of my claim and it shows exactly what the problem is
and where the precision is being lost.
-- Bill
[email protected]
------------------------------------------------------------------------------
clear
/*
We will explore
abs((x-y)/z)<1
when x = y+z. In that case, the expression becomes
abs(((y+z)-y)/z) < 1)
The correct answer is 0 because in infinite precision
abs(((y+z)-y)/z) < 1)
abs(z/z) < 1)
abs(1 < 1)
Below, we will look at two expressions:
x - y (which is (y+z)-y in our special case)
abs((x-y)/z)
*/
mata:
real scalar f(y, z)
{
real scalar x
x = y + z
return(abs((x-y)/z) < 1)
}
f(0, .1) // no problem
f(1, .1)
f(1000, .1) // no problem
f(1e+20, .1) // problem
real scalar findprob(y, z, mult)
{
while (f(y,z)==0) y = mult*y
return(y)
}
findprob(1, .1, 2)
f(4, .1) // problem at 4
void fdetail(y, z)
{
real scalar x
printf(" y = %21x\n", y)
printf(" z = %21x\n", z)
printf("-----------------------------------\n")
printf(" x = y+z = %21x\n", x = y+z)
printf(" x - y = %21x (%s)\n", x-y,
(x-y)==z ? "okay, x-y == z" : "problem, x-y != z")
printf("\n")
printf("error:\n")
printf(" (x-y)-z = %g\n", (x-y)-z)
printf(" |(x-y)/z| - 1 = %g\n", abs((x-y)/z)-1)
printf(" note, e(y)=%g, e(x)=%g e(z)=%g e(1)=%g\n",
epsilon(y), epsilon(x), epsilon(z), epsilon(1))
}
fdetail(4, .1)
fdetail(16384, .1)
fdetail(32768, .1)
fdetail(65536, .1)
fdetail(131072, .1)
fdetail(262144, .1)
fdetail(1e+10, .1)
fdetail(1e+20, .1) // complete loss of precision
/* ----- the point of all the above is that
if y>z, then (x-y) has error bounded by epsilon(x).
assuming numbers are such that we have any precision at all.
*/
/* ----------------------------------- */
// now consider y<z:
f(.1, 1)
f(.1, 2)
f(.1, 4)
fdetail(.1, 4)
fdetail(.1, 512)
fdetail(.1, 1024)
fdetail(.1, 8192)
fdetail(.1, 16384)
fdetail(.1, 32768)
fdetail(.1, 262144)
end
------------------------------------------------------------------------------
*
* 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/
