12
$\begingroup$

Since I'm a software engineer trying to learn more stats you'll have to forgive me before I even start, this is serious newb territory...

I've been learning PyMC and working through some really (really) simple examples. One problem I can't get to work (and can't find any related examples for) is fitting a model to data generated from two normal distributions.

Say I have 1000 values; 500 generated from a Normal(mean=100, stddev=20) and another 500 generated from a Normal(mean=200, stddev=20).

If I want to fit a model to them, ie determine the two means and the single standard deviation, using PyMC. I know it's something along the lines of ...

mean1 = Uniform('mean1', lower=0.0, upper=200.0)
mean2 = Uniform('mean2', lower=0.0, upper=200.0)
precision = Gamma('precision', alpha=0.1, beta=0.1)

data = read_data_from_file_or_whatever()

@deterministic(plot=False)
def mean(m1=mean1, m2=mean2):
    # but what goes here?

process = Normal('process', mu=mean, tau=precision, value=data, observed=True)

i.e., the generating process is Normal, but mu is one of two values. I just don't know how to represent the "decision" between whether a value comes from m1 or m2.

Perhaps I'm just completely taking the wrong approach to modeling this? Can anyone point me at an example? I can read BUGS and JAGS so anything is ok really.

Improve this question
$\endgroup$

2 Answers 2

Reset to default
12
$\begingroup$

Are you absolutely certain that half came from one distribution and the other half from the other? If not, we can model the proportion as a random variable (which is a very bayesian thing to do).

The following is what I would do, some tips are embedded.

from pymc import *

size = 10
p = Uniform( "p", 0 , 1) #this is the fraction that come from mean1 vs mean2

ber = Bernoulli( "ber", p = p, size = size) # produces 1 with proportion p.

precision = Gamma('precision', alpha=0.1, beta=0.1)

mean1 = Normal( "mean1", 0, 0.001 ) #better to use normals versus Uniforms (unless you are certain the value is  truncated at 0 and 200 
mean2 = Normal( "mean2", 0, 0.001 )

@deterministic
def mean( ber = ber, mean1 = mean1, mean2 = mean2):
    return ber*mean1 + (1-ber)*mean2


#generate some artificial data   
v = np.random.randint( 0, 2, size)
data = v*(10+ np.random.randn(size) ) + (1-v)*(-10 + np.random.randn(size ) )


obs = Normal( "obs", mean, precision, value = data, observed = True)

model = Model( {"p":p, "precision": precision, "mean1": mean1, "mean2":mean2, "obs":obs} )
Improve this answer
$\endgroup$
8
  • 2
    $\begingroup$ Shameless promotion: I just wrote a blog article about Bayes and pyMC literally 1 minute before you posted this, so I invite you to check it out. The Awesome Power of Bayes - Part 1 $\endgroup$
    – Cam.Davidson.Pilon
    Commented Dec 27, 2012 at 23:23
  • $\begingroup$ awesome! this approach to the mixing of the two means is exactly what i was trying to get my head around. $\endgroup$
    – mat kelcey
    Commented Dec 27, 2012 at 23:29
  • $\begingroup$ Not sure I fully understand the true modelling benefit of saying mean1 & mean2 are Normally distributed instead of Uniform (Same goes really for the precision to be honest, I've been using Gamma since "someone else did"). I've got a lot to learn :) $\endgroup$
    – mat kelcey
    Commented Dec 27, 2012 at 23:31
  • $\begingroup$ Using a Uniform, as in your original example, implies that you know with absolute certainty that the mean does not exceed some value. This is somewhat pathological. It is better to use a normal, as it allows all real numbers to be considered. $\endgroup$
    – Cam.Davidson.Pilon
    Commented Dec 27, 2012 at 23:34
  • 1
    $\begingroup$ The choice of gamma has a mathematical reason. The gamma is the conjugate prior of the precision, see table here $\endgroup$
    – Cam.Davidson.Pilon
    Commented Dec 27, 2012 at 23:36
7
$\begingroup$

A couple of points, related to the discussion above:

  1. The choice of diffuse normal vs. uniform is pretty academic unless (a) you are worried about conjugacy, in which case you would use the normal or (b) there is some reasonable chance that the true value could be outside the endpoints of the uniform. With PyMC, there is no reason to worry about conjugacy, unless you specifically want to use a Gibbs sampler.

  2. A gamma is actually not a great choice for an uninformative prior to a variance/precision parameter. It can end up being more informative that you think. A better choice is to put a uniform prior on the standard deviation, then transform it by an inverse square. See Gelman 2006 for details.

Improve this answer
$\endgroup$
6
  • 1
    $\begingroup$ ah fonnesbeck is one of the core developers of pymc! Can you show us an example of how to code point 2? $\endgroup$
    – Cam.Davidson.Pilon
    Commented Dec 28, 2012 at 16:29
  • $\begingroup$ thanks fonnesbeck and, yes please! to a quick eg of point 2 :) $\endgroup$
    – mat kelcey
    Commented Dec 29, 2012 at 4:41
  • 1
    $\begingroup$ in fact i'm guessing you mean something along the lines of ... gist.github.com/4404631 ? $\endgroup$
    – mat kelcey
    Commented Dec 29, 2012 at 4:51
  • $\begingroup$ Yes, exactly. You can do the transform a little more concisely: tau = std_dev**-2 $\endgroup$
    – fonnesbeck
    Commented Jan 3, 2013 at 18:40
  • $\begingroup$ what would be the right place to read about where this relation between precision and std_dev comes from? $\endgroup$
    – user979
    Commented May 22, 2013 at 9:37

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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