Speedup logcdf tests

[ricardoV94]

The logcdf tests were running very slow due to constantly rebuilding the logcdf function.

Also fixed a couple of failing tests on float32 (which were difficult to identify before via because setting n_samples=-1 took ages to run.

[ricardoV94]

I am getting a strange issue with some of the logcdf methods which seems to be behind the failing tests:

Setting the values to constants gives the correct values

with pm.Model() as model:
    value = pm.NegativeBinomial('value', mu=0.9, alpha=0.1)

logcdf = logpt(model['value'], cdf=True)
fn1 = model.fastfn(logcdf)

print(fn1({'value': 1.0}))  # -0.14408081390569186

But using shared variables does not

with pm.Model() as model:
    mu = aesara.shared(np.asarray(0.9))
    alpha = aesara.shared(np.asarray(0.1))
    value = pm.NegativeBinomial('value', mu=mu, alpha=alpha)

logcdf = logpt(model['value'], cdf=True)
fn2 = model.fastfn(logcdf)

mu.set_value(0.9)
alpha.set_value(0.1)
print(fn2({'value': 1.0}))  # -0.08777554398474162

This doesn't happen with all distributions, and unfortunately those that do have a crazy long graph to be able to compare easily side by side :/

Found the culprint: it's the incomplete_beta function:

incomplete_beta(0.1, 2.0, 0.1).eval()   # array(0.86581778)

alpha = aesara.shared(0.1, 'alpha')
beta = aesara.shared(2.0, 'beta')
incomplete_beta(alpha, beta, 0.1).eval()   # array(0.91596645)

import scipy.special
scipy.special.betainc(0.1, 2.0, 0.1)   # 0.8658177758494668

All failing tests (expect for the HyperGeometric, which is a different issue) rely on the incomplete_beta.

Maybe it's time I finish #4519

[ricardoV94]

Also found some issues with the recent initval changes, as it doesn't respect the parents initvals.

This snippet often leads to a ValueError:

with pm.Model() as model:
    mu = pm.Normal('mu', initval=100)
    alpha = pm.HalfNormal('alpha', initval=100, transform=None)
    value = pm.NegativeBinomial('value', mu=mu, alpha=alpha)

model.initial_values
# {mu: array(100., dtype=float32),
#  alpha: array(100., dtype=float32),
#  value: array(1)}

When it doesn't fail, initval is still far from the expected ~ 100:

pm.NegativeBinomial.dist(mu=100, alpha=100).eval()
# array(98)

[ricardoV94]

The HyperGeometric test is failing because of this issue in Aesara: pymc-devs/pytensor#461

[ricardoV94]

Closing in favor of #4736