MAINT: Use expm1(x) instead of exp(x) - 1

[gaurav1086]

np.exp(1e-99) - 1
0.0
np.expm1(1e-99)
1e-99

[eric-wieser]

@bashtage, can you comment on stream compatibility here?

[bashtage]

Good idea in distributions but legacy should not change unless something isn't compiling. Doesn't really affect the stream except that in very rare cases one might see a different number. I'm not sure that standard exponential would ever produce a value where the precision difference would appear since there are only 2**53 ish distinct values possible.

[bashtage]

Are there numerical difference in reasonable sample sizes? As in any different values generated?

[miccoli]

Doesn't really affect the stream except that in very rare cases one might see a different number. I'm not sure that standard exponential would ever produce a value where the precision difference would appear

I do not fully understand what this means, but in double precision, almost 50% of the samples drawn by pareto(a=1.1) would differ by 1 ULP or more, with respect to the old implementation.

[bashtage]

I do not fully understand what this means, but in double precision, almost 50% of the samples drawn by pareto(a=1.1) would differ by 1 ULP or more, with respect to the old implementation.

I didn't have any idea for which values of expm1(x) and exp(x)-1 differ. The example at the top (1e-99)doesn't seem useful for random pareto, although It could matter for large values of a. I did a quick check and for a=1 the difference is there and so this seems like a good idea to me in distributions (but not for legacy_distributions).

[seberg]

@gaurav1086 would it be possible to add a test that would notice the loss of numerical precision if we changed things back?

[gaurav1086]

@seberg , sure, I think that would be a good idea.

[seberg]

@gaurav1086 could you add the small test? Would be nice to finish this up, since it is an obvious improvement. Let me know if you need pointers for doing that.

[gaurav1086]

@seberg , sorry for the delay. Will add it today. Thank you.

[gaurav1086]

@seberg , added the test in numpy/numpy/random/tests/test_random.py

def test_pareto_expm1(self):
assert_(np.random.pareto(1e99) > 0.0)

Currently,

np.random.pareto(1e99)
0.0

After the change, it should be non-zero of very small magnitude. Please correct me if I am wrong. Thank you.

[seberg]

@gaurav1086 thanks, the test actually fails on master. On the up-side, it is supposed to fail, since np.random.pareto is the legacy implementation, which you hopefully did not change. You should test rng = np.random.default_rng() (or whatever is commonly used in the tests, it may be a different test file).
Also make sure that the test will never fail, so if there is a chance that 1e99 actually does return 0, please set a seed.

[gaurav1086]

@seberg , I changed random_pareto() in distribution.c . The legacy was another change which I reverted back.

[seberg]

@gaurav1086 yes, but your test tests the legacy one, so it fails. Can you fix the test?

[mattip]

You should try running tests locally before pushing. Your fix does not work:
SeedSequence expects int or sequence of ints for entropy not 1e+99
You can run tests locally with:

python3 -mvenv /tmp/py_to_test
/tmp/py_to_test/bin/python -mpip install -r test_requirements.txt
/tmp/py_to_test/bin/python runtests.py