[racket] Rounding

On 03/16/2013 05:33 PM, Danny Yoo wrote:
>>
>> Since when is round 0.5 not giving 1?!
>
> According to the documentation:
>
>      http://docs.racket-lang.org/reference/generic-numbers.html#%28def._%28%28quote._~23~25kernel%29._round%29%29
>
>      "Returns the integer closest to x, resolving ties in favor of an
> even number... "
>
> And zero is an even number.  (http://www.bbc.co.uk/news/magazine-20559052)
>
> So although it may be surprising, at least the behavior is documented.
>
>
>
> Mathworld and Wikipedia say that this is done to combat statistical biasing:
>
>      http://mathworld.wolfram.com/NearestIntegerFunction.html
>      http://en.wikipedia.org/wiki/Rounding#Round_half_to_even

This is the reason the *low-order bit* of a floating-point number (or 
any finite-precision number) is rounded with ties toward even. Assuming 
that the unrepresented bits are distributed uniformly, it tends to make 
sums of values that are rounded multiple times more accurate.

IMO, the reason floating-point hardware rounds toward even for integer 
rounding is code reuse. Its bias correction is miniscule to nonexistent 
at that magnitude, unless you round after nearly every operation.

Here's a demo:


#lang typed/racket

(require math)

(: round/ties-toward-zero (Real -> Real))
;; Like the rounding we learned in school
(define (round/ties-toward-zero x)
   (if (x . < . 0)
       (- (round/ties-toward-zero (- x)))
       (floor (+ x 1/2))))

(: double-rounding-error ((Real -> Real) Real (Listof Real) -> Real))
;; Round every x once to the nearest `low-bit', then compare the sum
;; of those with the sum of them rounded to the nearest integer
(define (double-rounding-error round low-bit xs)
   (let ([xs  (map (λ: ([x : Real])
                     (* low-bit (round (/ x low-bit))))
                   xs)])
     (absolute-error (sum xs) (sum (map round xs)))))

;; Get some random numbers
(define xs (sample (uniform-dist 0 10) 10000))

;; Wow, this really helps if the second rounding is near the magnitude
;; of the low-order bit
(double-rounding-error round 1/2 xs)
(double-rounding-error round/ties-toward-zero 1/2 xs)

;; If the low-order bit is much smaller, though, it doesn't help at all
(double-rounding-error round 1/512 xs)
(double-rounding-error round/ties-toward-zero 1/512 xs)

;; Also, if we violate assumptions, it doesn't help as much (in this
;; case, the violated assumption is that there are as many values near
;; odd numbers as near even numbers)
(define ys (sample (uniform-dist 0 3) 10000))
(double-rounding-error round 1/2 ys)
(double-rounding-error round/ties-toward-zero 1/2 ys)

;; Violated assumption: missing bits are approximately uniformly
;; distributed (still does okay, though)
(define zs (map flsqrt (sample (uniform-dist 0 100) 10000)))
(double-rounding-error round 1/2 zs)
(double-rounding-error round/ties-toward-zero 1/2 zs)


Neil ⊥