No subject
From: ()
Date: Mon Dec 3 19:58:15 EST 2012 |
|
|#
;; Assumes for z < 0
;; Use the reflection formula (Jeffrey 11.1.6):
;; digamma(-x) =3D digamma(x+1) + pi*cot(pi*x)
;; y =3D digamma(-x+1) + pi*cot(-pi*x);
;; This is related to the identity
;; digamma(-x) =3D digamma(x+1) - digamma(z) + digamma(1-z)
;; where z is the fractional part of x
;; For example:
;; digamma(-3.1) =3D 1/3.1 + 1/2.1 + 1/1.1 + 1/0.1 + digamma(1-0.1)
;; =3D digamma(4.1) - digamma(0.1) + digamma(1-0.1)
;; Then we use
;; digamma(1-z) - digamma(z) =3D pi*cot(pi*z)
;; cot(x) =3D tan (pi/2 - x)
(provide:
[digamma (Float -> Float)])
(require
racket/match
racket/pretty
(only-in racket/flonum
fltan
fllog)
(only-in racket/math
pi))
(define: PI : Float pi)
(define: NEG-PI : Float (- PI))
(define: PI-OVER-TWO : Float (/ PI 2.0))
(define: PI-SQUARED-OVER-SIX : Float (/ (* PI PI) 6))
(define: EULER-MASCHERONI : Float -0.5772156649015328606065121)
(define: DIGAMMA-1 : Float EULER-MASCHERONI) ;; digamma(1)
(define: DIGAMMA-2 : Float PI-SQUARED-OVER-SIX)
(define: DIGAMMA-COEF-2 : Float (- (/ 1.0 12.0)))
(define: DIGAMMA-COEF-3 : Float (/ 1.0 120.0))
(define: DIGAMMA-COEF-4 : Float (- (/ 1.0 252.0)))
(define: DIGAMMA-COEF-5 : Float (/ 1.0 240.0))
(define: DIGAMMA-COEF-6 : Float (- (/ 1.0 132.0)))
(define: DIGAMMA-COEF-7 : Float (/ 691.0 32760.0))
(define: DIGAMMA-COEF-8 : Float (- (/ 1.0 12.0)))
(define: DIGAMMA-LARGE : Float 9.5)
(define: DIGAMMA-SMALL : Float 0.000001)
;;anxn + =E2=80=A6 + a2x2 + a1x + a0 =3D a0 + x(a1 + x(a2 + x(=E2=80=A6 an)=
=E2=80=A6) )
(: digamma-poly (Float -> Float))
(define (digamma-poly z)
(* z (+ DIGAMMA-COEF-2
(* z (+ DIGAMMA-COEF-3
(* z (+ DIGAMMA-COEF-4
(* z (+ DIGAMMA-COEF-5
(* z (+ DIGAMMA-COEF-6
(* z (+ DIGAMMA-COEF-7
(* z (+
DIGAMMA-COEF-8)))))))))))))))
(: cotan (Float -> Float))
(define (cotan x)
(fltan (- PI-OVER-TWO x)))
(: digamma-negative (Float -> Float))
(define (digamma-negative z)
;; (assert (< z 0.0))
(+ (digamma (+ (- z) 1.0))
(* PI (cotan (* z NEG-PI)))))
(: digamma-shift (Float -> (Values Float Float)))
(define (digamma-shift z)
(let loop ((z z) (y 0.0))
(if (>=3D z DIGAMMA-LARGE)
(values z y)
(loop (+ z 1.0) (- y (/ z))))))
(: digamma (Float -> Float))
(define (digamma z)
(cond
((eq? z 0.0)
-inf.0)
((< z 0.0)
(digamma-negative z))
((eq? z 1.0)
DIGAMMA-1)
((<=3D z DIGAMMA-SMALL)
(/ (- DIGAMMA-1 1.0)
(+ z (* DIGAMMA-2 z))))
(else
(let-values (((z z-shift) (digamma-shift z))) ;; Reduce x < large to
digamma(x+n) where (x+n) >=3D large
(let ((r (/ z))) ;; Use de Moivre's
expansion if argument >=3D large.
(let ((A (+ z-shift (- (fllog z) (* 0.5 r)))))
(+ A (digamma-poly (* r r)))))))))
(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
(=CE=BB: ([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 30)
(plot-error digamma digamma* 0 30)
(plot-error flpsi0 digamma* -8 0)
(plot-error digamma digamma* -8 0)
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=3D227&spec=3D**svn227<http://code.google.com/p/pmtk3/so=
urce/browse/trunk/util/digamma.m?r=3D227&spec=3Dsvn227>
>> 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 =3D 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 <=3D 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 <=3D 0.5 ulps is *correct*. I call a function with =
<=3D
> 2 ulps friggin' fantastic, and one with <=3D 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
> (=CE=BB: ([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. <=3D 2 ulps, which is friggin' fantastic). If using 12 terms instea=
d 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 =E2=8A=A5
>
> ____________________
> Racket Users list:
> http://lists.racket-lang.org/**users <http://lists.racket-lang.org/users=
>
>
--047d7b3a9684d43f1e04dc27158c
Content-Type: text/html; charset=UTF-8
Content-Transfer-Encoding: quoted-printable
<div dir=3D"ltr"><div><div>Well since I gone this far, I took the unusual s=
tep of implementing the algorithm (more?) correctly this time.<br><br>On th=
e Toronto scale, we are mostly accurate for most of the domain=C2=A0 x >=
0.=C2=A0 <br>
<br>For x<0 the use of psi(1 - x) - psi (x) =3D pi * cot(pi * x), not so=
good.<br><br></div>For x > 0, outside the subdomain around the circle x=
=3D1.48 where things are a bit wobbly it is pretty solid.=C2=A0 I don't=
know if this is a zero of the polynomial or what the heck is going on ther=
e.<br>
<br></div><div>;; Improved Impl of digamma<br></div><div><div><br>#lang typ=
ed/racket/base<br><br>#|<br>Reference:<br><br>J Bernardo,<br>Psi ( Digamma =
) Function,<br>Algorithm AS 103,<br>Applied Statistics,<br>Volume 25, Numbe=
r 3, pages 315-317, 1976.<br>
<br>From <a href=3D"http://www.psc.edu/~burkardt/src/dirichlet/dirichlet.f"=
>http://www.psc.edu/~burkardt/src/dirichlet/dirichlet.f</a><br>|#<br><br><b=
r>;; Assumes for z < 0<br>;;=C2=A0 Use the reflection formula (Jeffrey 1=
1.1.6):<br>
;;=C2=A0 digamma(-x) =3D digamma(x+1) + pi*cot(pi*x)<br>;;=C2=A0=C2=A0 y =
=3D digamma(-x+1) + pi*cot(-pi*x);<br>;;=C2=A0 This is related to the ident=
ity<br>;;=C2=A0 digamma(-x) =3D digamma(x+1) - digamma(z) + digamma(1-z)<br=
>;;=C2=A0 where z is the fractional part of x<br>
;;=C2=A0 For example:<br>;;=C2=A0 digamma(-3.1) =3D 1/3.1 + 1/2.1 + 1/1.1 +=
1/0.1 + digamma(1-0.1)<br>;;=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=
=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0 =3D digamma(4.1) - digamma(0.=
1) + digamma(1-0.1)<br>;;=C2=A0 Then we use<br>;;=C2=A0 digamma(1-z) - diga=
mma(z) =3D pi*cot(pi*z)<br>
;;=C2=A0 cot(x) =3D tan (pi/2 - x)<br><br>(provide:<br>=C2=A0[digamma (Floa=
t -> Float)])<br><br>(require<br>=C2=A0racket/match<br>=C2=A0racket/pret=
ty<br>=C2=A0(only-in racket/flonum<br>=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=
=C2=A0=C2=A0=C2=A0 fltan<br>=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=
=A0=C2=A0 fllog)<br>=C2=A0(only-in racket/math<br>
=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0 pi))<br><br>(define:=
PI : Float pi)<br>(define: NEG-PI : Float (- PI))<br>(define: PI-OVER-TWO =
: Float (/ PI 2.0))<br>(define: PI-SQUARED-OVER-SIX : Float (/ (* PI PI) 6)=
)<br>(define: EULER-MASCHERONI : Float -0.5772156649015328606065121)<br>
<br>(define: DIGAMMA-1 : Float EULER-MASCHERONI)=C2=A0=C2=A0=C2=A0=C2=A0 ;;=
digamma(1)<br>(define: DIGAMMA-2 : Float PI-SQUARED-OVER-SIX)<br><br>(defi=
ne: DIGAMMA-COEF-2 : Float (- (/ 1.0 12.0)))<br>(define: DIGAMMA-COEF-3 : F=
loat (/ 1.0 120.0))<br>
(define: DIGAMMA-COEF-4 : Float (- (/ 1.0 252.0)))<br>(define: DIGAMMA-COEF=
-5 : Float (/ 1.0 240.0))<br>(define: DIGAMMA-COEF-6 : Float (- (/ 1.0 132.=
0)))<br>(define: DIGAMMA-COEF-7 : Float (/ 691.0 32760.0))<br>(define: DIGA=
MMA-COEF-8 : Float (- (/ 1.0 12.0)))<br>
<br>(define: DIGAMMA-LARGE : Float 9.5)<br>(define: DIGAMMA-SMALL : Float 0=
.000001)<br><br>;;anxn + =E2=80=A6 + a2x2 + a1x + a0 =3D a0 + x(a1 + x(a2 +=
x(=E2=80=A6 an) =E2=80=A6) )<br>(: digamma-poly (Float -> Float))<br>(d=
efine (digamma-poly z)<br>
=C2=A0 (* z (+ DIGAMMA-COEF-2 <br>=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=
=A0=C2=A0=C2=A0 (* z (+ DIGAMMA-COEF-3<br>=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=
=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0 (* z =
(+ DIGAMMA-COEF-4<br>=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=
=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=
=A0=C2=A0=C2=A0=C2=A0 (* z (+ DIGAMMA-COEF-5<br>=C2=A0=C2=A0=C2=A0=C2=A0=C2=
=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=
=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=
=A0=C2=A0=C2=A0=C2=A0 (* z (+ DIGAMMA-COEF-6<br>
=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=
=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=
=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=
=A0=C2=A0=C2=A0=C2=A0 (* z (+ DIGAMMA-COEF-7<br>=C2=A0=C2=A0=C2=A0=C2=A0=C2=
=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=
=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=
=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=
=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0 (* z (+ DIGAMMA-COEF-8))))))))))=
)))))<br><br>(: cotan (Float -> Float))<br>(define (cotan x)<br>=C2=A0 (=
fltan (- PI-OVER-TWO x)))<br>
<br>(: digamma-negative (Float -> Float))<br>(define (digamma-negative z=
)<br>=C2=A0 ;; (assert (< z 0.0))<br>=C2=A0 (+ (digamma (+ (- z) 1.0))<b=
r>=C2=A0=C2=A0=C2=A0=C2=A0 (* PI (cotan (* z NEG-PI)))))<br><br>(: digamma-=
shift (Float -> (Values Float Float)))<br>
(define (digamma-shift z)<br>=C2=A0 (let loop ((z z) (y 0.0))<br>=C2=A0=C2=
=A0=C2=A0 (if (>=3D z DIGAMMA-LARGE)<br>=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=
=C2=A0=C2=A0 (values z y)<br>=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0 (lo=
op (+ z 1.0) (- y (/ z))))))<br><br>(: digamma (Float -> Float))<br>(def=
ine (digamma z)<br>
=C2=A0 (cond <br>=C2=A0=C2=A0=C2=A0 ((eq? z 0.0)<br>=C2=A0=C2=A0=C2=A0=C2=
=A0 -inf.0) <br>=C2=A0=C2=A0=C2=A0 ((< z 0.0)<br>=C2=A0=C2=A0=C2=A0=C2=
=A0 (digamma-negative z))<br>=C2=A0=C2=A0=C2=A0 ((eq? z 1.0)<br>=C2=A0=C2=
=A0=C2=A0=C2=A0 DIGAMMA-1)<br>=C2=A0=C2=A0=C2=A0 ((<=3D z DIGAMMA-SMALL)=
<br>=C2=A0=C2=A0=C2=A0=C2=A0 (/ (- DIGAMMA-1 1.0)<br>=C2=A0=C2=A0=C2=A0=C2=
=A0=C2=A0=C2=A0=C2=A0 (+ z (* DIGAMMA-2 z))))<br>
=C2=A0=C2=A0=C2=A0 (else=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=
=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0 <br>=
=C2=A0=C2=A0=C2=A0=C2=A0 (let-values (((z z-shift) (digamma-shift z))) ;; R=
educe x < large to digamma(x+n) where (x+n) >=3D large=C2=A0=C2=A0=C2=
=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=
=C2=A0 <br>=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0 (let ((r (/ z)))=C2=A0=C2=
=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=
=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=
=A0 ;; Use de Moivre's expansion if argument >=3D large.<br>
=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0 (let ((A (+ z-shift (- (fl=
log z) (* 0.5 r)))))<br>=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=
=A0=C2=A0 (+ A (digamma-poly (* r r)))))))))<br><br><br>(require plot/typed=
<br>=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0 math)<br><br>(: digamm=
a* (Flonum -> Real))<br>(define (digamma* x)<br>
=C2=A0 (parameterize ([bf-precision=C2=A0 53])<br>=C2=A0=C2=A0=C2=A0 (bigfl=
oat->real (bfpsi0 (bf x)))))<br><br>(: plot-error ((Flonum -> Flonum)=
(Flonum -> Real) Real Real -> Any))<br>(define (plot-error f-appx f-=
truth x-min x-max)<br>
=C2=A0 (plot (function<br>=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0 =
(=CE=BB: ([x : Real])<br>=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=
=C2=A0=C2=A0 (let ([x=C2=A0 (fl x)])<br>=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=
=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0 (define err (flulp-error (f-appx x)=
(f-truth x)))<br>=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=A0=C2=
=A0=C2=A0=C2=A0 (if (rational? err) err -1.0)))<br>=C2=A0=C2=A0=C2=A0=C2=A0=
=C2=A0=C2=A0=C2=A0=C2=A0 x-min x-max)))<br>
<br>(plot-error flpsi0 digamma* 0 30)<br>(plot-error digamma digamma* 0 30)=
<br><br>(plot-error flpsi0 digamma* -8 0)<br>(plot-error digamma digamma* -=
8 0)<br><br></div></div></div><div class=3D"gmail_extra"><br><br><div class=
=3D"gmail_quote">
On Tue, May 7, 2013 at 11:39 AM, Neil Toronto <span dir=3D"ltr"><<a href=
=3D"mailto:neil.toronto at gmail.com" target=3D"_blank">neil.toronto at gmail.com=
</a>></span> wrote:<br><blockquote class=3D"gmail_quote" style=3D"margin=
:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
<div class=3D"im">On 05/06/2013 01:25 PM, Ray Racine wrote:<br>
<blockquote class=3D"gmail_quote" style=3D"margin:0 0 0 .8ex;border-left:1p=
x #ccc solid;padding-left:1ex">
Neil,<br>
<br>
Here in is an implementation of the digamma function.<br>
</blockquote>
<br></div>
Cool! Welcome to the High Priesthood of the... crazy people who do floating=
-point implementations. :)<br>
<br>
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.<div =
class=3D"im"><br>
<br>
<blockquote class=3D"gmail_quote" style=3D"margin:0 0 0 .8ex;border-left:1p=
x #ccc solid;padding-left:1ex">
Initially I tried<br>
Haskell's, then a Java impl from Mallet =C2=A0(LDA library) but they di=
dn't<br>
align (BTW I think the Haskell impl is flat wrong for digamma (2)) and<br>
sought alternate sources.<br>
<br>
I finally found<br>
<a href=3D"http://code.google.com/p/pmtk3/source/browse/trunk/util/digamma.=
m?r=3D227&spec=3Dsvn227" target=3D"_blank">http://code.google.com/p/<u>=
</u>pmtk3/source/browse/trunk/<u></u>util/digamma.m?r=3D227&spec=3D<u><=
/u>svn227</a><br>
from pmtk2 a probabilistic modeling toolkit for Matlab/Octave, which is<br>
MIT LICENSED so I think we are good there.<br>
</blockquote>
<br></div>
You've discovered two things I found very frustrating while looking for=
code to copy for `math/special-functions' and `math/distributions'=
:<br>
<br>
1. Most implementations have ambiguous or terrible license terms.<br>
<br>
2. Almost all implementations aren't good.<br>
<br>
#1 made me write a lot of new code. #2 let me feel okay with that. :D<div c=
lass=3D"im"><br>
<br>
<blockquote class=3D"gmail_quote" style=3D"margin:0 0 0 .8ex;border-left:1p=
x #ccc solid;padding-left:1ex">
My impl is based (well ok copied) from that impl.<br>
<br>
Octave also appears to have a Gnu Science Library FFI which has an impl<br>
for psi (digamma) as well which I've used as a sanity check for positiv=
e<br>
x values.<br>
</blockquote>
<br></div>
You can use `bfpsi0' from `math/bigfloat' to get arbitrarily accura=
te error measurements. I used it to produce error plots when implementing `=
flpsi0'. Demo below.<div class=3D"im"><br>
<br>
<blockquote class=3D"gmail_quote" style=3D"margin:0 0 0 .8ex;border-left:1p=
x #ccc solid;padding-left:1ex">
Since it works for x =3D 10 to a few decimal places, I'm sure it works<=
br>
everywhere ... perfectly.<br>
</blockquote>
<br></div>
:D<div class=3D"im"><br>
<br>
<blockquote class=3D"gmail_quote" style=3D"margin:0 0 0 .8ex;border-left:1p=
x #ccc solid;padding-left:1ex">
So far the Racket impl appears to be aligned to Octave+GSL, but I have<br>
no intuition as what is acceptable epsilon variance in floating point vs<br=
>
what is true error in implementation. =C2=A0The impl handles z < 0, yet<=
br>
Octave+GSL blows up.<br>
</blockquote>
<br></div>
I set some pretty high standards for the math library. I call a flonum func=
tion implementation acceptable when:<br>
<br>
=C2=A01. It gives sensible answers on the largest possible domain. (So<br>
=C2=A0 =C2=A0 Octave+GSL's digamma doesn't qualify.)<br>
<br>
=C2=A02. Its error in ulps is <=3D 10 over most of the domain.<br>
<br>
=C2=A03. When error in ulps is higher in some subdomain, it's because t=
here<br>
=C2=A0 =C2=A0 are no reasonably fast algorithms for computing the function<=
br>
=C2=A0 =C2=A0 accurately in that subdomain.<br>
<br>
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 tru=
e answer and the approximation.<br>
<br>
A function with error <=3D 0.5 ulps is *correct*. I call a function with=
<=3D 2 ulps friggin' fantastic, and one with <=3D 10 accurate.<b=
r>
<br>
<blockquote class=3D"gmail_quote" style=3D"margin:0 0 0 .8ex;border-left:1p=
x #ccc solid;padding-left:1ex"><div class=3D"im">
I added addition poly coefficients up to 12 but the impl currently only<br>
uses 7.<br>
<br></div>
impl follows: =C2=A0[...]<br>
</blockquote>
<br>
Paste this at the bottom to test it:<br>
<br>
<br>
(require plot/typed<br>
=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0math)<br>
<br>
(: digamma* (Flonum -> Real))<br>
(define (digamma* x)<br>
=C2=A0 (parameterize ([bf-precision =C2=A053])<br>
=C2=A0 =C2=A0 (bigfloat->real (bfpsi0 (bf x)))))<br>
<br>
(: plot-error ((Flonum -> Flonum) (Flonum -> Real) Real Real -> An=
y))<br>
(define (plot-error f-appx f-truth x-min x-max)<br>
=C2=A0 (plot (function<br>
=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0(=CE=BB: ([x : Real])<br>
=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0(let ([x =C2=A0(fl x)])<br>
=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0(define err (flulp-error (f=
-appx x) (f-truth x)))<br>
=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0(if (rational? err) err -1.=
0)))<br>
=C2=A0 =C2=A0 =C2=A0 =C2=A0 =C2=A0x-min x-max)))<br>
<br>
(plot-error flpsi0 digamma* 0 32)<br>
(plot-error digamma digamma* 0 32)<br>
<br>
(plot-error flpsi0 digamma* -8 0)<br>
(plot-error digamma digamma* -8 0)<br>
<br>
<br>
It looks like your `digamma' is as accurate as `flpsi0' for x > =
21 or so (i.e. <=3D 2 ulps, which is friggin' fantastic). If using 1=
2 terms instead of 7 in your polynomial approximation reduced the error for=
smaller x, it would be worth using that polynomial approximation in the ma=
th library. (I can't test this because I don't know the coefficient=
s' signs.) In the subdomain it would likely be accurate in, `flpsi0'=
; uses 12- to 15-degree Chebyshev polynomials, which are about twice as slo=
w to evaluate as similar-degree polynomials using Horner's rule (as in =
your `digamma-poly').<br>
<br>
Neil =E2=8A=A5<br>
<br>
____________________<br>
=C2=A0Racket Users list:<br>
=C2=A0<a href=3D"http://lists.racket-lang.org/users" target=3D"_blank">http=
://lists.racket-lang.org/<u></u>users</a><br>
</blockquote></div><br></div>
--047d7b3a9684d43f1e04dc27158c--
![[logo]](/logo.png){kind=logo}