Fitting model for two normal distributions in PyMC 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.
 A: A couple of points, related to the discussion above:


*

*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.

*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.
A: 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} )

