[Date Prev][Date Next][Thread Prev][Thread Next][Date index][Thread index]

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/

- Prev by Date:
**st: Multinomial logits with proportions for dependent variables** - Next by Date:
**Re: st: Heterogeneity of variances** - Previous by thread:
**RE: st: Another precision problem with Mata** - Next by thread:
**st: Calculating Standard Errors of Predicted Probabilities in ProbitModels** - Index(es):

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