Re: [statistics] Release 1.0

[Alex Herbert <al...@gmail.com>]

I have finished working on a revised test suite for [statistics]. It
is ready to merge into master. I will do this when it has passed the
CI checks [1].

The test suite uses properties files to contain reference values for a
parameterisation of a distribution. Ideally all distributions should
be tested using these files as it ensures that the distribution is run
through all the standard tests. In the case of some old legacy tests
for issues discovered in Commons Math the tests have been preserved.

I updated most test data to higher accuracy and now most of the test
suite passes using a *relative* tolerance of 1e-12. Using a relative
tolerance is helpful for assessing the inverse CDF as the output value
x may not lie in the interval [0, 1]. For testing probabilities an
absolute accuracy is OK but if the probabilities have a large range of
magnitudes then a high tolerance needed for some values will not
identify problems in small values. Absolute accuracy is supported by
JUnit assertions, but not relative accuracy. For relative accuracy I
created a wrapper interface (DoubleTolerance) that uses Commons
Numbers Precision to compute equivalence. This can use absolute,
relative or ULP distance tolerances; and combinations of them with And
and Or operations.

I fixed a lot of issues that were missed due to sparsity of the
previous test suite. I found additional cases for distributions which
had bugs in functions that were not tested. I found domain issues with
overflow of integers to negative. There are so many fixes that it is
better to merge this into master and work on specific outstanding
issues.

I will create tickets for the following unresolved issues that may
need to be addressed before a release:

1. The PoissonDistribution can be parameterised with any mean and yet
the PoissonSampler in RNG excludes a mean above 2^30 to avoid
truncation. This change was made to avoid distribution truncation and
because the algorithm uses (int) floor(mean) (see [5]). Commons RNG
now has a LongSampler interface and it may be possible to create a
PossionSampler that outputs a long with the same algorithm and avoid
truncation to int values. The final PoissonDistribution is still
bounded by int values but the sampler can be less restrictive if it
uses long. At large means the Poisson sampler mainly uses a Gaussian
sampler and runtime performance should not be impacted.

2. Many distributions only compute the inverse CDF to 1e-9 relative
error. Some fail to compute when the input p value maps to a very
small x value (e.g. 1e-50). This is related to the defaults for the
brent solver used in the default implementation. I tried to adjust the
defaults but found a bug in the brent solver, which I have fixed [6].
The brent solver requires a release to be used with this fix so work
can use a local copy until all issues with it have been resolved, and
the changes ported back to numbers.

3. The Normal distribution uses InverseErfc for the inverse CDF. This
is a wrapper for InverseErf and as such it is limited in precision as
p -> 0 since the method computes using (2p - 1). When p is < 1e-16
then the inverse CDF returns -infinity. The smallest z is around -8.2
as computed by matlab:

norminv(1e-16)
-8.222082216130437

scipy, R, matlab, octave, Mathematica all have methods to invert small
p. Here is output from matlab for norminv(5e-324):
-38.472346342768788. It is accurate down to the minimum value for a
64-bit float.

There are two free libraries that seem to be used for this: Boost and
Cephes (function ntdri). The Boost licence is compatible with Apache
and there are some Boost derived works already in Commons Numbers.
Incorporating the Boost Erf functions into numbers would be useful.

The inability to invert the Erfc for all p affects the Normal,
Truncated normal and Levy distributions.

4. The Poisson distribution currently allows user configured
tolerances for computing the Gamma functions for the CDF. These may
not be necessary and it should be investigated. The version of the
Gamma function used by the Poisson distribution may not be the same as
those that required these tolerances in Commons Math 2 (use of the
tolerances go as far back as the commit history in the git repo).

5. The Gamma distribution density is not accurate when shape < 1.
There are extensive tests related to MATH-753 [7] for accuracy of the
Gamma distribution. These test shape >= 1. The work requires extension
to smaller shapes where the density -> infinite as X -> 0. When shape
>= 1 then density -> 0 as X -> 0 (see [4]). This results in inaccuracy
of the Chi-squared distribution which uses Gamma(df/2, 2) when the
degrees of freedom (df) are < 2. I've read the code and the
documentation for the Boost functions which the code is based on and I
think the switch to the alternative computation to avoid overflow is
currently not working for all cases.

Regards,

Alex

[1] https://github.com/apache/commons-statistics/pull/31
[2] Boost error function inverses:
https://www.boost.org/doc/libs/1_77_0/libs/math/doc/html/math_toolkit/sf_erf/error_inv.html
[3] Cephes (ntdri in the probability functions) : https://www.netlib.org/cephes/
[4] https://en.wikipedia.org/wiki/Gamma_distribution
[5] https://issues.apache.org/jira/browse/RNG-52
[6] https://issues.apache.org/jira/browse/NUMBERS-168
[7] https://issues.apache.org/jira/browse/MATH-753

---------------------------------------------------------------------
To unsubscribe, e-mail: dev-unsubscribe@commons.apache.org
For additional commands, e-mail: dev-help@commons.apache.org

Re: [statistics] Release 1.0

[Alex Herbert <al...@gmail.com>]

On Wed, 6 Oct 2021 at 23:20, Gilles Sadowski <gi...@gmail.com> wrote:
>
>
> Are there other issues with "BrentSolver"?
> I'd suggest that we release v1.1 (or v1.0.1) of Commons Numbers
> (with the fix for NUMBERS-167 and the changes related to Erfc).

I don't think so. The bug I found was in other brent solver
implementations I compared it to. It is a standard way to look for a
sign change:

double x, y;

if (x * y < 0) {
    // sign change
}

But this does not work if x*y underflows to -0.0. The sign change is
not detected.

I stopped working on the brent solver used in the inverse CDF as I
moved on to other things. I just need a bit more time to look at it.
It may just require some updates to the tolerances that are being
used.

Alex

---------------------------------------------------------------------
To unsubscribe, e-mail: dev-unsubscribe@commons.apache.org
For additional commands, e-mail: dev-help@commons.apache.org

Re: [statistics] Release 1.0

[Gilles Sadowski <gi...@gmail.com>]

Le mer. 6 oct. 2021 à 23:16, Alex Herbert <al...@gmail.com> a écrit :
>
> I have finished working on a revised test suite for [statistics]. It
> is ready to merge into master. I will do this when it has passed the
> CI checks [1].

+1

>
> [...]
> I fixed a lot of issues [...]

Thanks!

>
> I will create tickets for the following unresolved issues that may
> need to be addressed before a release: [...]
>
>
> 2. Many distributions only compute the inverse CDF to 1e-9 relative
> error. Some fail to compute when the input p value maps to a very
> small x value (e.g. 1e-50). This is related to the defaults for the
> brent solver used in the default implementation. I tried to adjust the
> defaults but found a bug in the brent solver, which I have fixed [6].
> The brent solver requires a release to be used with this fix so work
> can use a local copy until all issues with it have been resolved, and
> the changes ported back to numbers.

Are there other issues with "BrentSolver"?
I'd suggest that we release v1.1 (or v1.0.1) of Commons Numbers
(with the fix for NUMBERS-167 and the changes related to Erfc).

Best,
Gilles

---------------------------------------------------------------------
To unsubscribe, e-mail: dev-unsubscribe@commons.apache.org
For additional commands, e-mail: dev-help@commons.apache.org