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

From |
"Jann Ben" <ben.jann@soz.gess.ethz.ch> |

To |
<statalist@hsphsun2.harvard.edu> |

Subject |
RE: st: Another precision problem with Mata |

Date |
Mon, 22 Aug 2005 21:42:50 +0200 |

Many thanks for your explanations, Bill. The problem arose while implementing a kernel density estimator in Mata. There are kernel functions such as real matrix rd_kernel_rectangle(real matrix R) { return(1/2*(abs(R):<1)) } that are called from within another function as (*kern)((X[i]:-Y):/h) where Y is the data vector, h is the bandwidth (fixed or variable), and X[i] is the value at which the density be estimated. (Actually, the rectangular kernel is the only kernel that really has the problem I described because it is, well, rectangular. Most other kernels tend to zero the closer Y is to the border of the local data window.) It is common practice to choose min(Y)-h as the smallest X-value and max(Y)+h as the largest (see help -kdensity-) (note that this is somewhat at odds with the description of the problem in my last message: I should have written x = y + z instead of y = x + z). A common rule for choosing the bandwidth is h = 0.9 * sd / N^(1/5) or h = 0.9 * (IQR/1.349) / N^(1/5) where sd is the standard deviation of Y, IQR is the interquartile range and N is the sample size. So, this should give a rough idea of the relative magnitudes the values. Would this be enough information to "bound the error"? ben Bill wrote: > 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/

- Prev by Date:
**st: RE: RE: Turning Graph Legend Off** - Next by Date:
**st: unbalanced panel and nested logit** - Previous by thread:
**Re: st: Another precision problem with Mata** - Next by thread:
**Re: RE: st: Another precision problem with Mata** - Index(es):

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