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   wgould@stata.com (William Gould, Stata)
To   statalist@hsphsun2.harvard.edu
Subject   Re: RE: st: Another precision problem with Mata
Date   Tue, 23 Aug 2005 09:55:31 -0500

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

------------------------------------------------------------------------------
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/



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