Groups | Search | Server Info | Keyboard shortcuts | Login | Register [http] [https] [nntp] [nntps]

Groups > comp.lang.forth > #23328

GammaLn, Complete Beta function, Probability Density function

Show all headers | View raw

Hi!

I designed the following functions, partly not available in FSL:
- GammaLn
- Complete Beta function
- Beta Probability Density function

AFAIK the GammaLn function is also in the Gamma code of FSL, but this one uses a Lanczos approximation. It seems to be a bit slower, but a little more precise close to zero (about 13 digits of precision). The Beta function and PDF function build on that, but they also work fine with the loggamma of FSL.

Excuse the sometimes odd code, but this is a straight port from 4tH.

Hans Bezemer

---8<---
create _L 6 floats allot

76.18009172947146e    _L 0 floats + f!
-86.50532032941677e   _L 1 floats + f!
24.01409824083091e    _L 2 floats + f!
-1.231739572450155e   _L 3 floats + f!
0.1208650973866179e-2 _L 4 floats + f!
-0.5395239384953e-5   _L 5 floats + f!

: gammaln                              ( f1 -- f2)
  190015e-15 1e f+ fover
  6 0 do 1e f+ _L i floats + f@ fover f/ frot f+ fswap loop fdrop
  fover fdup 55e-1 f+ fdup fln
  frot 1e 2e f/ f+ f* f- fnegate fswap
  746310005e-16 25066282e-7 f+ f* frot f/ fln f+
;

fvariable _beta                        \ temporary variable
                                       \ be nice to the tiny FP stack ;-)
: fbeta                                ( f1 f2 -- f3) 
  fover fover f+ _beta f! gammaln fswap gammaln f+ _beta f@ gammaln f- fexp
;

fvariable _arg1                        \ first argument
fvariable _arg2                        \ second argument
                                       \ be nice to the tiny FP stack ;-)
: fbeta.pdf                            ( f1 f2 f3 -- f4)
  fover _arg1 f! fdup _arg2 f!         ( x a b)
  1e f- frot 1e fover f-               ( a b-1 x 1-x)
  frot f** frot frot                   ( 1-x**b-1 a x)
  fswap 1e f- f** f*                   ( 1-x**b-1*x**a-1)
  _arg1 f@ _arg2 f@ fbeta f/
;
---8<---

Back to comp.lang.forth | Previous | Next | Find similar

Thread

GammaLn, Complete Beta function, Probability Density function The Beez <the.beez.speaks@gmail.com> - 2013-06-08 09:36 -0700

csiph-web