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 <[email protected]> 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
> [email protected]
*
* 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/
