Discretization of Continuous Variables in Pymc

You need to replace np. operations with pt. operations, and remove .eval() calls. For example:

    plf_bin = np.digitize(PLF.eval(), bins) - 1

is equivalent to:

import pytensor.tensor as pt
plf_bin = pt.searchsorted(bins, PLF) - 1

(This equivalence between digitize and searchsorted is described in the documentation for np.digitize)

You will also not be able to use a dictionary for lookup_table_COM, because current_state will not be a tuple of integers, it will be a tuple of symbolic computations. Instead, you can make it an array, and you can make it an array:

lookup_table_COM = np.zeros((num_states_plf), num_states_slf))
for i, j in itertools.product(range(num_states_plf), range(num_states_slf)):
    idx = i * num_states_slf + j
    lookup_table_COM[i, j] = Combination_cpt[idx]

lookup_table_pt =pt.as_tensor_variable(lookup_table_COM)

which you can then symbolically index into just like a numpy array:

p_combination = lookup_table_pt[plf_bin, slf_bin]