Simple example, mixture of gaussians

[timini]

Sorry for bugging you with a question that is probably a fault of my modelling please feel free to delete this issue. I posted a question on cross validated about porting a simple model from pymc2 to pymc3:

http://stats.stackexchange.com/questions/81650/converting-a-mixture-of-gaussians-to-pymc3

[fonnesbeck]

I will answer you question on CV.

[fonnesbeck]

Actually, this may be a bug so I am reopening.

[fonnesbeck]

It turns out that your model is not working, so nodes like p have a trace populated with 10000 values of 0.5. So, when SciPy tries to calculate a KDE of the trace, it chokes. We probably need to add a more informative error message there.

[fonnesbeck]

One of the issues is that ber is specified to be a scalar (the default). It looks like it needs to be the same length as the data, no?

[fonnesbeck]

The other issue is that you are trying to use Metropolis to update a Bernoulli random variable. You should use BinaryMetropolis. I have taken the liberty of updating your model, including a few improvements on the specification. This one should work:

with pm.Model() as model:
    #priors
    p = pm.Uniform( "p", 0 , 1) #this is the fraction that come from mean1 vs mean2

    ber = pm.Bernoulli( "ber", p = p, shape=len(data)) # produces 1 with proportion p.

    sigma = pm.Uniform('sigma', 0, 100)
    precision = sigma**-2

    mean = pm.Normal( "mean", 0, 0.01, shape=2 )

    mu = pm.Deterministic('mu', mean[ber])

    process = pm.Normal('process', mu=mu, tau=precision, observed=data)

with model:
    step1 = pm.Metropolis([p, sigma, mean])
    step2 = pm.BinaryMetropolis([ber])
    trace = pm.sample(10000, [step1, step2])

[fonnesbeck]

BTW, the best place to put usage questions regarding PyMC is on StackOverflow; CrossValidated tends not to like software questions. Be sure to add the pymc tag.