Rewrite core of PyMC3 (upcoming v4) to rely on aesara (among many other improvements)

[twiecki]

No description provided.

[ricardoV94]

In case this helps:

1. Since a80cf9a, prior sampling returns the same value for every draw.

with pm.Model() as m:
    x = pm.Normal('x', 0, 1)
    prior = pm.sample_prior_predictive(10)
    print(prior)
    
{'x': array([0.019066, 0.019066, 0.019066, 0.019066, 0.019066, 0.019066,
       0.019066, 0.019066, 0.019066, 0.019066])}

2. For the HalfCauchy (and maybe other distributions) there may be a disconnect between the logp and the random arguments. Our logp considers only the scale argument, but the random method (from scipy) takes loc as well and that's the first optional argument. I don't think there is yet any logic to connect the arguments between the the random and logp, right?

[brandonwillard]

1. Since a80cf9a, prior sampling returns the same value for every draw.

Yes, that's because it's using the same RNG seed every time. The reason why: RandomVariable.inplace == False. This code was supposed to take care of that.

[brandonwillard]

1. For the HalfCauchy (and maybe other distributions) there may be a disconnect between the logp and the random arguments. Our logp considers only the scale argument, but the random method (from scipy) takes loc as well and that's the first optional argument. I don't think there is yet any logic to connect the arguments between the the random and logp, right?

Is there an example of this, or a test that's failing?

[ricardoV94]

1. For the HalfCauchy (and maybe other distributions) there may be a disconnect between the logp and the random arguments. Our logp considers only the scale argument, but the random method (from scipy) takes loc as well and that's the first optional argument. I don't think there is yet any logic to connect the arguments between the the random and logp, right?

Is there an example of this, or a test that's failing?

The random tests are still disabled, so this wouldn't show up right? I wanted to check with the prior predictive sampling but right now that's also unusable.

There was an issue in the logp of the HalfCauchy that lead to the addition of an unused parameter (not sure if related):
https://github.com/pymc-devs/pymc3/pull/4508/files/3db730e1b632e0fa023cd181c60d19de4aab6e48#r594550615

[ricardoV94]

Managed to do it with .dist

import pymc3 as pm
import scipy.stats as st

pymc_samples = pm.HalfCauchy.dist(beta=.01, size=10_000).eval()
print(pymc_samples.mean(), pymc_samples.std())
# (6.647935041418325, 81.64505597965584)

# Shoud match this 
scipy_samples = st.halfcauchy(loc=0, scale=.01).rvs(10_000)
print(scipy_samples.mean(), scipy_samples.std())
# (0.049328030173446065, 0.3650001369075893)

# But matches this
scipy_samples = st.halfcauchy(loc=0.01, scale=1).rvs(10_000)
print(scipy_samples.mean(), scipy_samples.std())
# (5.537612320455185, 73.96359661595363)

# V3
pymc_samples = pm.HalfCauchy.dist(beta=.01, shape=10_000).random()      
print(pymc_samples.mean(), pymc_samples.std())                          
# (0.05829908302484347, 1.122931711761449)

The HalfCauchy is notoriously hard to evaluate, but with these parameters it seems reliable.

[brandonwillard]

The random tests are still disabled, so this wouldn't show up right?

The tests in pymc3.tests.test_distributions perform some sampling and those are running for all the converted Distributions; otherwise, I believe you're talking about pymc3.tests.test_distributions_random. Those tests are not enabled because they're redundant (for the RandomVariables that already exist in Aesara, at least).

I wanted to check with the prior predictive sampling but right now that's also unusable.

It's not unusable; you only need to set/change the seed between calls to pm.sample_*_predictive. Regardless, I've pushed a commit that correctly sets RandomVariable.inplace = True, so you shouldn't need to do that anymore.

[ricardoV94]

Thanks Brandon.

I didn't mean it's unusable. I just didn't know how to make use of it :b

I confirmed my suspicion about the Half Cauchy though. There is a mismatch between the old logp and the new randomOp, as the later expects two arguments, while the former expects only one.

That's why the logp test would fail if the unused alpha argument was not added here: https://github.com/pymc-devs/pymc3/blob/a80cf9ac7dec48cc71a2924362884a4d8cc6e30d/pymc3/distributions/continuous.py#L2276

It's not too difficult to add this loc parameter to the logp, but just wanted to make sure we want to do that. Do we want to adapt to the number of arguments expected in aesara or just customize our own random ops in terms of the arguments we want to use?

[brandonwillard]

It's not too difficult to add this loc parameter to the logp, but just wanted to make sure we want to do that. Do we want to adapt to the number of arguments expected in aesara or just customize our own random ops in terms of the arguments we want to use?

If an adaptation is something as simple as inverting a parameter, then we'll always want to do that, but this case sounds like we'll need to do the latter (i.e. customize a RandomVariable). There's now an example of such a customization for Multinomial.

[ricardoV94]

On another front:

import pymc3 as pm
with pm.Model() as m:
    x = pm.Uniform('x', 0, 1)
    y = pm.Uniform('y', 0, x)
    
m.logp({'x_interval': 0, 'y_interval': 0})
# array(-2.83990034)
m.logp({'x_interval': 0, 'y_interval': 0})
# array(-2.09721045)
m.logp({'x_interval': 0, 'y_interval': 0})
# array(-2.23865389)

Also it ignores when I set transform=None, and creates the interval variables anyway.

[michaelosthege]

@ricardoV94 were the trailing __ on the transformed variables removed?

On master:

>>> m.logp({'x_interval': 0, 'y_interval': 0})
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "C:\Users\osthege\Repos\pymc3-dev\pymc3\model.py", line 1565, in __call__
    return self.f(**point)
  File "C:\Users\osthege\AppData\Local\Continuum\miniconda3\envs\pm3-dev\lib\site-packages\aesara\compile\function\types.py", line 956, in __call__
    raise TypeError(
TypeError: Missing required input: x_interval__ ~ TransformedDistribution
>>> m.logp({'x_interval__': 0, 'y_interval__': 0})
array(-2.77258872)
>>> m.logp({'x_interval__': 0, 'y_interval__': 0})
array(-2.77258872)
>>> m.logp({'x_interval__': 0, 'y_interval__': 0})
array(-2.77258872)

[ricardoV94]

Yes they were removed for the time being.

[brandonwillard]

Also it ignores when I set transform=None, and creates the interval variables anyway.

I'll add fixes for those in a minute.

[michaelosthege]

Made these comments a while back, but forgot to hit submit..