ValueError: Random variables detected in the logp graph

Questions v5
1

Hi all,

I’m trying to learn the latest version of pytensor/pymc by implementing the Event-Based Model (Event Based Model — Disease Progression Modelling) in it. It’s a straightforward generative process. But I’ve been finding it a bit tricky as it involves estimating a permutation.

First, here’s a summary of the data:

synth_data = pd.read_csv('./data/synthetic_400_3.csv')
synth_data_np = synth_data.to_numpy()
K = synth_data_np.shape[1] -1
N=synth_data_np.shape[0]
synth_data.shape```

(400, 4)```

The data are 400 data points with three features and the last dimension is a binary class label.

Here’s the model definition code (which runs fine).

my_mod = pm.Model()
with my_mod:
    musd = pm.Normal.dist( mu=np.array([0., 0., 0.]), sigma=np.array([1, 1, 1]))
    sigsd = pm.HalfNormal.dist( sigma=np.array([2., 2., 2.]))
    musc = pm.Normal.dist(mu=np.array([-1., -1., -1]), sigma=np.array([1,1,1]))
    sigsc = pm.HalfNormal.dist(sigma=np.array([2.,2.,2.]))
    S = pt.tensor.random.permutation(x=K, name='S')

    for n in np.arange(N):
        n_k = pm.DiscreteUniform(f'n{n}_k', lower=0, upper=K-1)
        cur_mus = pt.tensor.switch(pt.tensor.ge(S,n_k), musd, musc)
        cur_sigs = pt.tensor.switch(pt.tensor.ge(S,n_k), sigsd, sigsc)
        if synth_data_np[n,K] == 0:
            cur_mus = musc
            cur_sigs = sigsc
        x = pm.Normal(f'x_{n}', mu=cur_mus, sigma=cur_sigs, observed=synth_data_np[n,0:K])

Here’s the code for sampling and error I get:

with  my_mod:
    samps = pm.sample()

