Announcement

Collapse
No announcement yet.
X
  • Filter
  • Time
  • Show
Clear All
new posts

  • Machine precision strategies in Mata

    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.
    Code:
    mata
    
    p=.95
    q=.025
    r=.5*(1-p)
    
    q==r
    40*q==1
    40*r==1
    
    end
    Code:
    : p=.95
    
    : q=.025
    
    : r=.5*(1-p)
    
    :
    : q==r
      0
    
    : 40*q==1
      1
    
    : 40*r==1
      0

  • #2
    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

    Code:
    r = .5 * ( (100-(100*p))/100 )
    subtracting integers.

    Comment


    • #3
      I don't think it is possible always to maintain precision during calculations.
      So my suggestion is to use the absolute difference:
      Code:
      : 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
      Kind regards

      nhb

      Comment


      • #4
        Thanks Daniel and Niels for the helpful insights.

        Comment


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

          Code:
          . 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:

          Code:
          : x*J(6,1,1) == 6*.025
            1
          Last edited by daniel klein; 08 May 2022, 23:25.

          Comment


          • #6
            Originally posted by Niels Henrik Bruun View Post
            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:

            Code:
            . 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.
            Last edited by daniel klein; 08 May 2022, 23:39.

            Comment


            • #7
              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).
              Code:
              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:
              Code:
              :
              : 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

              Comment


              • #8
                You may also try:

                Code:
                . 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
                ------------------------------------------------------------------------------------------------------

                Comment

                Working...
                X