Closed zaxtax closed 9 years ago
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.
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?
Seems like elemwise_logp
is broken. model.factors
should be model.free_RVs (I think) and not sure what add
should be.
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 https://github.com/pymc-devs/pymc/commit/f8661d2ef425c3c86d2d5a224e3ad73db1b8d703.
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 :(
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.
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 https://github.com/pymc-devs/pymc3/issues/443#issuecomment-97566035.
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.
@hgbrian Thanks, that indeed solves that one problem. Still, it seems like the categorical is never changing. I suspect `ElemwiseCategoricalStep' isn't working correctly.
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 https://github.com/pymc-devs/pymc3/issues/443#issuecomment-97947331.
Is this with current master theano
?
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 https://github.com/pymc-devs/pymc3/issues/443#issuecomment-98331054.
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.
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
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 https://github.com/pymc-devs/pymc3/issues/443#issuecomment-98418546.
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.
I agree. So maybe we should change get_test_val in class Distribution. Is that the right place to make the change?
Yep, exactly. That should check for isnull
I guess.
The nan
is buried in Theano. Is there a nice way to propagate that up?
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 https://github.com/pymc-devs/pymc3/issues/443#issuecomment-98553242.
Awesome!
Unfortunately, even after fixing the nan issue, the model does not seem to work as no Metropolis step is ever accepted.
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.
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!
@aflaxman This is really neat! You solved quite a few problems, specifically:
I suppose we should start thinking if and how any of this could be productized for pymc3
.
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 https://github.com/pymc-devs/pymc3/issues/443#issuecomment-109911613.
@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 https://github.com/pymc-devs/pymc3/issues/443#issuecomment-109911613.
— Reply to this email directly or view it on GitHub https://github.com/pymc-devs/pymc3/issues/443#issuecomment-109941134.
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.
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?
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.
Can't find the mixture model notebook in the examples folder.
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 https://github.com/pymc-devs/pymc3/issues/443#issuecomment-230085260, or mute the thread https://github.com/notifications/unsubscribe/AApJmCRFJVAZNNb0FwCQN0Fc7o2q2UPFks5qRfoZgaJpZM4BXWJt .
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 = metropselect(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.