Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Mixture Modeling example not working in pymc3 #443

Closed
zaxtax opened this issue Jan 3, 2014 · 31 comments

Comments

Projects
None yet
7 participants
@zaxtax
Copy link
Contributor

commented Jan 3, 2014

I am trying to port the mixture model to pymc3, but running into problem

https://gist.github.com/zaxtax/7968670

I get the following error as present

File "finite_mixture.py", line 21, in
shape=ndata)
File "/usr/local/lib/python2.7/dist-packages/pymc-3.0-py2.7.egg/pymc/distributions/distribution.py", line 19, in new
return model.Var(name, dist, data)
File "/usr/local/lib/python2.7/dist-packages/pymc-3.0-py2.7.egg/pymc/model.py", line 144, in Var
var = FreeRV(name=name, distribution=dist, model=self)
File "/usr/local/lib/python2.7/dist-packages/pymc-3.0-py2.7.egg/pymc/model.py", line 328, in init
distribution.shape, distribution.dtype) * distribution.default()
ValueError: operands could not be broadcast together with shapes (500) (3)

If I replace Multinomial with Categorical I get

Traceback (most recent call last):
File "finite_mixture.py", line 27, in
tr = sample(3000, step=Metropolis(model.vars))
File "/usr/local/lib/python2.7/dist-packages/pymc-3.0-py2.7.egg/pymc/sample.py", line 48, in sample
random_seed=random_seed)):
File "/usr/local/lib/python2.7/dist-packages/pymc-3.0-py2.7.egg/pymc/sample.py", line 118, in iter_sample
point = step.step(point)
File "/usr/local/lib/python2.7/dist-packages/pymc-3.0-py2.7.egg/pymc/step_methods/arraystep.py", line 24, in step
apoint = self.astep(bij.map(point), inputs)
File "/usr/local/lib/python2.7/dist-packages/pymc-3.0-py2.7.egg/pymc/step_methods/metropolis.py", line 105, in astep
q_new = metrop_select(logp(q) - logp(q0), q, q0)
File "/usr/local/lib/python2.7/dist-packages/pymc-3.0-py2.7.egg/pymc/blocking.py", line 119, in call
return self.fa(self.fb(x))
File "/usr/local/lib/python2.7/dist-packages/pymc-3.0-py2.7.egg/pymc/model.py", line 290, in call
return self.f(
*state)
File "/usr/local/lib/python2.7/dist-packages/theano/compile/function_module.py", line 588, in call
self.fn.thunks[self.fn.position_of_error])
File "/usr/local/lib/python2.7/dist-packages/theano/compile/function_module.py", line 579, in call
outputs = self.fn()
IndexError: index 3 is out of bounds for size 3
Apply node that caused the error: AdvancedSubtensor1(dd, category)
Inputs shapes: [(3,), (500,)]
Inputs strides: [(8,), (8,)]
Inputs types: [TensorType(float64, vector), TensorType(int64, vector)]
Use the Theano flag 'exception_verbosity=high' for a debugprint of this apply node.

@twiecki

This comment has been minimized.

Copy link
Member

commented Jan 3, 2014

I played around a little and it seems to at least sample with this (some other unnecessary changes):

import numpy as np
import pymc as pm
from pymc import Model, Gamma, Normal, Dirichlet
from pymc import Categorical
from pymc import sample, Metropolis

k = 3
ndata = 500

v = np.random.randint(0, k, ndata)
data = ((v == 0)*(50 + np.random.randn(ndata))
        + (v == 1)*(-50 + np.random.randn(ndata))
        + (v == 2)*np.random.randn(ndata))

with Model() as model:
    dd = Dirichlet('dd', a=np.array([1., 1., 1.]), shape=k)
    sd = pm.Uniform('precs', lower=0, upper=20, shape=k)
    means = Normal('means', [-50, 0, 50], sd=15, shape=k)
    category = Categorical('category',
                           p=dd,
                           shape=ndata)

    points = Normal('obs',
                     means[category],
                     sd=sd[category],
                     observed=data)
    tr = sample(3000, step=Metropolis())

However, category is always set to 0 in the trace.

