I’m trying to implement a correlated random walk over the slope and intercept parameters of a linear model. I’m basing it loosely off the multivariate random walk tutorial on the pymc site. However, I keep running into some errors when trying to sample the prior predictive and posterior distributions. I’m relatively new to PyMC, so any help is greatly appreciated. Here is some sample code that throws the errors:
#simple linear model
x = np.random.uniform(-1, 1, 100)
x.sort(axis=0)
y = x * 10 + np.random.normal(0,.1,100)
The model I’m trying to fit looks like this:
with pm.Model() as mv_randomwalk:
# x = pm.MutableData('x', x_t)
# Sigma for mu:
sigma = pm.Normal('sigma', 0, 1)
# LKJCholeskyCov for the covariance between alpha and beta
chol_alpha_beta, *_ = pm.LKJCholeskyCov(
'chol_alpha_beta', n=2, eta=2, sd_dist=pm.HalfCauchy.dist(2.5), compute_corr=True
)
# Correlated random walks for alpha and beta
alpha_beta = pm.MvGaussianRandomWalk(
'alpha_beta', mu=np.zeros(2), chol=chol_alpha_beta, shape=(x.shape, 2),
)
#extract params
alpha = pm.Deterministic('alpha', alpha_beta[:, 0])
beta = pm.Deterministic('beta', alpha_beta[:, 1])
print(x.ndim)
#regression
mu = beta * x + alpha
#likelihood
y_hat = pm.Normal('y_hat', mu, sigma=sigma, observed=y_t)
with mv_randomwalk:
idata = pm.sample_prior_predictive()
with mv_randomwalk:
idata.extend(pm.sample(nuts_sampler='numpyro'))
Here is the error I get when sampling the prior predictive:
ValueError: scale < 0
Apply node that caused the error: normal_rv{0, (0, 0), floatX, True}(RandomGeneratorSharedVariable(<Generator(PCG64) at 0x2374384A0>), [100], 11, Composite{((i0 * i1) + i2)}.0, sigma)
Toposort index: 14
Inputs types: [RandomGeneratorType, TensorType(int64, shape=(1,)), TensorType(int64, shape=()), TensorType(float64, shape=(100,)), TensorType(float64, shape=())]
Inputs shapes: ['No shapes', (1,), (), (100,), ()]
Inputs strides: ['No strides', (8,), (), (8,), ()]
Inputs values: [Generator(PCG64) at 0x2374384A0, array([100]), array(11), 'not shown', array(-0.28429503)]
Outputs clients: [['output'], ['output']]
HINT: Re-running with most PyTensor optimizations disabled could provide a back-trace showing when this node was created. This can be done by setting the PyTensor flag 'optimizer=fast_compile'. If that does not work, PyTensor optimizations can be disabled with 'optimizer=None'.
HINT: Use the PyTensor flag `exception_verbosity=high` for a debug print-out and storage map footprint of this Apply node.
And when I try to run the NUTS sampler, I get this dimensions-related error:
ValueError: Input dimension mismatch: (input[0].shape[1] = 3, input[1].shape[1] = 75)
Apply node that caused the error: Composite{switch(i4, (((-0.5 * sqr(((i0 - i1) / i2))) - 0.9189385332046727) - i3), -inf)}(<Matrix(float64, shape=(?, ?))>, Add.0, ExpandDims{axes=[0, 1]}.0, Log.0, Gt.0)
Toposort index: 23
Inputs types: [TensorType(float64, shape=(None, None)), TensorType(float64, shape=(1, 75)), TensorType(float64, shape=(1, 1)), TensorType(float64, shape=(1, 1)), TensorType(bool, shape=(1, 1))]
Inputs shapes: [(75, 3), (1, 75), (1, 1), (1, 1), (1, 1)]
Inputs strides: [(24, 8), (600, 8), (8, 8), (8, 8), (1, 1)]
Inputs values: ['not shown', 'not shown', array([[0.77063881]]), array([[-0.26053548]]), array([[ True]])]
Outputs clients: [[Sum{axes=None}(sigma > 0)]]
HINT: Re-running with most PyTensor optimizations disabled could provide a back-trace showing when this node was created. This can be done by setting the PyTensor flag 'optimizer=fast_compile'. If that does not work, PyTensor optimizations can be disabled with 'optimizer=None'.
HINT: Use the PyTensor flag `exception_verbosity=high` for a debug print-out and storage map footprint of this Apply node.
I’m particularly confused about where this input[0].shape[1] = 3 is coming from. It seems like the result of MV_GaussianRandomwalk must have a shape that is different from what I would expect, (something like (100, 1, 3)). I would expect a shape of (100,) corresponding to one of the two (100, 2) random walk dimensions specified. When I naively try to inspect the shape of beta: print(beta.shape) I get Shape.0.
I have a feeling I have two related issues here- first, the model is probably mis-specified in some way. I suspect I am not indexing the alpha and beta randomwalks correctly. Secondly, as a result, the x array and the alpha and beta variables are not shaped the same. I’m relatively new to PyMC (and PyTensor), and could use some handholding here- how can I inspect the output of PyMC variables to see things like size, shape etc? What else am I doing wrong here? Thanks for the great package and for any help.