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

Well crap ...

I think I searched the Racket doc for digamma found nada, tried psi and
blew it as I never progressed beyond "implements a polygamma function" and
failed to recognize it for what it was.    In other words, I flat missed
it.  So I was sort of asking about the process of voir dire'ing my impl
with the idea of addition to the standard Racket math libraries, thinking
it was absent, but ... never mind :).

Worthwhile I think to get 'digamma' to show up as a hit in the Racket doc
search and maybe extending the note in the psi section, "Note that (psi 0
x) = (psi0 x)" to maybe "Note that (psi 0 x) = (psi0 x) is the digamma
function" for when someone searches on 'psi'.




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/c82068e0/attachment-0001.html>