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

Questions
8

Dear Adrian,

I got your model to work easily. However, the problem is not gone: it’s maybe just happening more rarely. See this trace (at ca. 4500):

bhglm_measurement1_ASsuggestion
bhglm_measurement1_ASsuggestion1920×948 568 KB

The code is below, I only loosened prior constraints.

If what Junpeng suggests is true, and it sounds plausible, then the reason that your reparametrization did not show any flat traces in the first place was maybe the narrow prior range. I changed this moderately (sd=10), because I think it is good to load as few assumptions as possible.

further questions

I did not understand the factor of 0.2 in your first model.

Could you please explain where it comes from?

Did I get it right, that your suggestion did not account for a potentially different standard deviation of the subjects? With the global intercept and the “subject_sd * subject”, only an individual intercept is incorporated, but no individual variation. This might be reasonable, yet I’m trying to understand if it might be another involuntary change that made sampling “easier”.

Finally, your “scale-mixture prior” could be a good alternative. In fact, our hypothesis is that the “condition” does nothing and that the outlier is a “weirdo”. However, I do not get the meaning of tau/teta/eta. Do you have any reference or example where I could read this up?

Best,

Falk

Here’s how I changed the priors (TT = theano.tensor, NP = numpy, PM = pymc3.1 master):

    # Model
    with PM.Model() as hierarchical_model:
        intercept = PM.Normal('intercept', mu=0., sd=10)

        subject_sd_raw = PM.Uniform('subject_sd_raw', lower = 0, upper=NP.pi/2)
        subject_sd = 0.2*TT.tan(subject_sd_raw)
        subject = PM.Normal('subject', mu=0, sd=10, shape=n_subjects)

        altered_mu = PM.Normal('altered_mu', mu=0, sd=10)
        altered_subject_sd_raw = PM.Uniform('altered_subject_sd_raw', lower = 0, upper=NP.pi/2)
        altered_subject_sd = TT.tan(altered_subject_sd_raw)

        altered_subject = PM.Normal('altered_subject', mu=0, sd=10, shape=n_subjects)

        is_altered = data['is_altered'].values
        estimator = (intercept
                     + subject_sd * subject[subject_idx]
                     + altered_mu * is_altered
                     + altered_subject_sd * altered_subject[subject_idx] * is_altered)

        residual = PM.HalfCauchy('residual', 5)
        nu = PM.Gamma('nu', 2, 0.1)
        likelihood = PM.StudentT('y', nu=nu, mu=estimator,
                                 sd=residual, observed=data[predicted_variable])
        trace = PM.sample( \
                          draws = 10000 \
                        , tune = 1000 \
                        , njobs = 4 \
                        , nuts_kwargs={'target_accept': 0.95} \
                        )