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

Re: RE: st: Another precision problem with Mata

From (William Gould, Stata)
Subject   Re: RE: st: Another precision problem with Mata
Date   Tue, 23 Aug 2005 09:55:31 -0500

Ben Jann <> 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

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

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


        We will explore 


        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)


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

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("   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("         (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)
*   For searches and help try:

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