@bjedwards

This comment has been minimized.

Copy link
Contributor

commented Jan 4, 2014

I think you need to use ElemwiseCategoricalStep for the category, however, I can't get it to work.

import numpy as np
import pymc as pm
from pymc import Model, Gamma, Normal, Dirichlet
from pymc import Categorical
from pymc import sample, Metropolis, ElemwiseCategoricalStep

k = 3
ndata = 500

v = np.random.randint(0, k, ndata)
data = ((v == 0)*(50 + np.random.randn(ndata))
        + (v == 1)*(-50 + np.random.randn(ndata))
        + (v == 2)*np.random.randn(ndata))

with Model() as model:
    dd = Dirichlet('dd', a=np.array([1., 1., 1.]), shape=k)
    sd = pm.Uniform('precs', lower=0, upper=20, shape=k)
    means = Normal('means', [-50, 0, 50], sd=15, shape=k)
    category = Categorical('category',
                           p=dd,
                           shape=ndata)

    points = Normal('obs',
                     means[category],
                     sd=sd[category],
                     observed=data)
    step1 = Metropolis(vars = [dd,sd,means])
    step2 = ElemWiseCategoricalStep(var=category,vals=[0,1,2])
    tr = sample(3000, step=[step1,step2])

I get an weird attribute error I haven't run down yet

---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
<ipython-input-60-ea16950056f3> in <module>()
      1 with model:
      2     step1 = pm.Metropolis(vars=[dd,sd,means])
----> 3     step2 = pm.ElemwiseCategoricalStep(var=category,values=[0,1,2])

/nfs/adaptive/bedwards/Documents/Dev/pymc/pymc/step_methods/gibbs.pyc in __init__(self, var, values, model)
     29         self.var = var
     30 
---> 31         ArrayStep.__init__(self, [var], [elemwise_logp(model, var)])
     32 
     33     def astep(self, q, logp):

/nfs/adaptive/bedwards/Documents/Dev/pymc/pymc/step_methods/gibbs.pyc in elemwise_logp(model, var)
     37 
     38 def elemwise_logp(model, var):
---> 39     terms = [term for term in model.factors if var in inputs([term])]
     40     return add(*terms)
     41 

AttributeError: 'Model' object has no attribute 'factors'

This is wierd, as I thought model used to have factors. Was this something that escaped the Pyclasses2 merge?

@twiecki

This comment has been minimized.

Copy link
Member

commented Jan 4, 2014

Seems like elemwise_logp is broken. model.factors should be model.free_RVs (I think) and not sure what add should be.

@jsalvatier

This comment has been minimized.

Copy link
Member

commented Jan 5, 2014

To clarify, the original problem was the missing shape=ndata on Categorical.

The second problem with elemwise_logp @bjedwards found should be fixed now in pymc-devs/pymc@f8661d2.

Fixing that makes bjedwards example work for me.

There's an additional problem with traceplot that I'm trying to fix.

Also the issue with metropolis in https://github.com/pymc-devs/pymc/issues/358 seems to be present here as well :(

@twiecki

This comment has been minimized.

Copy link
Member

commented Apr 29, 2015

This is the updated code:

import numpy as np
import pymc3 as pm
from pymc3 import Model, Gamma, Normal, Dirichlet
from pymc3 import Categorical
from pymc3 import sample, Metropolis, ElemwiseCategoricalStep

k = 3
ndata = 500

v = np.random.randint(0, k, ndata)
data = ((v == 0)*(50 + np.random.randn(ndata))
        + (v == 1)*(-50 + np.random.randn(ndata))
        + (v == 2)*np.random.randn(ndata))

a = pm.constant(np.array([1., 1., 1.]))

with Model() as model:
    p, p_m1 = model.TransformedVar(
        'p', 
        Dirichlet.dist(a, shape=k, testval=[0.34,.33,.33]),
        pm.simplextransform)

    sd = pm.Uniform('sd', lower=0, upper=20, shape=k)
    means = Normal('means', [-50, 0, 50], sd=15, shape=k)

    category = Categorical('category',
                           p=p,
                           shape=ndata)

    points = Normal('obs',
                     means[category],
                     sd=sd[category],
                     observed=data)

    step1 = Metropolis(vars = [p_m1, sd, means])
    step2 = ElemwiseCategoricalStep(var=category, values=[0,1,2])
    #start = pm.find_MAP()
    tr = sample(300, step=[step1,step2])

