Implement incomplete-beta family of functions
The beta function is easily implemented in terms of gamma:
(defun beta (p q)
"Compute beta with gamma function"
(/ (* (NET.COMMON-LISP.OCT::gamma p) (NET.COMMON-LISP.OCT::gamma q))
(NET.COMMON-LISP.OCT::gamma (+ p q))))
More useful however are the integral functions of beta. There are two BSD/MIT licensed version based on gamma, incb.c and an implementation in the CRISP genomic project. Both use Lentz's algorithm to evaluate the continued fraction, which OCT also uses for gamma.
It seems it should be possible to add the incomplete beta function in a similar manner. This isn't a type of programming I do often, and it's not entirely clear how to use the OCT lentz function to port the code from one of these implementations. Perhaps if someone could sketch out how a a 'proper' (from an OCT perspective) implementation of upper and lower incomplete beta might work I can run with it.
To use in a statistical context, we need:
- incomplete-beta-lower
- incomplete-beta-upper
- beta
- inverse-incomplete-beta (to compute quantiles)