ricardoV94
commented
11 months ago The problem is that you're sampling a discrete variable with a generic mcmc sampler, and the sampler will sometimes propose values for p that sum to more or less than n. They would get rejected because they have 0 probability, but the logp graph errors out in the meantime. This is a known problem with using discrete variables structurally such as in indexing operations.
One proposed solution was to have discrete transforms: https://github.com/pymc-devs/pymc/pull/6102#issuecomment-1531558429 I guess we could also have the Metropolis sampler catch errors and treat them as 0 probability, but that could mask bugs in the user code.
For your case you can post-process p so bad values never make it to the indexes. It shouldn't bias your results because those draws would never be accepted in the first place. Something like:
import arviz, pymc, numpy, pytensor
data=numpy.concatenate([numpy.random.randn(300)/10-0.1,numpy.random.randn(450)/10+0.2,numpy.random.randn(250)/10+0.6])
with pymc.Model() as model:
K, N = 3, 1000
r = pymc.Normal("r", mu=0, sigma=1, shape=K)
w = pymc.Dirichlet("w", numpy.ones(K))
p = pymc.DirichletMultinomial("p", N, w)
# Avoid index errors: if p does not sum to N, replace by [N, 0, 0, ...]
p = pymc.math.switch(pymc.math.eq(p.sum(), N), p, [N] + [0]*(K-1))
idx = pymc.Deterministic("idx", pymc.math.concatenate(
[pytensor.tensor.tile(i,p[i]) for i in range(K)]))
obs = pymc.Normal("obs", mu=r[idx], sigma=1e-3, observed=data)
trace = pymc.sample(1000) The model will probably sample horribly, because the DirichletMultinomial is very constrained and the Metropolis sampler doesn't know about these contraints. It may be better to transform the Dirichlet deterministically into a valid p
Describe the issue:
I wanted to test the detection for several breaking points in a time series, so I came up with the following model (cf. code example) to detect 3 switching points. As long as I don't add any observation data, the model runs fine and generates all the variables at sampling time. The idx tensor variable seems to be of the proper shape (1 dimension and size of 1000). However, adding them generates an error (cf. message). The dimension mismatch is also variable, depending on the execution, it can say input[1].shape[0] = 998, or input[1].shape[0] = 1004, it seems pretty random.
The problem seems to be with the idx index, but generating it in another way works. If I do for instance
everything ends up being fine.
This index is basically used to map the proper value of the mean on a given interval. Based on the outputs from the DirichletMultinomial, which is a tuple of K numbers ($a_0$, $a1$,..., $a{K-1}$) such as $$\sum_{i=0}^{K-1}a_i=N$$ I create an array that would look like [0,0,0,1,1,1,1,1,........,K-1,K-1,K-1]. 0 is repeated $a_0$ times, 1 is repeated $a_1$ times, etc...
Reproduceable code example:
Error message:
PyMC version information:
I am running a conda environment, set up from miniforge.
Context for the issue:
No response