Machine precision strategies in Mata

[John Mullahy]

I'm curious what might be your favorite strategies for handling machine precision issues like those illustrated by the following example. The ultimate goal would be to work with quantities like r but, once computed, "manage" them so they equal .025. Thanks for sharing any suggestions, hints, or tricks that work for you.

mata

p=.95
q=.025
r=.5*(1-p)

q==r
40*q==1
40*r==1

end

: p=.95

: q=.025

: r=.5*(1-p)

:
: q==r
  0

: 40*q==1
  1

: 40*r==1
  0

[daniel klein]

If I remember correctly, Bill Gould talks about rounding errors in calculations in his Mata Book and how they arise only(?) in addition and subtraction of non-integers. In your example, you could calculate r as

r = .5 * ( (100-(100*p))/100 )

subtracting integers.

[Niels Henrik Bruun]

I don't think it is possible always to maintain precision during calculations.
So my suggestion is to use the absolute difference:

: p=.95
: q=.025
: r=.5*(1-p)

: q==r
  0
: abs(q-r), abs(q-r) < 1e-10
                 1             2
    +-----------------------------+
  1 |  2.08167e-17             1  |
    +-----------------------------+
: 40*q==1
  1
: 40*r==1
  0
: abs(40*r-1), abs(40*r-1) < 1e-10
                 1             2
    +-----------------------------+
  1 |  8.88178e-16             1  |
    +-----------------------------+
: end

[John Mullahy]

Thanks Daniel and Niels for the helpful insights.

[daniel klein]

Some functions have quad*() variations that increase precision during the calculation; might be another useful tool:

. mata
------------------------------------------------- mata (type end to exit) -----------------------------------------------------------------------------------------------------------------
: x = J(1, 6, .025)

: rowsum(x) == 6*.025
  0

: quadrowsum(x) == 6*.025
  1

: end
-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

Edit:

Matrix-algebra might also help:

: x*J(6,1,1) == 6*.025
  1

[daniel klein]

So my suggestion is to use the absolute difference:

I guess the approach is fine. There might be a reason related to the same underlying problems of precision to prefer relative difference:

. mata
------------------------------------------------- mata (type end to exit) -----------------------------------------------------------------------------------------------------------------
: p=.95

: q=.025

: r=.5*(1-p)

:
: abs(q-r)
  2.08167e-17

: reldif(q, r)
  2.03090e-17

: end
-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

I am not sure whether that ever makes a difference in real-world problems.

Edit:

Looking at the formula, the relative difference, as defined in reldif(), appears to be relevant when values get close to 0. But I am really guessing here and I will stop to do so.

[John Mullahy]

I just came upon the undocumented _signif() function (help mf__signif; note two underscores). Note its use in the following (an extension of the code in #1).

mata

p=.95
q=.025
r=.5*(1-p)
ra=.5*(1-p)
_signif(ra,15)

q==r
q==ra

40*q==1
40*r==1
40*ra==1

end

Results:

:
: p=.95

: q=.025

: r=.5*(1-p)

: ra=.5*(1-p)

: _signif(ra,15)

:
: q==r
  0

: q==ra
  1

:
: 40*q==1
  1

: 40*r==1
  0

: 40*ra==1
  1

[Hua Peng (StataCorp)]

You may also try:

. mata:
------------------------------------------------- mata (type end to exit) ----------------------------
: p = 0.95

: q = 1 - p

: printf("%21x", q)
+1.99999999999a0X-005
: r = 0.05

: printf("%21x", r)
+1.999999999999aX-005
: end
------------------------------------------------------------------------------------------------------