Unfortunately, the dirichlet always samples nans and I'm not sure why.

@zaxtax

This comment has been minimized.

Copy link
Contributor Author

commented Apr 29, 2015

Interesting. Let me poke around.
On 29 Apr 2015 16:04, "Thomas Wiecki" notifications@github.com wrote:

This is the updated code:

import numpy as npimport pymc3 as pmfrom pymc3 import Model, Gamma, Normal, Dirichletfrom pymc3 import Categoricalfrom pymc3 import sample, Metropolis, ElemwiseCategoricalStep

k = 3
ndata = 500

v = np.random.randint(0, k, ndata)
data = ((v == 0)(50 + np.random.randn(ndata))
+ (v == 1)
(-50 + np.random.randn(ndata))
+ (v == 2)*np.random.randn(ndata))

a = pm.constant(np.array([1, 1, 1]))
with Model() as model:
p, p_m1 = model.TransformedVar(
'p',
Dirichlet.dist(np.array([1, 1, 1]), shape=k),
pm.simplextransform)

sd = pm.Uniform('sd', lower=0, upper=20, shape=k)
means = Normal('means', [-50, 0, 50], sd=15, shape=k)

category = Categorical('category',
                       p=p,
                       shape=ndata)

points = Normal('obs',
                 means[category],
                 sd=sd[category],
                 observed=data)

step1 = Metropolis(vars = [p_m1, sd, means])
step2 = ElemwiseCategoricalStep(var=category, values=[0,1,2])
#start = pm.find_MAP()
tr = sample(300, step=[step1,step2])

Unfortunately, the dirichlet always samples nans and I'm not sure why.


Reply to this email directly or view it on GitHub
#443 (comment).

@hgbrian

This comment has been minimized.

Copy link
Contributor

commented Apr 30, 2015

Just to note that replacing

Dirichlet.dist(a, shape=k)

with

Dirichlet.dist(a, shape=k, testval=[0.34,.33,.33])

appears to address the nans problem.

@twiecki

This comment has been minimized.

Copy link
Member

commented Apr 30, 2015

@hgbrian Thanks, that indeed solves that one problem. Still, it seems like the categorical is never changing. I suspect `ElemwiseCategoricalStep' isn't working correctly.

@zaxtax

This comment has been minimized.

Copy link
Contributor Author

commented May 2, 2015

Weird,

I get the following:

Cannot compute test value: input 0 (<TensorType(float32, matrix)>) of Op
Dot22(<TensorType(float32, matrix)>, <TensorType(float32, matrix)>) missing
default value
zsh: floating point exception (core dumped)  ipython

When I run this under gdb I get:

gdb python
run finite_mixture.py

Program received signal SIGFPE, Arithmetic exception.
0x00007fffd112ecc2 in run (this=0x15874f0) at /home/zv/.theano/compiledir_Linux-3.13--generic-      x86_64-with-Ubuntu-14.04-trusty-x86_64-2.7.6-64/tmpOIM6Kn/mod.cpp:461
461         V1_i = ((-V3_i) / (-V5_i));

On Thu, Apr 30, 2015 at 3:50 PM, Thomas Wiecki notifications@github.com
wrote:

@hgbrian https://github.com/hgbrian Thanks, that indeed solves that one
problem. Still, it seems like the categorical is never changing. I suspect
`ElemwiseCategoricalStep' isn't working correctly.


Reply to this email directly or view it on GitHub
#443 (comment).

@twiecki

This comment has been minimized.

Copy link
Member

commented May 2, 2015

Is this with current master theano?

@zaxtax

This comment has been minimized.

Copy link
Contributor Author

commented May 2, 2015

Yep. Also I get the floating point exception with the theano in pypi.
On 2 May 2015 04:37, "Thomas Wiecki" notifications@github.com wrote:

