Distribution variable API difficult to use for basic debugging: logp, logcdf, random samples, ppf

[hottwaj]

Description of your problem

This is more of a proposal - please let me know if this would be better put elsewhere.

One thing I often want to do while experimenting with or debugging models is to individually explore the RVs I'm using, but I find this difficult with the current API.

Maybe I'm missing something in the docs, but I find the current API awkward if I then want to: evaluate or plot pdf/cdf/ppf or draw RVs. I find it really helpful to do this especially for distributions where pymc (very usefully) provides parameterisations that are not provided by e.g. scipy (e.g. Beta has mu and sigma parameters). If I can't "see" what pymc is doing with my basic priors then I either have to find the formulae to convert the parameterisation into scipy-compatible parameters, or resort to trial and error (the horror!).

Currently you have to do things like this:

# generate 1000 samples of Bernoulli RV
pymc.Bernoulli.dist(p = 0.04, shape = [1000]).eval()  

# evaluate PDF of Binomial RV over interval 0-1000
pymc.logp(pymc.Binomial.dist(p = 0.04, n=1000), numpy.arange(0, 1000)).eval()  

# I'm not sure of any way to evaluate PPF but maybe I'm missing it?

Is e.g. pymc.Continuous.dist(...) wrapping basic scipy distributions somewhere (I guess not), and if so could those be exposed in a more obvious way so that we can access the scipy API e.g. https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.rv_continuous.html

Alternatively, could a more straightforward API for rvs/pdf/cdf/ppf etc be provided?

See also: #5032

Thanks for the awesome library :)

[hottwaj]

#5308 is also somewhat relevant

[ricardoV94]

Two short notes:

1. We don't have ppfs implemented (yet) for our distributions. We do have cdf via pymc.cdf(dist, value) in case you care about that.
2. You should you pymc.draw(dist, draws) instead of dist.eval() for proper seeding in next releases.

---

For the main question. What do you imagine a more straightforward API would be?

Keep in mind that PyMC defines everything in terms of Aesara logp/random graphs, that's why you always get back a symbolic expression that needs to be compiled before evaluation (what your .eval() is doing).

[hottwaj]

Thanks for your reply! Noted on pymc.draw, thanks

Did you mean pymc.logcdf instead of pymc.cdf? Latter is not available in the version I'm using (4.0.0b6).

As for a more straightforward API I suppose I am biased having had some experience of scipy, but I find that API fairly intuitive
e.g.

my_beta = scipy.stats.beta(a = 1, b = 2)

# Can then call things like my_beta.rvs, or pdf/cdf/ppf/moment/mean etc

My opinion (please feel free to disagree :) ) is that it is easier for someone to "discover" the rvs/pdf/cdf etc methods in this case because the distribution object itself provides them. If you are experimenting in jupyter lab or some IDE, they will come up as suggestions as soon as you write my_rv.[ completion keyboard shortcut ]

With pymc discovery is a bit harder because relevant methods are of form pymc.logcdf, pymc.logp, pymc.draw. It's not obvious to me (as a user of scipy and earlier versions of pymc) that the pymc module is the place to look to find those methods.

Thanks!

[ricardoV94]

Yes I meant logcdf. I understand what you mean now.

It was a decision choice to use functions instead of methods, because we are using vanilla Aesara objects under the hood that do not have such methods.

One longer term reason, is that we are doing some work that will bring logp methods not only for pure Distributions but also graphs based on those, so that you can do:

import pymc as pm
import aesara.tensor as at

censored_normal = at.clip(pm.Normal.dist(), -1, 1)
pm.logp(censored_normal, value).eval()   # It doesn't work just yet, but will hopefully in a future release

The point being that distributions are no longer the only things that have logp (and perhaps one day logcdf) methods. This is already true for draws.

censored_rvs = pm.draw(censored_normal, draws=5)

As such, any work towards adding methods that work with arbitrary Aesara objects and not just the RandomVariables that are returned directly by vanilla PyMC distributions has larger benefits.

We definitely want to make the basic functions logp, logcdf, draw easy to find. Is there any chance that the uneasiness you feel comes from experience with the older versions of PyMC3?

[hottwaj]

Yes you're right I'm sure it makes a difference that I'm more familiar with a previous version.

Thanks for the extra detail and I understand your point about using functions rather than methods. I feel like (again just my opinion) a python API that has lots of module.function(obj) methods for dealing with heterogenous objects is a bit more unusual than an API that has obj.method() methods for dealing with those objects. I still feel like it is harder to know that the methods are "there", and it is not obvious which methods are intended to work with which objects (whereas with obj.method() either the method is there or not, so you know whether it should work on that basis). The more module.function(obj) methods you have, the more time a user will need to sift through a list of functions to find the right one to do the desired thing with some object.

OTOH I guess the team chose the current approach to avoid having to write thin wrappers around aesara objects and all their possible operations? As a user rather than a developer of pymc, I obviously am not aware of how complicated things might be behind the scenes, so I understand if this discussion is not helpful and has already been thought through in much more detail elsewhere :)

Thanks!

[ricardoV94]

The discussion is definitely useful, thanks for bringing it up. Yet the biggest issue is with wrapping everything while not loosing the ability to manipulate those objects freely.

[hottwaj]

One alternative option might be to group the "distribution related" functions into their own module e.g. pymc.dist? That way it might be more obvious what the full list of those methods is (in contrast if everything is in pymc module it takes longer for a user to enumerate and get their head around what is there), and that they all work in a similar way i.e. on some "distribution like" object

Thanks!