Changes to logprob computation between v4 and v5

[michaelosthege]

Under v4 this was perfectly fine to evaluate a likelihood:

x_true = np.array([1, 2, 3, 4, 5])
y_obs = np.random.normal(2 * x_true)

with pm.Model() as pmodel:
    x_hat = pm.Uniform("x_hat", 0, 1, shape=x_true.shape, transform=None)
    # L = pm.StudentT("L", mu=x_hat, observed=y_obs)  # the first version of the code sample was this,
                                                        but this was not actually what I ran.
    rv = pm.StudentT("L", mu=x_hat, observed=y_obs)
    L = pm.logprob.joint_logprob({rv: y_obs}, sum=True)  # This is what I actually had
assert isinstance(L, at.TensorVariable)
assert L.ndim == 0

# PyMC v4 returns the RV, but for .eval() we need the RV-value-variable
vvar = pmodel.rvs_to_values[x_hat]

# compare the loglikelihood from PyMC with the one from SciPy
x_test = numpy.random.normal(x_true, scale=0.1)
actual = L.eval({vvar: x_test})
expected = scipy.stats.t.logpdf(loc=x_test, loc=y_obs, scale=1)
numpy.testing.assert_almost_equal(actual, expected, 6)

But with v5 I get an unused input error:

pytensor.compile.function.types.UnusedInputError: pytensor.function was asked to create a function computing outputs given certain inputs, but the provided input variable at index 0 is not part of the computational graph needed to compute the outputs: x_hat.

And compiling a function with on_unused_input="ignore" gives different results:

f = pytensor.function(inputs=[vvar], outputs=[L], on_unused_input="ignore")
actual = f(x_test)[0]
numpy.testing.assert_almost_equal(actual, expected, 6)

E       AssertionError: 
E       Arrays are not almost equal to 6 decimals
E       
E       Mismatched elements: 1 / 1 (100%)
E       Max absolute difference: 14.72333561
E       Max relative difference: 2.29985815
E        x: array(-21.125181)
E        y: array(-6.401845)

I'm asking this because it's currently blocking my PR to upgrade our calibr8 package for PyMC v5 compatibility: https://github.com/JuBiotech/calibr8/pull/61/commits

[ricardoV94]

L is a RV not a logp expression. Theres no place where the value variable enters the graph of the RV, hence the correct unused input error.

If you want the logp of L call model.logp(vars=[L], sum=False) (or compile_logp if you want a compiled function)

[michaelosthege]

I don't want the logp which already includes the prior's contribution! I want only the log-likelihood so I can check if it's identical to the one I get from SciPy.

Both model.logp and compile_logp include the prior or am I misunderstanding something this fundamental?

[ricardoV94]

Yes, you only get the terms of the variables in the vars list. If you say vars=[L], you only get the likelihood.

Also this was exactly the same in V4. You couldn't get the likelihood with your original snippet then either.

[ricardoV94]

You can check this is what model.observedlogp does:

pymc/pymc/model.py

Lines 879 to 882 in 23a1636

	@property
	def observedlogp(self) -> Variable:
	"""PyTensor scalar of log-probability of the observed variables"""
	return self.logp(vars=self.observed_RVs)

[michaelosthege]

So first off, a massive thank you, for helping me out here!

Your answer is actually correct here, because the code snippet I posted was actually incorrect/incomplete.
I updated it.

Essentially the change boiled down to the following:

Beforehand I used

# if in modelcontext
    rv = pm.Distribution(..., observed=y_obs)
# else
    rv = pm.Distribution.dist(...)
pm.joint_logp(rv, y_obs, sum=True)

and this 👆 worked in v4.

For v5 I attempted this 👇 which does NOT work.

# if in modelcontext
    rv = pm.Distribution(..., observed=y_obs)
# else
    rv = pm.Distribution.dist(...)
pm.logprob.joint_logprob({rv, y_obs}, sum=True)

Using different ways to construct the logplikelihood graph is what fixed the issue:

# if in modelcontext
    rv = pm.Distribution(..., observed=y_obs)
    pmodel.logp([rv], sum=True)
# else
    rv = pm.Distribution.dist(...)
    pm.logprob.joint_logprob({rv, y_obs}, sum=True)

If someone is interested in the complete picture, the changes are part of this PR: JuBiotech/calibr8#61