Sampler occasionally gets static in pymc3.1, but not in pymc3.0

Questions
13

Dear all,

I managed to solve this problem (code below for reference). Yet some questions remain.
It turns out that the divergent periods disappear if I remove the hyperprior for standard deviation (hypersigma_slope in my original model). The working model now still uses non-centered reparametrization and has moderately wide prior distributions. I still get spikes, yet they are not persistent.

I guess the data sample size is just at the edge of being insufficient for the model. Feel free to close the thread, but see below.

Sorry for the lack of structure in my previous post. There are two general questions about which I would appreciate a brief discussion:
(1) Is it possible in pymc3.1 master, via a sampler parameter or general setting, to fix step size (mimic v3.0 behavior), ignoring the funnel at own risk? (assuming we currently lack a more sophisticated way to avoid that the sampler gets trapped in some weird local geometry)
(2) Is there a general method to evaluate the relation of data coverage and model complexity? Can we test if data is sufficient to distinguish between two models? (see Boettiger paper from my previous comment)

Thanks for everyone’s input. I learned a lot!
Best,

Falk

Model

with PM.Model() as hierarchical_model:
    ### population intercept
    intercept_population = PM.Normal('intercept_population', mu = 0., sd = 10)
    sd_population = PM.HalfCauchy('sd_population', beta = 5)

    ### individual means
    intercept_subject = PM.Normal('intercept_subject', mu = intercept_population, sd = sd_population, shape = n_subjects)

    ### priors for "altered"
    # population level
    slope_altered_population = PM.Normal('slope_altered_population', mu = 0., sd = 10)

    # individual level
    offset_altered_subject = PM.Normal('offset_altered_subject', mu = slope_altered_population, sd = 10, shape = n_subjects)
    sd_altered_subject = PM.HalfCauchy('sd_altered_subject', beta = 5)

    slope_altered_subject = PM.Deterministic('slope_altered_subject', offset_altered_subject * sd_altered_subject)
    

    ## estimator
    is_altered = data['is_altered'].values
    estimator = intercept_subject[subject_idx] \
              + slope_altered_subject[subject_idx] * is_altered \

    residual = PM.HalfCauchy('residual', 5)
    dof = PM.Gamma('dof', 2, 0.1)
    likelihood = PM.StudentT( 'likelihood' \
                            , nu = dof \
                            , mu = estimator \
                            , sd = residual \
                            , observed = data[predicted_variable] \
                            )

    return hierarchical_model