Is this with current master theano?


Reply to this email directly or view it on GitHub
#443 (comment).

@zaxtax

This comment has been minimized.

Copy link
Contributor Author

commented May 2, 2015

Update: I don't get any errors on python 3.4 or on OS X, so I'm going to assume I just borked my python 2.7 stack.

@zaxtax

This comment has been minimized.

Copy link
Contributor Author

commented May 2, 2015

I can also make the nans go away also by changing in pymc3/distributions/multivariate.py:

    def __init__(self, a, *args, **kwargs):
        super(Dirichlet, self).__init__(*args, **kwargs)
        self.a = a
        self.k = a.shape[0]
        self.mean = a / sum(a)

        self.mode = switch(all(a > 1),
                           (a - 1) / sum(a - 1),
                           nan)

to:

def __init__(self, a, *args, **kwargs):
        super(Dirichlet, self).__init__(*args, **kwargs)
        self.a = a
        self.k = a.shape[0]
        self.mean = a / sum(a)

        self.mode = self.mean
@twiecki

This comment has been minimized.

Copy link
Member

commented May 3, 2015

That indeed looks wrong and explains what we've been seeing.
On May 3, 2015 1:53 AM, "zaxtax" notifications@github.com wrote:

I can also make the nans go away also by changing in
pymc3/distributions/multivariate.py:

def __init__(self, a, *args, **kwargs):
    super(Dirichlet, self).__init__(*args, **kwargs)
    self.a = a
    self.k = a.shape[0]
    self.mean = a / sum(a)

    self.mode = switch(all(a > 1),
                       (a - 1) / sum(a - 1),
                       nan)

to:

def init(self, a, _args, *_kwargs):
super(Dirichlet, self).init(_args, *_kwargs)
self.a = a
self.k = a.shape[0]
self.mean = a / sum(a)

    self.mode = self.mean


Reply to this email directly or view it on GitHub
#443 (comment).

@twiecki

This comment has been minimized.

Copy link
Member

commented May 3, 2015