--------------------------------------------------------------------------- ValueError Traceback (most recent call last) Cell In[15], line 2 1 with my_mod: ----> 2 samps = pm.sample() File [~/miniconda3/envs/pymc_env/lib/python3.11/site-packages/pymc/sampling/mcmc.py:653](https://vscode-remote+wsl-002bubuntu.vscode-resource.vscode-cdn.net/home/jausterw/work/playjuly2023/~/miniconda3/envs/pymc_env/lib/python3.11/site-packages/pymc/sampling/mcmc.py:653), in sample(draws, tune, chains, cores, random_seed, progressbar, step, nuts_sampler, initvals, init, jitter_max_retries, n_init, trace, discard_tuned_samples, compute_convergence_checks, keep_warning_stat, return_inferencedata, idata_kwargs, nuts_sampler_kwargs, callback, mp_ctx, model, **kwargs) 650 auto_nuts_init = False 652 initial_points = None --> 653 step = assign_step_methods(model, step, methods=pm.STEP_METHODS, step_kwargs=kwargs) 655 if nuts_sampler != "pymc": 656 if not isinstance(step, NUTS): File [~/miniconda3/envs/pymc_env/lib/python3.11/site-packages/pymc/sampling/mcmc.py:211](https://vscode-remote+wsl-002bubuntu.vscode-resource.vscode-cdn.net/home/jausterw/work/playjuly2023/~/miniconda3/envs/pymc_env/lib/python3.11/site-packages/pymc/sampling/mcmc.py:211), in assign_step_methods(model, step, methods, step_kwargs) 209 methods_list: List[Type[BlockedStep]] = list(methods or pm.STEP_METHODS) 210 selected_steps: Dict[Type[BlockedStep], List] = {} --> 211 model_logp = model.logp() 213 for var in model.value_vars: 214 if var not in assigned_vars: 215 # determine if a gradient can be computed File [~/miniconda3/envs/pymc_env/lib/python3.11/site-packages/pymc/model.py:764](https://vscode-remote+wsl-002bubuntu.vscode-resource.vscode-cdn.net/home/jausterw/work/playjuly2023/~/miniconda3/envs/pymc_env/lib/python3.11/site-packages/pymc/model.py:764), in Model.logp(self, vars, jacobian, sum) 762 rv_logps: List[TensorVariable] = [] 763 if rvs: --> 764 rv_logps = transformed_conditional_logp( 765 rvs=rvs, 766 rvs_to_values=self.rvs_to_values, 767 rvs_to_transforms=self.rvs_to_transforms, 768 jacobian=jacobian, 769 ) 770 assert isinstance(rv_logps, list) 772 # Replace random variables by their value variables in potential terms File [~/miniconda3/envs/pymc_env/lib/python3.11/site-packages/pymc/logprob/basic.py:385](https://vscode-remote+wsl-002bubuntu.vscode-resource.vscode-cdn.net/home/jausterw/work/playjuly2023/~/miniconda3/envs/pymc_env/lib/python3.11/site-packages/pymc/logprob/basic.py:385), in transformed_conditional_logp(rvs, rvs_to_values, rvs_to_transforms, jacobian, **kwargs) 383 rvs_in_logp_expressions = _find_unallowed_rvs_in_graph(logp_terms_list) 384 if rvs_in_logp_expressions: --> 385 raise ValueError(RVS_IN_JOINT_LOGP_GRAPH_MSG % rvs_in_logp_expressions) 387 return logp_terms_list ValueError: Random variables detected in the logp graph: {S, normal_rv{0, (0, 0), floatX, False}.out, normal_rv{0, (0, 0), floatX, False}.out, halfnormal_rv{0, (0, 0), floatX, False}.out, halfnormal_rv{0, (0, 0), floatX, False}.out, normal_rv{0, (0, 0), floatX, False}.out, halfnormal_rv{0, (0, 0), floatX, False}.out, normal_rv{0, (0, 0), floatX, False}.out, halfnormal_rv{0, (0, 0), floatX, False}.out, halfnormal_rv{0, (0, 0), floatX, False}.out, normal_rv{0, (0, 0), floatX, False}.out, normal_rv{0, (0, 0), floatX, False}.out, halfnormal_rv{0, (0, 0), floatX, False}.out, halfnormal_rv{0, (0, 0), floatX, False}.out, normal_rv{0, (0, 0), floatX, False}.out, halfnormal_rv{0, (0, 0), floatX, False}.out, normal_rv{0, (0, 0), floatX, False}.out, normal_rv{0, (0, 0), floatX, False}.out, halfnormal_rv{0, (0, 0), floatX, False}.out, normal_rv{0, (0, 0), floatX, False}.out, halfnormal_rv{0, (0, 0), floatX, False}.out, halfnormal_rv{0, (0, 0), floatX, False}.out, normal_rv{0, (0, 0), floatX, False}.out, halfnormal_rv{0, (0, 0), floatX, False}.out, halfnormal_rv{0, (0, 0), floatX, False}.out, S, halfnormal_rv{0, (0, 0), floatX, False}.out, halfnormal_rv{0, (0, 0), floatX, False}.out, normal_rv{0, (0, 0), floatX, False}.out, halfnormal_rv{0, (0, 0), floatX, False}.out, halfnormal_rv{0, (0, 0), floatX, False}.out, halfnormal_rv{0, (0, 0), floatX, False}.out, normal_rv{0, (0, 0), floatX, False}.out, normal_rv{0, (0, 0), floatX, False}.out, halfnormal_rv{0, (0, 0), floatX, False}.out, normal_rv{0, (0, 0), floatX, False}.out, normal_rv{0, (0, 0), floatX, False}.out, normal_rv{0, (0, 0), floatX, False}.out, S, normal_rv{0, (0, 0), floatX, False}.out, halfnormal_rv{0, (0, 0), floatX, False}.out, halfnormal_rv{0, (0, 0), floatX, False}.out, normal_rv{0, (0, 0), floatX, False}.out, halfnormal_rv{0, (0, 0), floatX, False}.out, halfnormal_rv{0, (0, 0), floatX, False}.out, normal_rv{0, (0, 0), floatX, False}.out, S, S, halfnormal_rv{0, (0, 0), floatX, False}.out, normal_rv{0, (0, 0), floatX, False}.out, halfnormal_rv{0, (0, 0), floatX, False}.out, S, S, halfnormal_rv{0, (0, 0), floatX, False}.out, halfnormal_rv{0, (0, 0), floatX, False}.out, halfnormal_rv{0, (0, 0), floatX, False}.out, normal_rv{0, (0, 0), floatX, False}.out, halfnormal_rv{0, (0, 0), floatX, False}.out, halfnormal_rv{0, (0, 0), floatX, False}.out, halfnormal_rv{0, (0, 0), floatX, False}.out, normal_rv{0, (0, 0), floatX, False}.out, normal_rv{0, (0, 0), floatX, False}.out, halfnormal_rv{0, (0, 0), floatX, False}.out, S, normal_rv{0, (0, 0), floatX, False}.out, normal_rv{0, (0, 0), floatX, False}.out, normal_rv{0, (0, 0), floatX, False}.out, normal_rv{0, (0, 0), floatX, False}.out, halfnormal_rv{0, (0, 0), floatX, False}.out, normal_rv{0, (0, 0), floatX, False}.out, halfnormal_rv{0, (0, 0), floatX, False}.out, normal_rv{0, (0, 0), floatX, False}.out, halfnormal_rv{0, (0, 0), floatX, False}.out, normal_rv{0, (0, 0), floatX, False}.out, halfnormal_rv{0, (0, 0), floatX, False}.out, halfnormal_rv{0, (0, 0), floatX, False}.out, normal_rv{0, (0, 0), floatX, False}.out, halfnormal_rv{0, (0, 0), floatX, False}.out, halfnormal_rv{0, (0, 0), floatX, False}.out, halfnormal_rv{0, (0, 0), floatX, False}.out, normal_rv{0, (0, 0), floatX, False}.out, halfnormal_rv{0, (0, 0), floatX, False}.out, normal_rv{0, (0, 0), floatX, False}.out, normal_rv{0, (0, 0), floatX, False}.out, normal_rv{0, (0, 0), floatX, False}.out, halfnormal_rv{0, (0, 0), floatX, False}.out, normal_rv{0, (0, 0), floatX, False}.out, halfnormal_rv{0, (0, 0), floatX, False}.out, halfnormal_rv{0, (0, 0), floatX, False}.out, halfnormal_rv{0, (0, 0), floatX, False}.out, normal_rv{0, (0, 0), floatX, False}.out, normal_rv{0, (0, 0), floatX, False}.out, halfnormal_rv{0, (0, 0), floatX, False}.out, halfnormal_rv{0, (0, 0), floatX, False}.out, halfnormal_rv{0, (0, 0), floatX, False}.out, halfnormal_rv{0, (0, 0), floatX, False}.out, normal_rv{0, (0, 0), floatX, False}.out, normal_rv{0, (0, 0), floatX, False}.out, S, halfnormal_rv{0, (0, 0), floatX, False}.out, halfnormal_rv{0, (0, 0), floatX, False}.out, halfnormal_rv{0, (0, 0), floatX, False}.out, normal_rv{0, (0, 0), floatX, False}.out, normal_rv{0, (0, 0), floatX, False}.out, normal_rv{0, (0, 0), floatX, False}.out, halfnormal_rv{0, (0, 0), floatX, False}.out, normal_rv{0, (0, 0), floatX, False}.out, normal_rv{0, (0, 0), floatX, False}.out, halfnormal_rv{0, (0, 0), floatX, False}.out, halfnormal_rv{0, (0, 0), floatX, False}.out, normal_rv{0, (0, 0), floatX, False}.out, normal_rv{0, (0, 0), floatX, False}.out, S, halfnormal_rv{0, (0, 0), floatX, False}.out, normal_rv{0, (0, 0), floatX, False}.out, S, normal_rv{0, (0, 0), floatX, False}.out, normal_rv{0, (0, 0), floatX, False}.out, S, normal_rv{0, (0, 0), floatX, False}.out, normal_rv{0, (0, 0), floatX, False}.out, S, halfnormal_rv{0, (0, 0), floatX, False}.out, S, halfnormal_rv{0, (0, 0), floatX, False}.out, halfnormal_rv{0, (0, 0), floatX, False}.out, normal_rv{0, (0, 0), floatX, False}.out, halfnormal_rv{0, (0, 0), floatX, False}.out, halfnormal_rv{0, (0, 0), floatX, False}.out, halfnormal_rv{0, (0, 0), floatX, False}.out, normal_rv{0, (0, 0), floatX, False}.out, normal_rv{0, (0, 0), floatX, False}.out, normal_rv{0, (0, 0), floatX, False}.out, S, halfnormal_rv{0, (0, 0), floatX, False}.out, halfnormal_rv{0, (0, 0), floatX, False}.out, halfnormal_rv{0, (0, 0), floatX, False}.out, normal_rv{0, (0, 0), floatX, False}.out, normal_rv{0, (0, 0), floatX, False}.out, halfnormal_rv{0, (0, 0), floatX, False}.out, normal_rv{0, (0, 0), floatX, False}.out, halfnormal_rv{0, (0, 0), floatX, False}.out, S, halfnormal_rv{0, (0, 0), floatX, False}.out, halfnormal_rv{0, (0, 0), floatX, False}.out, S, S, normal_rv{0, (0, 0), floatX, False}.out, normal_rv{0, (0, 0), floatX, False}.out, halfnormal_rv{0, (0, 0), floatX, False}.out, normal_rv{0, (0, 0), floatX, False}.out, normal_rv{0, (0, 0), floatX, False}.out, halfnormal_rv{0, (0, 0), floatX, False}.out, halfnormal_rv{0, (0, 0), floatX, False}.out, normal_rv{0, (0, 0), floatX, False}.out, halfnormal_rv{0, (0, 0), floatX, False}.out, S, halfnormal_rv{0, (0, 0), floatX, False}.out, halfnormal_rv{0, (0, 0), floatX, False}.out, S, halfnormal_rv{0, (0, 0), floatX, False}.out, halfnormal_rv{0, (0, 0), floatX, False}.out, normal_rv{0, (0, 0), floatX, False}.out, normal_rv{0, (0, 0), floatX, False}.out, normal_rv{0, (0, 0), floatX, False}.out, normal_rv{0, (0, 0), floatX, False}.out, normal_rv{0, (0, 0), floatX, False}.out, S, normal_rv{0, (0, 0), floatX, False}.out, halfnormal_rv{0, (0, 0), floatX, False}.out, halfnormal_rv{0, (0, 0), floatX, False}.out, normal_rv{0, (0, 0), floatX, False}.out, halfnormal_rv{0, (0, 0), floatX, False}.out, halfnormal_rv{0, (0, 0), floatX, False}.out, halfnormal_rv{0, (0, 0), floatX, False}.out, normal_rv{0, (0, 0), floatX, False}.out, halfnormal_rv{0, (0, 0), floatX, False}.out, halfnormal_rv{0, (0, 0), floatX, False}.out, normal_rv{0, (0, 0), floatX, False}.out, normal_rv{0, (0, 0), floatX, False}.out, halfnormal_rv{0, (0, 0), floatX, False}.out, normal_rv{0, (0, 0), floatX, False}.out, normal_rv{0, (0, 0), floatX, False}.out, normal_rv{0, (0, 0), floatX, False}.out, normal_rv{0, (0, 0), floatX, False}.out, normal_rv{0, (0, 0), floatX, False}.out, halfnormal_rv{0, (0, 0), floatX, False}.out, normal_rv{0, (0, 0), floatX, False}.out, S, halfnormal_rv{0, (0, 0), floatX, False}.out, halfnormal_rv{0, (0, 0), floatX, False}.out, halfnormal_rv{0, (0, 0), floatX, False}.out, halfnormal_rv{0, (0, 0), floatX, False}.out, halfnormal_rv{0, (0, 0), floatX, False}.out, normal_rv{0, (0, 0), floatX, False}.out, normal_rv{0, (0, 0), floatX, False}.out, halfnormal_rv{0, (0, 0), floatX, False}.out, normal_rv{0, (0, 0), floatX, False}.out, normal_rv{0, (0, 0), floatX, False}.out, S, normal_rv{0, (0, 0), floatX, False}.out, S, halfnormal_rv{0, (0, 0), floatX, False}.out, normal_rv{0, (0, 0), floatX, False}.out, normal_rv{0, (0, 0), floatX, False}.out, normal_rv{0, (0, 0), floatX, False}.out, normal_rv{0, (0, 0), floatX, False}.out, S, normal_rv{0, (0, 0), floatX, False}.out, normal_rv{0, (0, 0), floatX, False}.out, halfnormal_rv{0, (0, 0), floatX, False}.out, halfnormal_rv{0, (0, 0), floatX, False}.out, S, normal_rv{0, (0, 0), floatX, False}.out, S, normal_rv{0, (0, 0), floatX, False}.out, halfnormal_rv{0, (0, 0), floatX, False}.out, halfnormal_rv{0, (0, 0), floatX, False}.out, halfnormal_rv{0, (0, 0), floatX, False}.out, halfnormal_rv{0, (0, 0), floatX, False}.out, normal_rv{0, (0, 0), floatX, False}.out, halfnormal_rv{0, (0, 0), floatX, False}.out, halfnormal_rv{0, (0, 0), floatX, False}.out, normal_rv{0, (0, 0), floatX, False}.out, halfnormal_rv{0, (0, 0), floatX, False}.out, halfnormal_rv{0, (0, 0), floatX, False}.out, normal_rv{0, (0, 0), floatX, False}.out, S, halfnormal_rv{0, (0, 0), floatX, False}.out, halfnormal_rv{0, (0, 0), floatX, False}.out, normal_rv{0, (0, 0), floatX, False}.out, halfnormal_rv{0, (0, 0), floatX, False}.out, normal_rv{0, (0, 0), floatX, False}.out, S, halfnormal_rv{0, (0, 0), floatX, False}.out, normal_rv{0, (0, 0), floatX, False}.out, normal_rv{0, (0, 0), floatX, False}.out, halfnormal_rv{0, (0, 0), floatX, False}.out, normal_rv{0, (0, 0), floatX, False}.out, halfnormal_rv{0, (0, 0), floatX, False}.out, normal_rv{0, (0, 0), floatX, False}.out, normal_rv{0, (0, 0), floatX, False}.out, normal_rv{0, (0, 0), floatX, False}.out, normal_rv{0, (0, 0), floatX, False}.out, halfnormal_rv{0, (0, 0), floatX, False}.out, halfnormal_rv{0, (0, 0), floatX, False}.out, halfnormal_rv{0, (0, 0), floatX, False}.out, halfnormal_rv{0, (0, 0), floatX, False}.out, normal_rv{0, (0, 0), floatX, False}.out, halfnormal_rv{0, (0, 0), floatX, False}.out, halfnormal_rv{0, (0, 0), floatX, False}.out, halfnormal_rv{0, (0, 0), floatX, False}.out, normal_rv{0, (0, 0), floatX, False}.out, halfnormal_rv{0, (0, 0), floatX, False}.out, S, S, S, halfnormal_rv{0, (0, 0), floatX, False}.out, halfnormal_rv{0, (0, 0), floatX, False}.out, normal_rv{0, (0, 0), floatX, False}.out, normal_rv{0, (0, 0), floatX, False}.out, normal_rv{0, (0, 0), floatX, False}.out, halfnormal_rv{0, (0, 0), floatX, False}.out, halfnormal_rv{0, (0, 0), floatX, False}.out, halfnormal_rv{0, (0, 0), floatX, False}.out, S, halfnormal_rv{0, (0, 0), floatX, False}.out, halfnormal_rv{0, (0, 0), floatX, False}.out, halfnormal_rv{0, (0, 0), floatX, False}.out, halfnormal_rv{0, (0, 0), floatX, False}.out, S, halfnormal_rv{0, (0, 0), floatX, False}.out, halfnormal_rv{0, (0, 0), floatX, False}.out, halfnormal_rv{0, (0, 0), floatX, False}.out, S, halfnormal_rv{0, (0, 0), floatX, False}.out, normal_rv{0, (0, 0), floatX, False}.out, normal_rv{0, (0, 0), floatX, False}.out, halfnormal_rv{0, (0, 0), floatX, False}.out, S, normal_rv{0, (0, 0), floatX, False}.out, S, normal_rv{0, (0, 0), floatX, False}.out, normal_rv{0, (0, 0), floatX, False}.out, halfnormal_rv{0, (0, 0), floatX, False}.out, normal_rv{0, (0, 0), floatX, False}.out, normal_rv{0, (0, 0), floatX, False}.out, normal_rv{0, (0, 0), floatX, False}.out,... This can happen when DensityDist logp or Interval transform functions reference nonlocal variables, or when not all rvs have a corresponding value variable.

(Note I cut a bit of the redundant error for space)

Thanks for your time
Joe

2

In case it is helpful, here’s the model definition:

EBM Math

Note this roughly follows [Redirecting](Fonteign 2012)

Definitions

The model assumes that a disease progresses according to a set of events E_1, \ldots, E_N, where N is the number of biomarkers. Our goal is to estimate an ordering S \in \sigma(N) over the events, which is a permutation of size N. The value of biomarker marker n for patient j is X_{nj} \in \mathbb{R}, which is a real-valued number. Each person j has a corresponding Bernoulli random variable d_j \in \{0,1\}, which denotes whether they have the disease or not (Note this is generally observed for these models, but I am including it for completeness). For person j (assumed to have the disease to simplify notation), k_n \in \{0, 1, \ldots, N\} denotes their current disease stage. Let \theta_n denote the parameters for the distribution of biomarker n when it is diseased and \phi_n be the corresponding parameters for when it is healthy.

Distributions

S is uniform over permutations.

p(X_{j} | S , z_j = 1, k_j) = \prod_{i=1}^{k_j}{p(X_{S(i)j} \mid \theta_n )} \prod_{i=k_j+1}^N{p(X_{S(i)j} \mid \phi_n)}

Note each term in a product is the likelihood of the biomarker value for participant j from either the diseased or healthy biomarker distributions.

As we don’t know what stage participant j is in the disease (we don’t know the value of k_j), we integrate over the possible values (0 to N).

P(X_{j} | z_j=1, S) = \sum_{k_j=0}^N{P(k_j) p(X_{j} \mid S, k_j)}

For simplicity, a uniform distribution over disease points is used.

Generative Process

We can write the EBM as the following generative process

S \sim {\rm UniformPermutation}(\cdot)

k_n \sim {\rm DiscreteUniform}(K)

X_{nj} | S, k_n \sim I(z_j == 1)I(j \lt k_n ) p(X_{nj} \mid \theta_n ) + \left(1-I(z_j==1) \right)I(j \geq k_n) p(X_{nj} \mid \phi_n)

3

This is adding a random variable in your generative graph. The logp will include it, since it’s not a model variable to be inferred. That’s where the error is coming from.

4

Thanks for the speedy response! I made a quick proper RV for the permutation rv, and took a subset of the data so that I could visualize the model graph. It’s clear something is wrong with the definition as there should be one S for all biomarkers, whereas the model graph has one S per biomarker (S should be outside of the plate in the figure rather than in the figure). I’m not sure why that’s happening as S is defined outside of the loop over participants. Do you have any suggestions?

Here’s the inferred model graph by pymc:

Screenshot 2023-07-11 111155
Screenshot 2023-07-11 111155965×321 17 KB

Here’s the code for defining the Permutation RV:

class PermutationRV(pm.distributions.distribution.Discrete):
    rv_op = pt.tensor.random.basic.permutation
    @classmethod
    def dist(cls, k, *args, **kwargs):
        k = pt.tensor.as_tensor_variable(k)
        return super().dist([k], *args, **kwargs)
    
    def moment(rv, k):
        # TODO: just return one permutation because they are all equally likely
        return pt.shared(np.range(k))

    def logp(value, k):
        #1 over number of permutations of size k -> 1/k factorial
        return -pt.tensor.gammaln(k+1)

and the updated code defining the model:


my_mod = pm.Model()
with my_mod:
    musd = pm.Normal.dist( mu=np.array([0., 0., 0.]), sigma=np.array([1, 1, 1]))
    sigsd = pm.HalfNormal.dist( sigma=np.array([2., 2., 2.]))
    musc = pm.Normal.dist(mu=np.array([-1., -1., -1]), sigma=np.array([1,1,1]))
    sigsc = pm.HalfNormal.dist(sigma=np.array([2.,2.,2.]))
    S = PermutationRV(name='S',k=K)

    for n in np.arange(N):
        n_k = pm.DiscreteUniform(f'n{n}_k', lower=0, upper=K-1)
        cur_mus = pt.tensor.switch(pt.tensor.ge(S,n_k), musd, musc)
        cur_sigs = pt.tensor.switch(pt.tensor.ge(S,n_k), sigsd, sigsc)
        if synth_data_np[n,K] == 0:
            cur_mus = musc
            cur_sigs = sigsc
        x = pm.Normal(f'x_{n}', mu=cur_mus, sigma=cur_sigs, observed=synth_data_np[n,0:K])

PS: I’m aware there’s probably a better way to parallelize this over people, but I wanted to try to get this to work first before adding in that extra step

PPS Sorry for the confusing notation – the code uses n/N to index participants and k/K for biomarkers/features, which is the reverse of the text/math description above.

5

Permutation returns a vector of size K if K is a scalar, so it seems plausible it would have shape=(3,)?

6

More importantly I’m afraid PyMC will not sample this variable well. For instance your logp does not do anything to stop repeated values or values out of bounds (just try to sample that variable alone with pm.sample and you’ll see it).

You may want to implement your custom step sampler for PermutationRVs which should be pretty simple, just return one random permutation.

7

They all do have shape 3 for this example. It’s a bit confusing because the last column of the data matrix is the classification, which isn’t considered a feature/biomarker here.

The RV samples using draw actually. The class method sets rv_op to pytensor.tensor.random.basic.permutation. If it has trouble once I get the model running, I’ll try and write a custom sampler for it.

Here is my updated logp for checking whether it is a valid permutation:

 def logp(val, k):
        #permutation sum needs to be 0.5*k*(k-1) 
        res = pt.tensor.switch(pt.tensor.and_(pt.tensor.eq(pt.tensor.max(val),k),
                                              pt.tensor.and_(pt.tensor.eq(pt.tensor.min(val),0),
                                                             pt.tensor.eq(pt.tensor.sum(val), 0.5 * k * (k-1)))),
                                              -pt.tensor.gammaln(k+1),
                                              -np.inf)
        #1 over number of permutations of size k -> 1/k factorial
        return res

The only thing it needs to check that it hasn’t is that all the values are integer (if the values are integers, the min is zero and the max is K-1, then it must be a permutation vector). I’m not sure how to check if an element is an integer though. Also, it’s possible I made a mistake in the code above, but it compiles and it isn’t where I’m getting an error yet.

I can draw values of S fine from the model (unconditioned), but I I’m still stuck on the same ValueError: Random variables detected in the logp graph. Something is off in the model definition as S should be outside of the plate (at least according to traditional graphical model figure notation).

I really appreciate your time and thanks for helping. I will post the notebook once it’s complete (and variants on it – I’m planning to write a Dirichlet process prior on permutations to find clusters of sequences).

8

To be a valid permutation you also must not have repeated values.

Metropolis will be incredibly inneficient because everything has to go right for it to propose different valid permutations.

I am conviced you need a custom step sampler. This would actually be a nice addition to PyMC.

The integer check shouldn’t be needed. I think PyMC knows this is a discrete variable in this case and will only propose discrete values. Worth printing value.dtype in the logp function to check.

On the shape thing, I am not sure what is your confusion, so at the risk of not helping at all I suggest you have a look at Distribution Dimensionality — PyMC 5.6.0 documentation