MvGaussianRandomWalk with multilevel

Questions
1

Hi,

I would like to fit only the bias of a multilevel model (level 1 with 20 categories and level 2 with 350) over time (15 periods)

To do that I did the following.

import theano.tensor as tt
import pymc3 as pm

model_randomwalk = pm.Model()
with model_randomwalk:
    # std of random walk
    sigma_bias_level_1 = pm.Exponential('sigma_alpha', 50., shape=20)
    cov_alpha = tt.diag(sigma_bias_level_1)

    sigma_bias_level_2 = pm.Exponential('sigma_beta', 50., shape=350)
    cov_beta = tt.diag(sigma_bias_level_2)
    
    bias_level_1 = pm.MvGaussianRandomWalk('alpha', cov=cov_alpha, shape=(15, 20)
    bias_level_2 = pm.MvGaussianRandomWalk('beta', cov=cov_beta, shape=(15, 350)

with model_randomwalk:
    regression = bias_level_1[index_ts, index_level_1] + bias_level_2[index_ts, index_level_2]

    sd = pm.HalfNormal('sd', sd=.1)
    likelihood = pm.Normal('y',
                           mu=regression,
                           sd=sd,
                           observed=y)

with model_randomwalk:
    trace_rw = pm.sample(tune=500, cores=2, init='auto')

But the sample method gets stuck at the beginning using jitter+adapt_diag. Is there something I could do better ?

Thank you

2

Is sampling works with init=‘adapt_diag’? The default jittering sometimes gives bad initial positions

3

Thank you for the quick answer.

Yes, it is working ! but with 40K observations it takes up to 4 hours (and I would like to add covariates in a second phase).

Is there a way to help the sampling by redefining the model ?

Thx

4

Looking solely on the model code you provide, I would suggest much more informative prior than pm.Exponential('sigma_alpha', 50., shape=20) - your current prior is likely too diffuse

5

Ok thanks, I am going to try that.

6

Hi,

This is an interesting model. I’m just learning about Pymc3 and trying to understand the syntax. Can you explain the line:

regression = bias_level_1[index_ts, index_level_1] + bias_level_2[index_ts, index_level_2]

I don’t understand how the indexing variables are defined. Are the defined somewhere not shown here? If so, how do you cycle through them?

Alternately, are these the names somehow assumed by Pymc3 internals? If so, how do you know how to name/define them?

Thanks,

Dorian