AR and AR1 models fail with default sampling

[phil-mcdowall]

AR and AR1 models fail with default sampling (and NUTS/Hamiltonian in general), raising a value error
ValueError: Mass matrix contains zeros on the diagonal. The derivative of RV beta.ravel()[0] is zero. Metropolis sampling seems to work for these models. Not sure if this is a bug or unsuitable sampling, but I've seen a few examples of AR models fit with default sampling, so this may just be the result of changes in default sampling.

Examples below...

# AR1 model with NUTS
beta_true = 0.8
y_true = np.zeros(100)
for i in range(1,len(y_true)):
    y_true[i] = beta_true * y_true[i-1] + np.random.normal(0,1)

with pymc3.Model() as AR1_model:
    beta = pymc3.Normal("beta",mu=0, sd=2)
    obs = pymc3.AR1("obs", k=beta, tau_e=1,observed=y_true)
    trace = pymc3.sample()

# AR model with NUTS
alpha_true = 0.2
beta1_true = 0.7
beta2_true = 0.05
y_true = np.zeros(100)
for i in range(2,len(y_true)):
    y_true[i] = alpha_true + beta1_true * y_true[i-1] + beta2_true * y_true[i-2] + np.random.normal(0,1)
y_true = np.array(y_true[1:])

with pymc3.Model() as AR_model:
    beta = pymc3.Uniform("beta",lower=0, upper=1, shape=2)
    obs = pymc3.AR("obs", rho=beta, sd=0.01, constant=True, observed=y_true)
    trace = pymc3.sample()

pymc3 version: '3.4.1'
theano version: '1.0.2'

EDIT:
Initializing with ADVI allows sampling to complete, but results in an acceptance rate of 1, and incorrect estimates.
The acceptance probability does not match the target. It is 1.0, but should be close to 0.8. Try to increase the number of tuning steps.

[junpenglao]

Anybody else reproduces the same problem? Runs fine locally for me.

[phil-mcdowall]

Tried a fresh install using pip rather than the master and I'm still getting the same error, minus the updated error message on the master that gives the RV causing the issue.
Full traceback:

Traceback (most recent call last):
  File "<stdin>", line 4, in <module>
  File "/home/phil/pymc3_fresh/lib/python3.6/site-packages/pymc3/sampling.py", line 462, in sample
    trace = _sample_many(**sample_args)
  File "/home/phil/pymc3_fresh/lib/python3.6/site-packages/pymc3/sampling.py", line 510, in _sample_many
    step=step, random_seed=random_seed[i], **kwargs)
  File "/home/phil/pymc3_fresh/lib/python3.6/site-packages/pymc3/sampling.py", line 554, in _sample
    for it, strace in enumerate(sampling):
  File "/home/phil/pymc3_fresh/lib/python3.6/site-packages/tqdm/_tqdm.py", line 930, in __iter__
    for obj in iterable:
  File "/home/phil/pymc3_fresh/lib/python3.6/site-packages/pymc3/sampling.py", line 652, in _iter_sample
    point, states = step.step(point)
  File "/home/phil/pymc3_fresh/lib/python3.6/site-packages/pymc3/step_methods/arraystep.py", line 222, in step
    apoint, stats = self.astep(array)
  File "/home/phil/pymc3_fresh/lib/python3.6/site-packages/pymc3/step_methods/hmc/base_hmc.py", line 115, in astep
    self.potential.raise_ok()
  File "/home/phil/pymc3_fresh/lib/python3.6/site-packages/pymc3/step_methods/hmc/quadpotential.py", line 190, in raise_ok
    raise ValueError('Mass matrix contains zeros on the diagonal. '
ValueError: Mass matrix contains zeros on the diagonal. Some derivatives might always be zero.

[fonnesbeck]

The examples above work just fine for me, but the first example (ar1) on the demo notebook in the repo fails for me, but with a different error:

---------------------------------------------------------------------------
IndexError                                Traceback (most recent call last)
<ipython-input-3-6ad2e5b82206> in <module>()
      2 with pm.Model() as ar1:
      3     beta = pm.Normal('beta', mu=0, sd=tau)
----> 4     data = pm.AR('y', beta, sd=1.0, observed=y)
      5     trace = pm.sample(1000, cores=4)
      6 

~/anaconda3/envs/dev/lib/python3.6/site-packages/pymc3/distributions/distribution.py in __new__(cls, name, *args, **kwargs)
     39                 raise TypeError("observed needs to be data but got: {}".format(type(data)))
     40             total_size = kwargs.pop('total_size', None)
---> 41             dist = cls.dist(*args, **kwargs)
     42             return model.Var(name, dist, data, total_size)
     43         else:

~/anaconda3/envs/dev/lib/python3.6/site-packages/pymc3/distributions/distribution.py in dist(cls, *args, **kwargs)
     50     def dist(cls, *args, **kwargs):
     51         dist = object.__new__(cls)
---> 52         dist.__init__(*args, **kwargs)
     53         return dist
     54 

~/anaconda3/envs/dev/lib/python3.6/site-packages/pymc3/distributions/timeseries.py in __init__(self, rho, sd, tau, constant, init, *args, **kwargs)
    105         else:
    106             try:
--> 107                 p = rho.shape.tag.test_value[0]
    108             except AttributeError:
    109                 p = rho.shape[0]

IndexError: index 0 is out of bounds for axis 0 with size 0

The ar2 models run just fine, however.

[junpenglao]

Oh that one is on me - will send a fix.

[phil-mcdowall]

Some further experimentation reveals that the issue is in no way related to AR1 models, and is in fact a complete failure of NUTS. The following model also fails with a mass matrix related ValueError...

y = np.random.normal(10,1,1000)

with pymc3.Model() as model:
    mean = pymc3.Normal('mean',0,100)
    obs = pymc3.Normal('obs',mu=mean,sd=1,observed=y)
    trace = pymc3.sample()

Will try a fresh install of pymc3 and all dependencies on a different system.

[phil-mcdowall]

Installs on windows, and Ubuntu 16 working fine. Couldn't get any install of pymc3 working on Ubuntu 18.04 but it seems like its probably something weird and specific to my system that resulted in all samples being identical, leading to zeros on the mass matrix diagonal.
With the sampling working the estimates are great and NUTS is really fast. Thanks!