Variable Lag and slicing

Questions version agnostic
7

The problem has to do with how discrete proposals are attempted by PyMC. The metropolis sampler will sometimes propose lags of 0 or -1, and check the logp delta to decide whether to accept or reject that proposal. Under “normal” models, the logp of the proposal will be -inf, and it will be promptly rejected.

But in this case it first leads p to have the wrong shape and the model logp to fail to evaluate at all. Tthe acceptance ratio takes the whole model logp in consideration, even if in this case the prior probability would be enough to rule out such proposal.

You can observe something similar by wrapping lag in a PyTensor assert:

import pymc as pm
import numpy as np
from pytensor.raise_op import assert_op

n_obs = 10000
true_lag = 2
lead = np.random.normal(loc=0, scale=0.5, size=n_obs + true_lag)
y_obs = -2.0 * lead[:-true_lag] + np.random.normal(loc=0, scale=0.4, size=n_obs)
lead = lead[true_lag:]

with pm.Model() as model:
    lead = pm.ConstantData("lead", value=lead)
    alpha = pm.Normal("alpha")
    lag = pm.DiscreteUniform("lag", 1, 3)
    lag = assert_op(lag, lag > 0)
    p = pm.Deterministic("p", lead[12-lag:-lag])
    sigma = pm.HalfNormal("sigma")
    y = pm.Normal("y", alpha * p, sigma=sigma)
    
    trace = pm.sample(chains=1)

AssertionError: PyTensor Assert failed!
Apply node that caused the error: Assert{msg=PyTensor Assert failed!}(Reshape{0}.0, Elemwise{gt,no_inplace}.0)
Toposort index: 21
Inputs types: [TensorType(int64, ()), TensorType(bool, ())]
Inputs shapes: [(), ()]
Inputs strides: [(), ()]
Inputs values: [array(0), array(False)]

Check the Assert failing with the proposed value of 0 array(0), which would otherwise lead to a weird shape in the indexing of lead.

np.arange(100)[12-0:-0]  # [ ]

A solution is to build your model to be robust to stupid proposals:

import pymc as pm
import numpy as np
import pytensor.tensor as pt

n_obs = 10000
true_lag = 2
lead = np.random.normal(loc=0, scale=0.5, size=n_obs + true_lag)
y_obs = -2.0 * lead[:-true_lag] + np.random.normal(loc=0, scale=0.4, size=n_obs)
lead = lead[true_lag:]

with pm.Model() as model:
    lead = pm.ConstantData("lead", value=lead)
    alpha = pm.Normal("alpha")
    lag = pm.DiscreteUniform("lag", 1, 3)
    p = pm.Deterministic("p", lead[12 - pt.max((1, lag)):-pt.max((1, lag))])
    sigma = pm.HalfNormal("sigma")
    y = pm.Normal("y", alpha * p, sigma=sigma)

    trace = pm.sample(chains=1)

We have considered other options, such as implementing specialized samplers for bounded discrete variables or even discrete variable transformations: Implement transforms for discrete variables by ricardoV94 · Pull Request #6102 · pymc-devs/pymc · GitHub

Actually if you had used a Categorical prior you wouldn’t have seen this issue, because it uses a less dumb sampler.

1 Like