[racket] ATTN: Neil T. Racket/Math digamma (psi) function implementation

On 05/06/2013 01:25 PM, Ray Racine wrote:
> Neil,
>
> Here in is an implementation of the digamma function.

Cool! Welcome to the High Priesthood of the... crazy people who do 
floating-point implementations. :)

BTW, I'm not clear on what you're requesting from me. Nearest I can come 
up with is how to measure error, so I'll address that below.

> Initially I tried
> Haskell's, then a Java impl from Mallet  (LDA library) but they didn't
> align (BTW I think the Haskell impl is flat wrong for digamma (2)) and
> sought alternate sources.
>
> I finally found
> http://code.google.com/p/pmtk3/source/browse/trunk/util/digamma.m?r=227&spec=svn227
> from pmtk2 a probabilistic modeling toolkit for Matlab/Octave, which is
> MIT LICENSED so I think we are good there.

You've discovered two things I found very frustrating while looking for 
code to copy for `math/special-functions' and `math/distributions':

1. Most implementations have ambiguous or terrible license terms.

2. Almost all implementations aren't good.

#1 made me write a lot of new code. #2 let me feel okay with that. :D

> My impl is based (well ok copied) from that impl.
>
> Octave also appears to have a Gnu Science Library FFI which has an impl
> for psi (digamma) as well which I've used as a sanity check for positive
> x values.

You can use `bfpsi0' from `math/bigfloat' to get arbitrarily accurate 
error measurements. I used it to produce error plots when implementing 
`flpsi0'. Demo below.

> Since it works for x = 10 to a few decimal places, I'm sure it works
> everywhere ... perfectly.

:D

> So far the Racket impl appears to be aligned to Octave+GSL, but I have
> no intuition as what is acceptable epsilon variance in floating point vs
> what is true error in implementation.  The impl handles z < 0, yet
> Octave+GSL blows up.

I set some pretty high standards for the math library. I call a flonum 
function implementation acceptable when:

  1. It gives sensible answers on the largest possible domain. (So
     Octave+GSL's digamma doesn't qualify.)

  2. Its error in ulps is <= 10 over most of the domain.

  3. When error in ulps is higher in some subdomain, it's because there
     are no reasonably fast algorithms for computing the function
     accurately in that subdomain.

Look up `flulp' and `flulp-error' in the math docs for definitions. 
Shortly, an ulp is the magnitude of the smallest change you can make to 
a flonum. Error in ulps is more or less the number of flonums between 
the true answer and the approximation.

A function with error <= 0.5 ulps is *correct*. I call a function with 
<= 2 ulps friggin' fantastic, and one with <= 10 accurate.

> I added addition poly coefficients up to 12 but the impl currently only
> uses 7.
>
> impl follows:  [...]

Paste this at the bottom to test it:


(require plot/typed
          math)

(: digamma* (Flonum -> Real))
(define (digamma* x)
   (parameterize ([bf-precision  53])
     (bigfloat->real (bfpsi0 (bf x)))))

(: plot-error ((Flonum -> Flonum) (Flonum -> Real) Real Real -> Any))
(define (plot-error f-appx f-truth x-min x-max)
   (plot (function
          (λ: ([x : Real])
            (let ([x  (fl x)])
              (define err (flulp-error (f-appx x) (f-truth x)))
              (if (rational? err) err -1.0)))
          x-min x-max)))

(plot-error flpsi0 digamma* 0 32)
(plot-error digamma digamma* 0 32)

(plot-error flpsi0 digamma* -8 0)
(plot-error digamma digamma* -8 0)


It looks like your `digamma' is as accurate as `flpsi0' for x > 21 or so 
(i.e. <= 2 ulps, which is friggin' fantastic). If using 12 terms instead 
of 7 in your polynomial approximation reduced the error for smaller x, 
it would be worth using that polynomial approximation in the math 
library. (I can't test this because I don't know the coefficients' 
signs.) In the subdomain it would likely be accurate in, `flpsi0' uses 
12- to 15-degree Chebyshev polynomials, which are about twice as slow to 
evaluate as similar-degree polynomials using Horner's rule (as in your 
`digamma-poly').

Neil ⊥