Interpolated PDF in a Mixture Model

Questions
1

Hi,

I’m trying to use an interpolated PDF in a mixture model with some Gaussian PDFs to fit a spectrum. I’ve attached a MWE.

import pymc3 as pm
import numpy as np

x = np.array([0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, 10.5, 11.5, 12.5, 13.5, 14.5, 15.5, 16.5, 17.5, 18.5])
y = np.array([1.72e-02, 1.99e-02, 2.06e-02, 2.04e-02, 1.95e-02, 1.81e-02, 1.66e-02, 1.48e-02, 1.30e-02, 1.11e-02, 9.22e-03, 7.41e-03, 5.72e-03, 4.18e-03, 2.84e-03, 1.72e-03, 8.58e-04, 2.78e-04, 1.54e-05])

with pm.Model() as model:
    mu = pm.Normal("mu0", mu=5., sd=1., shape=(1,))
    sd = pm.HalfCauchy("sd0", beta=1., shape=(1,))
    gaus = pm.Normal.dist(mu=mu, sd=sd, shape=(1,))
    interp_pdf = pm.Interpolated.dist(x_points=x, pdf_points=y, testval=2.5)
    w = pm.Dirichlet('w', a=np.array([1., 1.]))
    like = pm.Mixture('like', w=w, comp_dists=[gaus, interp_pdf], testval=2.5)

I’m getting the error:

Traceback (most recent call last):
  File "/Users/brianzhu/anaconda/envs/py36/lib/python3.6/site-packages/pymc3/distributions/mixture.py", line 110, in _comp_logp
    return comp_dists.logp(value_)
AttributeError: 'list' object has no attribute 'logp'

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "PyMCInterpolateTest.py", line 41, in <module>
    like = pm.Mixture('like', w=w, comp_dists=[gaus, interp_pdf], testval=2.5)
  File "/Users/brianzhu/anaconda/envs/py36/lib/python3.6/site-packages/pymc3/distributions/distribution.py", line 37, in __new__
    return model.Var(name, dist, data, total_size)
  File "/Users/brianzhu/anaconda/envs/py36/lib/python3.6/site-packages/pymc3/model.py", line 752, in Var
    total_size=total_size, model=self)
  File "/Users/brianzhu/anaconda/envs/py36/lib/python3.6/site-packages/pymc3/model.py", line 1130, in __init__
    self.logp_elemwiset = distribution.logp(self)
  File "/Users/brianzhu/anaconda/envs/py36/lib/python3.6/site-packages/pymc3/distributions/mixture.py", line 141, in logp
    return bound(logsumexp(tt.log(w) + self._comp_logp(value), axis=-1).sum(),
  File "/Users/brianzhu/anaconda/envs/py36/lib/python3.6/site-packages/pymc3/distributions/mixture.py", line 113, in _comp_logp
    axis=1)
  File "/Users/brianzhu/anaconda/envs/py36/lib/python3.6/site-packages/theano/tensor/basic.py", line 4709, in stack
    return join(axis, *[shape_padaxis(t, axis) for t in tensors])
  File "/Users/brianzhu/anaconda/envs/py36/lib/python3.6/site-packages/theano/tensor/basic.py", line 4709, in <listcomp>
    return join(axis, *[shape_padaxis(t, axis) for t in tensors])
  File "/Users/brianzhu/anaconda/envs/py36/lib/python3.6/site-packages/theano/tensor/basic.py", line 4601, in shape_padaxis
    raise IndexError(msg)
IndexError: axis 1 is out of bounds [-1, 1)

Thanks!

2

It should by removing the shape in the first mixture component:

with pm.Model() as model:
    # mixture component 1
    mu = pm.Normal("mu0", mu=5., sd=1.)
    sd = pm.HalfCauchy("sd0", beta=1.)
    gaus = pm.Normal('Norm', mu=mu, sd=sd)
    # mixture component 2
    interp_pdf = pm.Interpolated(
        'interp', x_points=x, pdf_points=y, testval=2.5)
    # mixture weight
    w = pm.Dirichlet('w', a=np.array([1., 1.]))
    # liklihood
    like = pm.Mixture(
        'like',
        w=w,
        comp_dists=[gaus.distribution, interp_pdf.distribution],
        testval=2.5)