# Mixture of normal and AR(p) process?

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)

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.

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.