I am currently experimenting with PYMC and I am trying out a simple example so that I start learning how things work (I am also a Python beginner, but an experienced machine learner).

I have set myself the following very simple model selection task:

Imagine we have a dataset that has been entirely generated either from model 0 or model 1 (not a mixture).

  • Model 0 is a normal distribution with mean m0 and standard deviation s0.

  • Likewise, model 1 is also a normal distribution with mean m1 and standard deviation s1.

  • We assume that we know all parameters m0,s0,m1,s1.

  • We postulate a hidden binary random variable V which indicates the model that really generated the dataset.

What we would like to do is infer the posterior of V.

I have made an attempt in the code below, but unfortunately after trying for a couple of days, I capitulate and throw myself at your feet for assistance.

The example is of course very simple and may appear contrived. What I would like to do in the future, is use PYMC for model selection. Hence, how one correctly handles a hidden indicator variable like V is important to me. Unfortunately, in my code I keep getting weird results for the posterior of V.

I hope the code below makes sense. I thank you all for your time.

from pymc  import *
import matplotlib 
from pylab import hist, show, plot, title, figure
import numpy as np
from numpy import *

# Generate some data

# Define means for gaussians
mu0 =  5.0
mu1 = -5.0

# Define standard deviations for gaussians 
s0 = 1.0
s1 = 1.0

# This variable chooses from which model we generate the observed data.
# It can be either 0 or 1
true_model = 0

# number of data items
numData = 100

if true_model==0:
    # Generate data from first model
    data = np.random.randn(numData,1)*s0 + mu0
elif true_model==1:
    # Generate data from second model
    data = np.random.randn(numData,1)*s1 + mu1

# Define variables
PR = Dirichlet("priorOnModels", theta=(0.5,0.5)) # both models are equally likely
V  = Categorical('V', p=PR)

# Define model log-likelihood
def mymodel_loglikel(sw=V):

    if sw==0:        
        loglikel = distributions.normal_like(data, mu0, 1.0/s0**2)

    elif sw==1:        
        loglikel = distributions.normal_like(data, mu1, 1.0/s1**2)

    return loglikel

# Run sampler
simplemodel = Model([mymodel_loglikel, PR, V])
mc = MCMC(simplemodel)

print "expectation of V variable is %.3f" % np.mean(mc.trace('V')[:])


EDIT 2: Updated after Tom's suggestions

In response to a comment below, the weird behaviour I get concerning V, is that its mean posterior is always 0.0. Judging from the trace, it seems that V is not properly sampled. I suspect a minor programming error somewhere, but can't pinpoint it.

I attach below a figure of the trace plots of V.

enter image description here

  • $\begingroup$ Could you show the results you get for the posterior of $V$ and explain why you think they are weird. $\endgroup$ – Juho Kokkala Oct 29 '14 at 16:25

When defining a potential, the function should return the log-probability. So instead of

return np.exp(loglikel)

You want

return loglikel

Another thing to watch out for is the expression 1/s0**2. Make sure that s0 is float in this expression, or change it to 1.0/s0**2.

PyMC also has a bug in its choice of proposals for binary variables (documented in section 5.7.3), which I first noticed in another answer. The fix is to add the line:

mc.use_step_method(DiscreteMetropolis, V, proposal_distribution='Prior')
| cite | improve this answer | |
  • $\begingroup$ Thanks for taking a look at the code Tom. I will try out your suggestions and get back to you.Thanks, Nikos $\endgroup$ – ngiann Oct 30 '14 at 18:50
  • $\begingroup$ I tried out both suggestions and I updated the code above in the original post. However, I still get weird behaviour for for variable V. I have accordingly updated the figure in the original post above. The weird thing is that variable V seems to be somehow stuck? In the trace plot above, it seems to always receive the value 0 during sampling? I imagine that a minor programming error is at fault somewhere here... Thanks, N. $\endgroup$ – ngiann Oct 31 '14 at 9:40
  • $\begingroup$ Try using a different step method for V, e.g. mc.use_step_method(DiscreteMetropolis, V) before mc.sample $\endgroup$ – Abraham D Flaxman Oct 31 '14 at 14:48
  • $\begingroup$ Thanks for the suggestion, but unfortunately it didn't change a thing.... $\endgroup$ – ngiann Oct 31 '14 at 15:06
  • $\begingroup$ I don't understand your updated question. When true_model = 0, it estimates V=0. When true_model = 1, it estimates V=1. Isn't that what you want? $\endgroup$ – Tom Minka Oct 31 '14 at 15:16

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.