# RE: st: Another precision problem with Mata

 From "Jann Ben" To 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:
>
> > 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 subtraction
>
> 		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
> 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/
```