oct issueshttps://gitlab.common-lisp.net/oct/oct/-/issues2020-12-28T06:15:02Zhttps://gitlab.common-lisp.net/oct/oct/-/issues/21Fix incomplete-gamma test failures2020-12-28T06:15:02ZRaymond ToyFix incomplete-gamma test failuresWith the tests added for #18, we need to fix the tests that are failing with errors. At least some of the tests are failing because the continued fraction isn't converging. This probably means we need to add a series implementation as ...With the tests added for #18, we need to fix the tests that are failing with errors. At least some of the tests are failing because the continued fraction isn't converging. This probably means we need to add a series implementation as well.Raymond ToyRaymond Toyhttps://gitlab.common-lisp.net/oct/oct/-/issues/15Error rates for log-gamma2021-01-28T05:58:28ZSteveError rates for log-gammaI've ported the test data from Boost for [log-gamma](https://www.boost.org/doc/libs/1_74_0/libs/math/doc/html/math_toolkit/sf_gamma/lgamma.html) and [gamma](https://www.boost.org/doc/libs/1_74_0/libs/math/doc/html/math_toolkit/sf_gamma/t...I've ported the test data from Boost for [log-gamma](https://www.boost.org/doc/libs/1_74_0/libs/math/doc/html/math_toolkit/sf_gamma/lgamma.html) and [gamma](https://www.boost.org/doc/libs/1_74_0/libs/math/doc/html/math_toolkit/sf_gamma/tgamma.html) and thought the results might be helpful. It seems that in some cases, OCT does better than GSL and Rmath, on other intervals worse, sometimes much worse. See, for example near 0, -10 and -55.
| Test | OCT | Boost |
| ------ | ------ | ------ |
| Factorial | Max = 2.10e+1ε<br/>Mean = 1.32e+0ε | Boost: Max = 0ε (Mean = 0ε)<br/>(GSL 2.1: Max = 33.6ε (Mean = 2.78ε))<br/>(Rmath 3.2.3: Max = 1.55ε (Mean = 0.592ε)) |
| near 0 | Max = 8.35e+15ε<br/>Mean = 2.42e+15ε| Boost: Max = 0ε (Mean = 0ε)<br/>(GSL 2.1: Max = 5.21ε (Mean = 1.57ε))<br/>(Rmath 3.2.3: Max = 0ε (Mean = 0ε)) |
| near 1 | Max = 3.91e+1ε<br/>Mean = 1.45e+1ε| Boost Max = 0ε (Mean = 0ε)<br/>(GSL 2.1: Max = 442ε (Mean = 88.8ε))<br/>(Rmath 3.2.3: Max = 7.99e+04ε (Mean = 1.68e+04ε)) |
| near 2 | Max = 3.91e+1ε<br/>Mean = 1.46e+1ε| Boost: Max = 0ε (Mean = 0ε)<br/>(GSL 2.1: Max = 1.17e+03ε (Mean = 274ε))<br/>(Rmath 3.2.3: Max = 2.63e+05ε (Mean = 5.84e+04ε)) |
| near -10 | Max = 6.91e+15ε<br/>Mean = 1.70e+15ε| Boost: Max = 0ε (Mean = 0ε)<br/>(GSL 2.1: Max = 24.9ε (Mean = 4.6ε))<br/>(Rmath 3.2.3: Max = 4.22ε (Mean = 1.26ε)) |
| near -55 | Max = 1.82e+14ε<br/>Mean = 8.74e+13ε| Boost: Max = 0ε (Mean = 0ε)<br/>(GSL 2.1: Max = 7.02ε (Mean = 1.47ε))<br/>(Rmath 3.2.3: Max = 250ε (Mean = 60.9ε)) |
This is with double-floats. Could this be related to Common Lisp having 62 bits to work with rather than the 64 of the C implementations? If that's true, I'd have expected some manner of difference, but not the orders of magnitude that some of the intervals have.
I've attached the log file from the SBCL test run here, which has all the details.
[gamma-test-run.log](/uploads/cb5ce9002ccd31bb47f8686995230663/gamma-test-run.log)https://gitlab.common-lisp.net/oct/oct/-/issues/12Implement inverse-incomplete-gamma2020-11-24T22:11:43ZSteveImplement inverse-incomplete-gammainverse-incomplete-gamma is used for computing quantiles of the gamma distribution and is a fundamental computation in a statistical setting. This issue proposes to develop this, and possibly other missing gamma related functions, to 'ro...inverse-incomplete-gamma is used for computing quantiles of the gamma distribution and is a fundamental computation in a statistical setting. This issue proposes to develop this, and possibly other missing gamma related functions, to 'round out' the gamma family.https://gitlab.common-lisp.net/oct/oct/-/issues/7OSX and linux test suite has different results.2020-11-22T16:09:54ZRaymond ToyOSX and linux test suite has different results.The test `BESSEL-J-1/2.Q.1` had a accuracy threshold of 169.45 bits. But when running on osx, the test passes. But on linux it fails with an actual accuracy of 164.9.
Don't understand why there's a difference, but it should be investi...The test `BESSEL-J-1/2.Q.1` had a accuracy threshold of 169.45 bits. But when running on osx, the test passes. But on linux it fails with an actual accuracy of 164.9.
Don't understand why there's a difference, but it should be investigated.
The output from linux for this test is:
```
Actual value:
(((7
#q1.0050537292583339756757221701852075697640307680435648743808955094q0)
(164.9639202995993d0 169.45
#q(#q0.67187119306031800410962218259898937011015226586845701354319806663q0
#q-9.1022800468382580449905590670828511598385571251839047903692211539q-67)
#q0.67187119306031800410962218259898937011015226586844228371069762899q0)))
```Raymond ToyRaymond Toyhttps://gitlab.common-lisp.net/oct/oct/-/issues/4Implement incomplete-beta family of functions2020-11-28T03:09:14ZSteveImplement incomplete-beta family of functionsThe 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))))
```...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](https://github.com/codeplea/incbeta) and [an implementation in the CRISP genomic project](https://github.com/vibansal/crisp/blob/e6f296df58171d86e0bd85ee645656b1bcacb595/FET/kfunc.c#L118). 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)