Error when observing on uniform with sampled parameters in PyMC

I’m new to PyMC and am trying to model a situation where you are rolling marbles at a wall and trying to find the block. The data is only for the values where the marble hits the block.

I’m first sampling an x location and size and then calculating a point from those with a uniform, but I’m getting an error.

import pymc3 as pm
import theano.tensor as tt

basic_model = pm.Model()

with basic_model:

    # We are assuming independence of these.
    x = pm.Uniform("x", lower=1, upper=30)
    l = pm.Uniform("l", lower=1, upper=30)

    lower = pm.Deterministic('lower', x-0.5*l)
    upper = pm.Deterministic('upper', x+0.5*l)

    point_x = pm.Uniform('point_x', lower=lower, upper=upper, observed=x_vals)

    pm.sample()

With the error:

SamplingError: Initial evaluation of model at starting point failed!
 Starting values:
{'x_interval__': array(0.), 'l_interval__': array(0.)}

 Initial evaluation results:
 x_interval__   -1.39
 l_interval__   -1.39
 point_x         -inf
 Name: Log-probability of test_point, dtype: float64

Clearly the issue is with point_x. I’m guessing the error has to do with the fact that the observed data may potentially fall outside the lower-upper range depending on the value of x and l sampled. But how might I fix this?

Solution

The sampler doesn't know how to handle starting off in an invalid region of the parameter space. A quick and dirty fix is to provide testval arguments that ensure the sampling begins in a logically valid solution. For example, we know the minimum block must have:

l_0 = np.max(x_vals) - np.min(x_vals)
x_0 = np.min(x_vals) + 0.5*l_0

and could use those:

x = pm.Uniform("x", lower=1, upper=30, testval=x_0)
l = pm.Uniform("l", lower=1, upper=30, testval=l_0)

Also, the nature of this model leads to many rejections due to impossibility, so you may want to use Metropolis for sampling, which almost always needs more steps and tuning

pm.sample(tune=10000, draws=10000, step=pm.Metropolis())

Alternative Models

Otherwise, consider reparameterizing so that only valid solutions are in the parameter space. One approach would be to sample l and then use that to constrain x. Something like:

other_model = pm.Model()

x_min = np.min(x_vals)
x_max = np.max(x_vals)
l_0 = x_max - x_min

with other_model:

    # these have logical constraints from the data
    l = pm.Uniform("l", lower=l_0, upper=30)
    x = pm.Uniform("x", lower=x_max - 0.5*l, upper=x_min + 0.5*l)
    
    lower = pm.Deterministic('lower', x - 0.5*l)
    upper = pm.Deterministic('upper', x + 0.5*l)

    point_x = pm.Uniform('point_x', lower=lower, upper=upper, observed=x_vals)

    res = pm.sample(step=pm.NUTS(), return_inferencedata=True)

Another approach would be to sample lower and upper directly, and compute the x and l as deterministic variables from those.