Stata The Stata listserver
[Date Prev][Date Next][Thread Prev][Thread Next][Date index][Thread index]

RE: st: Another precision problem with Mata


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/



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