[Date Prev][Date Next][Thread Prev][Thread Next][Date index][Thread index]
RE: st: Another precision problem with Mata
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)
that are called from within another function as
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)
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
Would this be enough information to "bound the error"?
> Ben Jann <email@example.com> 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)
> 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
* For searches and help try: