How to write complicated Deterministic?

Questions v5
1

I’m trying to translate examples in “Bayesian methods for hackers” into pyMC v5. Actually it is really painful since many documents from google are mixed with v2, v3, v4. Even I have read the migration guide of v4, there’s still a lot I don’t understand.

To be specific, in the SMS frequency example, how can I construct lambda_ that involves array index manipulation?

lambda_1 = pm.Exponential("lambda_1", alpha)
lambda_2 = pm.Exponential("lambda_2", alpha)
tau = pm.DiscreteUniform("tau", lower=0, upper=n_count_data)

@pm.deterministic
def lambda_(tau=tau, lambda_1=lambda_1, lambda_2=lambda_2):
    out = np.zeros(n_count_data)
    out[:tau] = lambda_1  # lambda before tau is lambda1
    out[tau:] = lambda_2  # lambda after tau is lambda2
    return out

observation = pm.Poisson("obs", lambda_, value=count_data, observed=True)
model = pm.Model([observation, lambda_1, lambda_2, tau])

I tried the following:

def mixed_lambda(lambda1=lambda1, lambda2=lambda2, tau=tau):
    out = np.zeros(len(real_pv))
    t = tau.eval()
    out[:t] = lambda1.eval()
    out[t:] = lambda2.eval()
    return out

with model:
    mixed_lambda = pm.Deterministic("mixed_lambda", mixed_lambda())

However pyMC blames incompatible between np.array and aesara.tensor. Could you tell me the right way to do it? Any related tip and best practice are welcome. Thank you in advance.

2

The book has an official PyMC3 version: GitHub - CamDavidsonPilon/Probabilistic-Programming-and-Bayesian-Methods-for-Hackers: aka "Bayesian Methods for Hackers": An introduction to Bayesian methods + probabilistic programming with a computation/understanding-first, mathematics-second point of view. All in pure Python ;)

That API should be much more similar to V4/V5.

In general you can write complex deterministics just like you would write Numpy (replacing np by pt or at), except for some special operations like branching, looping and writing to arrays (as you are doing with indexing). Writing to arrays is supported via set_subtensor, but it’s probably easier to use where which works just like in numpy.

You should not use eval. That’s for debugging only, as it tries to evaluate the expression you have. You should always give PyMC an unevaluated expression.

import pytensor.tensor as pt

with pm.Model() as m:
  lambda_1 = pm.Exponential("lambda_1", alpha)
  lambda_2 = pm.Exponential("lambda_2", alpha)
  tau = pm.DiscreteUniform("tau", lower=0, upper=n_count_data)

  lambda_ = pt.where(tau <= pt.arange(len(count_data)), lambda_1, lambda_2)
  observation = pm.Poisson("obs", lambda_, observed=count_data)
1 Like
3

Thank you for the reply, it works! pt.where and .eval concepts helps me a lot.

However, when I change pytensor to aesara, it blames at tau <= at.arange(len(real_pv)) part for:

NotImplementedError: Cannot convert ARange{dtype='int64'}.0 to a tensor variable.

File ~/pyenv/lib/python3.10/site-packages/pytensor/tensor/var.py:45, in _tensor_py_operators.__le__(self, other)
     44 def __le__(self, other):
---> 45     rval = at.math.le(self, other)
     46     rval._is_nonzero = False
     47     return rval

File ~/pyenv/lib/python3.10/site-packages/pytensor/graph/op.py:296, in Op.__call__(self, *inputs, **kwargs)
    254 r"""Construct an `Apply` node using :meth:`Op.make_node` and return its outputs.
    255 
    256 This method is just a wrapper around :meth:`Op.make_node`.
   (...)
    293 
    294 """
    295 return_list = kwargs.pop("return_list", False)
--> 296 node = self.make_node(*inputs, **kwargs)
    298 if config.compute_test_value != "off":
    299     compute_test_value(node)
...

It seems the tensor in pytensor and aesara are slightly different. The backend of pyMC changes a lot and really complicated. Why does pyMC mix two or more tensor structure? Is there any preference or idiom about tensor types for daily usage?

4

You cannot mix Aesara and Pytensor!!!

PyMC 4.x uses Aesara
PyMC 5.x uses PyTensor

Otherwise they work the same (PyTensor is a fork of the Aesara project maintained by the PyMC developers)

1 Like