Mixture of normal and AR(p) process?

Questions
1

I’m trying to estimate a model that is a pm.Mixture of a normal distribution and an AR process. How would I do that in PyMC3? My current attempt is:

        with pm.Model() as self.model:
            W = np.array([1., 1.])
            w = pm.Dirichlet('w', W)

            intercept = pm.Normal('intercept', mu=-5, sd=5., testval=-5.)
            rhos = pm.Uniform('rhos', lower=-1., upper=1., shape=self.num_lags)

            comp1 = pm.AR.dist(rho=rhos, sd=pm.math.exp(intercept))
            comp2 = pm.Normal.dist(mu=0., sd=0.01)

            pm.Mixture('obs', w=w, comp_dists=[comp1, comp2], observed=data)

But I run into the following error:

File “…”, line 111, in fit_model
pm.Mixture(‘obs’, w=w, comp_dists=[comp1, comp2], observed=data)
File “/usr/local/lib/python2.7/dist-packages/pymc3/distributions/distribution.py”, line 42, in new
return model.Var(name, dist, data, total_size)
File “/usr/local/lib/python2.7/dist-packages/pymc3/model.py”, line 838, in Var
total_size=total_size, model=self)
File “/usr/local/lib/python2.7/dist-packages/pymc3/model.py”, line 1319, in init
self.logp_elemwiset = distribution.logp(data)
File “/usr/local/lib/python2.7/dist-packages/pymc3/distributions/mixture.py”, line 145, in logp
return bound(logsumexp(tt.log(w) + self._comp_logp(value), axis=-1),
File “/usr/local/lib/python2.7/dist-packages/pymc3/distributions/mixture.py”, line 115, in _comp_logp
axis=1))
File “/usr/local/lib/python2.7/dist-packages/theano/tensor/basic.py”, line 4709, in stack
return join(axis, *[shape_padaxis(t, axis) for t in tensors])
File “/usr/local/lib/python2.7/dist-packages/theano/tensor/basic.py”, line 4601, in shape_padaxis
raise IndexError(msg)
IndexError: axis 1 is out of bounds [-1, 1)

2

AR§ is a vector, so you need to specific the shape in both comp_dist. It might be easier to write it down using a potential or DensityDist.

3

I get the same error when I specify the shape: here is what I tried (maybe I’m not doing it correctly):

    with pm.Model() as self.model:

        W = np.array([1., 1.])
        w = pm.Dirichlet('w', W)

        intercept = pm.Normal('intercept', mu=-5, sd=5., testval=-5.)
        rhos = pm.Uniform('rhos', lower=-1., upper=1., shape=self.num_lags)

        comp1 = pm.AR.dist(rho=rhos, sd=pm.math.exp(intercept), shape=len(data))
        comp2 = pm.Normal.dist(mu=0., sd=0.01, shape=len(data))

        pm.Mixture('obs', w=w, comp_dists=[comp1, comp2], observed=data)

I still get the same axis 1 is out of bounds [-1, 1) error. A custom likelihood with DensityDist is the next thing I’m going to try, I think.