It looks mathematically correct as it was before (http://en.wikipedia.org/wiki/Dirichlet_distribution). I think the problem is rather that when setting a starting value we should test if the mode is nan and if it is, try the mean.

@zaxtax

This comment has been minimized.

Copy link
Contributor Author

commented May 3, 2015

I agree. So maybe we should change get_test_val in class Distribution. Is that the right place to make the change?

@twiecki

This comment has been minimized.

Copy link
Member

commented May 3, 2015

Yep, exactly. That should check for isnull I guess.

@zaxtax

This comment has been minimized.

Copy link
Contributor Author

commented May 3, 2015

The nan is buried in Theano. Is there a nice way to propagate that up?

@jsalvatier

This comment has been minimized.

Copy link
Member

commented May 4, 2015

Yup, the value is in val.tag.test_value.

On Sun, May 3, 2015 at 3:49 PM, zaxtax notifications@github.com wrote:

The nan is buried in Theano. Is there a nice way to propagate that up?


Reply to this email directly or view it on GitHub
#443 (comment).

@zaxtax

This comment has been minimized.

Copy link
Contributor Author

commented May 4, 2015

Awesome!

@twiecki

This comment has been minimized.

Copy link
Member

commented May 26, 2015

Unfortunately, even after fixing the nan issue, the model does not seem to work as no Metropolis step is ever accepted.

@twiecki

This comment has been minimized.

Copy link
Member

commented Jun 5, 2015

I think the issue is that we do blocked updating here. Especially for the categorical this will result in 0 probability of ever accepting a jump -- the clusters are so far away that wrongly proposing a single data points' cluster will screw things up. I think we need Gibbs sampling here.

@aflaxman

This comment has been minimized.

Copy link
Contributor

commented Jun 7, 2015

I took a look at this today, and came up with an approach that seems to work. I haven't read through the whole discussion here, but I thought you might find my notebook on this useful: http://nbviewer.ipython.org/gist/aflaxman/64f22d07256f67396d3a

It seems like PyMC3 is really coming along. Cool!

@twiecki

This comment has been minimized.

Copy link
Member

commented Jun 8, 2015

@aflaxman This is really neat! You solved quite a few problems, specifically:

  • Gibbs sampling for the cluster assignments (outlined as an issue above)
  • Breaking symmetry in the posterior by enforcing index ordering
  • Ensure all clusters are non-empty.

I suppose we should start thinking if and how any of this could be productized for pymc3.

@zaxtax

This comment has been minimized.

Copy link
Contributor Author

commented Jun 8, 2015

I think some form of these step methods would be great to have in pymc3. As
another note, it might help to document the step methods already there and
some of their strengths and weaknesses.
On Jun 8, 2015 4:46 AM, "Thomas Wiecki" notifications@github.com wrote:

@aflaxman https://github.com/aflaxman This is really neat! You solved
quite a few problems, specifically:

  • Gibbs sampling for the cluster assignments (outlined as an issue
    above)
  • Breaking symmetry in the posterior by enforcing index ordering
  • Ensure all clusters are non-empty.

I suppose we should start thinking if and how any of this could be
productized for pymc3.


Reply to this email directly or view it on GitHub
#443 (comment).

@jsalvatier

This comment has been minimized.

Copy link
Member

commented Jun 9, 2015

@thomas, I agree. We should think about how to productize this.

On Mon, Jun 8, 2015 at 3:09 AM, zaxtax notifications@github.com wrote:

I think some form of these step methods would be great to have in pymc3. As
another note, it might help to document the step methods already there and
some of their strengths and weaknesses.
On Jun 8, 2015 4:46 AM, "Thomas Wiecki" notifications@github.com wrote:

@aflaxman https://github.com/aflaxman This is really neat! You solved
quite a few problems, specifically:

  • Gibbs sampling for the cluster assignments (outlined as an issue
    above)
  • Breaking symmetry in the posterior by enforcing index ordering
  • Ensure all clusters are non-empty.

I suppose we should start thinking if and how any of this could be
productized for pymc3.


Reply to this email directly or view it on GitHub
#443 (comment).


Reply to this email directly or view it on GitHub
#443 (comment).

@twiecki

This comment has been minimized.

Copy link
Member

commented Aug 17, 2015

Looking back at this example, it seems to work just as well with:

# fit model
with model:
    step1 = pm.Metropolis(vars=[p, sd, means])
    step2 = pm.ElemwiseCategoricalStep(var=category, values=[0, 1, 2])
    tr = pm.sample(10000, step=[step1, step2])

That isn't to say that we shouldn't add Gibbs sampling but it seems like we already have functioning categorical Gibbs sampling.

@twiecki

This comment has been minimized.

Copy link
Member

commented Aug 17, 2015

I added Abe's NB with a modification to use the pymc3 categorical sampler to the examples folder: https://github.com/pymc-devs/pymc3/blob/master/pymc3/examples/gaussian_mixture_model.ipynb

Hope that's OK @aflaxman?

@twiecki

This comment has been minimized.

Copy link
Member

commented Aug 17, 2015

Come to think of it, the automatic transforms (and the stick breaking transform for the dirichlet) might have helped as well. Anyway, seems like this is working now, so closing.

@twiecki twiecki closed this Aug 17, 2015

@parashardhapola

This comment has been minimized.

Copy link

commented Jul 2, 2016

Can't find the mixture model notebook in the examples folder.

@twiecki

This comment has been minimized.

Copy link
Member

commented Jul 2, 2016

https://github.com/pymc-devs/pymc3/blob/master/docs/source/notebooks/gaussian-mixture-model-advi.ipynb

On Sat, Jul 2, 2016 at 7:41 AM, Parashar notifications@github.com wrote:

Can't find the mixture model notebook in the examples folder
https://github.com/pymc-devs/pymc3/tree/master/pymc3/examples.


You are receiving this because you modified the open/close state.
Reply to this email directly, view it on GitHub
#443 (comment),
or mute the thread
https://github.com/notifications/unsubscribe/AApJmCRFJVAZNNb0FwCQN0Fc7o2q2UPFks5qRfoZgaJpZM4BXWJt
.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.