Establish new test framework for the "random" behavior of Distributions

[ricardoV94]

#4508 and #4548 started the refactoring of Distributions whose random Ops were already defined in Aesara, and we should be able to start converting the remaining ones anytime now.

It seems that the logp/ logcdf tests in test_distributions.py should stay more or less unchanged, but we might want to change or remove most of the tests in test_distributions_random.py.

Right now, in test_distribtions_random.py, we are generating a couple hundred samples via pymc3 and scipy and checking they match by doing a Kolmogorov-Smirnov test.

@brandonwillard said that since we are basically just calling numpy / scipy random routines under the hood, it should be enough to check we are converting properly between our parametrizations and the scipy ones.

For the shapes, we can assume that if things are working on Aesara side, they should be working here?

---

This is the minimalist example of what we need to check:

exp = pm.Exponential.dist(lam=10)
exp.eval()  # Make sure we are calling the scipy.stats.expon(1/lam) under the hood

---

Anything I am missing?

[junpenglao]

+1 to this, we should set up using https://hypothesis.readthedocs.io/en/latest/

See for example https://github.com/tensorflow/probability/blob/master/tensorflow_probability/python/internal/hypothesis_testlib.py

[matteo-pallini]

If it's still free, I would be keen to work on this

[ricardoV94]

It is still free. Do you have any ideas on how to tackle this?

[matteo-pallini]

My understanding is that given that random in pymc is being deprecated (aesara is currently handling that) all the tests in test_distribtions_random using pymc3_random and pymc3_random_discrete can be deleted (ad per comment here). Theoretically the ones for which the refactoring has been done in PyMC.

> we are basically just calling numpy / scipy random routines under the hood

I trust this, but I haven't understood what's based on (maybe a discussion somewhere else?). It would be great to get more context.

This is the minimalist example of what we need to check:

Still not sure how to do this , but to be fair I haven't spent too much time on it, and haven't read the implementation on the aesara side. I should get a better idea after having done so (tomorrow/day after). Still, If you already have a good approach in mind, please share it.

[matteo-pallini]

I am still having a hard time understanding where exactly in aesera the RandomVariable distribution to be used is being specified. My current guess is that it relies on the class name (similarly to rng_fn in RandomVariable).

Despite that, I would expect the two samples below to converge, given a large enough n.

import numpy as np
from pymc3.distributions.continuous import Exponential

exp_np = np.random.RandomState().exponential
exp_vals_np = [exp_np(scale=10) for _ in range(10_000)]

exp_pymc = Exponential.dist(lam=10)
exp_vals_pymc = [exp_pymc.eval().tolist() for _ in range(10_000)]

So, possibly the Kolmogorov–Smirnov and the Chi-Square currently ran respectively by pymc3_random_discrete and pymc3_random could be simply ran on samples generated through the above? I think that this change would effectively move the existing tests to perform their original task on the refactored version of the distributions.

However, given that

brandonwillard said that since we are basically just calling numpy / scipy random routines under the hood, it should be enough to check we are converting properly between our parametrizations and the scipy ones.

This approach may be overkill, I still haven't considered how to test the right parameterization. I have the feeling that one way or another it will still be necessary to rely on the distribution generated...

Let me know what you think

[brandonwillard]

I am still having a hard time understanding where exactly in aesera the RandomVariable distribution to be used is being specified. My current guess is that it relies on the class name (similarly to rng_fn in RandomVariable).

This is the default RandomVariable.rng_fn; it uses self.name.

[ricardoV94]

The exponential is not working right now (#4576) that's why you see those failing tests (if you try 1/10 in one or the other they should match). Trying with another distribution like normal that can be given in terms of tau and sigma should give you a better idea of what is needed.

[brandonwillard]

This approach may be overkill, I still haven't considered how to test the right parameterization. I have the feeling that one way or another it will still be necessary to rely on the distribution generated...

This should be straightforward; for example, if the underlying RandomVariable uses a different parameterization than the corresponding Distribution, create an instance of the Distribution and make sure that the resulting RandomVariable's parameters are correct for their parameterization.

Here's an illustration based on your example:

import numpy as np

import aesara

from pymc3.distributions.continuous import Exponential


test_lambda = 10.

rng_at = aesara.shared(np.random.RandomState(23092))
exp_pymc = Exponential.dist(lam=test_lambda, rng=rng_at)

# We only need to confirm that the parameterization from `Exponential` to
# `ExponentialRV` is correct.  We can do this directly, exactly, and without
# resorting to costly statistical tests.
#
# `ExponentialRV` uses the `scale` argument to `np.random.exponential`, 
# so we only need to check that the `scale` parameter is equal to one over 
# the `lam` value we gave to `Exponential.`
#
# All the `Apply` node inputs after the third one are distribution parameters,
# so we get and compare those.
rv_scale, = exp_pymc.owner.inputs[3:]

assert rv_scale.eval() == 1/test_lambda

Otherwise, are the NumPy results and the RandomVariable results the same?

np.random.seed(23092)
exp_vals_np = np.array([np.random.exponential(scale=10) for _ in range(10_000)])
exp_vals_pymc = np.array([exp_pymc.eval() for _ in range(10_000)])

>>> exp_vals_np.mean() - exp_vals_pymc.mean()
0.0

They are.

[ricardoV94]

@brandonwillard I think you are missing the fact that the exponential scale parameter is not the pymc3.Exponential lambda parameter. A lambda of 10 should result in an exponential of 1/10 = 0.1. This confusion is exactly why we need to be careful with these tests.

In [4]: np.random.exponential(scale=10, size=1000).mean()                             
Out[4]: 10.024648140707738

# pymc3 master
In [7]: pm.Exponential.dist(lam=10).random(size=1000).mean()                    
Out[7]: 0.10322241829781785

[brandonwillard]

I think you are missing the fact that the exponential scale parameter is not the pymc3.Exponential lambda parameter.

No, we're talking about how to make tests for bugs like this, and not that there are bugs like this.

FYI: I updated my example above to highlight that the current parameterizations are the same, although they shouldn't be. Instead, I changed the example to look like the exact kind of test we would want for Exponential.

[matteo-pallini]

Thanks to both, in particular for the detailed illustration and pointing out that Exponential is broken.

I created a WIP PR with a first possible approach to test the parameterization. I still need to add more distributions and also check whether there are other tests which is worth removing.

I am still having a hard time understanding where exactly in aesera the RandomVariable distribution to be used is being specified. My current guess is that it relies on the class name (similarly to rng_fn in RandomVariable).

Fundamentally I simply needed to spend more time on the docs. At that point, I only really read the aesara code, but I got lost into the aesara.compile.function.function. So, it still wasn't still clear to me that Op.perform was the method responsible for computing the output variables (in this case, instantiate the right distribution and get the random variable).

[brandonwillard]

So, it still wasn't still clear to me that Op.perform was the method responsible for computing the output variables (in this case, instantiate right distribution and get random variable).

In general, Op.perform isn't necessarily responsible for a computation containing an Op. When function compiles a graph, it may remove the Op during optimization, and, when the Op also has a C implementation it may use that instead of the Python implementation provided by Op.perform.

RandomVariables only have Python implementations right now, so the C implementation thing won't get in the way. Also, I don't believe there are any optimizations that will remove a RandomVariable unless it simply isn't being used. That said, in RandomVariable's case, you can expect that RandomVariable.perform will be called.

[matteo-pallini]

Closed the original PR #4580 after messing up the branch, and opened a new one #4608

I just need to check whether 1 test can be removed and then it should be ready for review (probably tomorrow)

[ricardoV94]

Closed via #4608