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

From: Ray Racine (ray.racine at gmail.com)
Date: Tue May 7 13:18:43 EDT 2013

And thanks for the overview on FP implementation and evaluation, very
instructive.  Saving this email for reference down the road.


On Tue, May 7, 2013 at 11:39 AM, Neil Toronto <neil.toronto at gmail.com>wrote:

> 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<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 ⊥
>
> ____________________
>  Racket Users list:
>  http://lists.racket-lang.org/**users <http://lists.racket-lang.org/users>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.racket-lang.org/users/archive/attachments/20130507/a87c9514/attachment.html>

Posted on the users mailing list.