Dear developers,
I’ve noticed a problem that is most likely sampler related and was introduced with the upgrade to pymc3.1. This happens on sampling a hierarchical (two level) GLM.
Description
I observe that the sampler occasionally gets stuck in parameter space (not: temporal slowdown), exploring only a very restricted fraction of possible values, before returning to the “hairy caterpillar” behavior. This temporary local “attractor” can be within or outside the variance range of the rest of the trace. If it happens, all model parameters are affected together.
In addition, I get an “acceptance probability” warning, even if I choose extensive tuning and a high “target_accept”.
Note that in Version 3.0, there are also occasional “spikes” on the trace. However, the sampler returns to the usual range much more quickly.
Please find here some additional documents:

minimal code example LocomotorAnalysis.py (8.7 KB)

data example [20170613_data.csv (2.7 KB)
output.txt] 
and also:
 traceplot from pymc3 version 3.0: http://amor.cms.huberlin.de/~mielkefa/bhglm_measurement1_6666steps_pymc30.png  traceplot from pymc3 version 3.1: http://amor.cms.huberlin.de/~mielkefa/bhglm_measurement1_6666steps_pymc31.png  console outputs on my machine http://amor.cms.huberlin.de/~mielkefa/output.txt
(sorry for text links and external storage, this is to circumvent the “two links” restriction for newbies)
Model configuration
My model follows introductory examples from Thomas Wiecki and others (“bayesianglms3”, “bayesianhierchicalnoncentered”).
When using an “offset” (noncentered) reparametrization suggested by Thomas Wiecki in the latter notebook, the problem seems to occur less often and only for shorter periods (unquantified). Hence I suspect this might be related to sampler coverage. Also, when I switch to Slice sampler (default is NUTS) this still happens, but less frequently. The problem persists through the very recent switch do a default NUTS initialization via “jitter+adapt_diag”.
Data
The data is limited to four subjects and twenty repeats per condition for each of them; I’ve reduced it to only one test condition and one measured value. Regarding model specification, I’ve tried many variants, but I am happy if anyone can point me at a conceptional bug. My knowledge leaves me when it comes to sampler properties and step choice. The reference model definition in the attached script can also be found below.
Expectation
If somebody with a deeper knowledge of sampling procedures could take an educated guess on what’s going wrong, that would already help.
I would hope to get help with the model definition to get this working with the little data we have.
Alternatively, maybe someone knows a configuration parameter in “sample” that can restore the old behavior?
System
I install the different pymc versions as follows:
 pymc3.1 latest via “git clone”/“git pull” and “pip install .”
 pymc3.0 via "pip install pymc3==3.0"
Python is python 3.6; all required packages are up to date.
My system is Manjaro Linux; Kernel 4.9.39.
Acknowledgements
Thomas Wiecki was very supportive when I sent him my problem previously, but suggested to bring this here. Any help is much appreciated.
Thanks to everyone in advance!
Falk
with PM.Model() as hierarchical_model:
hypermu_intercept = PM.Normal(‘hypermu intercept’, mu=0., sd=100)
hypersigma_intercept = PM.HalfNormal(‘hypersigma intercept’, sd=100)
hypermu_slope = PM.Normal(‘hypermu slope’, mu=0., sd=100)
hypersigma_slope = PM.HalfNormal(‘hypersigma slope’, sd=100)residual = PM.HalfCauchy('residual', 5) intercept = PM.Normal('intercept', mu=hypermu_intercept, sd=hypersigma_intercept, shape=n_subjects) slope = PM.Normal('slope', mu=hypermu_slope, sd=hypersigma_slope, shape=n_subjects) estimator = intercept[subject_idx] \ + slope[subject_idx] * data.is_altered.values likelihood = PM.Normal( 'likelihood' \ , mu = estimator \ , sd = residual \ , observed = data[predicted_variable]\ )