# balls-and-urns in PyMC

In trying to wrap my head around Bayesian statistics and PyMC I've been trying to program up a balls-and-urns example. Sadly, I'm failing miserably.

Assume that I'm drawing balls from an urn and know that the composition of the urn is either 80 black/20 red balls or 20 black/80 red balls. I have a uniform prior over the probability that I'm drawing from the first urn and now get to draw repeatedly and with replacement from the urn to figure it out. I wrote this up as a PyMC model, like so:

import pymc as pymc
import numpy as np

# we're actually drawing from the first urn
balls = np.random.binomial(1,0.8,size=1000)

black_probs = [0.8,0.2]

urn_0_prob = pymc.Uniform('urn_0_prob',lower=0,upper=1)
assignment = pymc.Categorical("assignment", p=[urn_0_prob,1-urn_0_prob])

@pymc.deterministic
def black_prob(assignment=assignment):
return black_probs[assignment]

observations = pymc.Bernoulli("observations", p=black_prob, value=balls, observed=True)

model = pymc.Model([urn_0_prob,assignment,black_prob,observations])


This looks sensible to me but I'm obviously doing something wrong because this doesn't work:

mcmc = pymc.MCMC(model)
mcmc.sample(10000,0)

mcmc.stats()

{'assignment': {'95% HPD interval': array([ 0.,  0.]),
'mc error': 0.0,
'mean': 0.0,
'n': 10000,
'quantiles': {2.5: 0.0, 25: 0.0, 50: 0.0, 75: 0.0, 97.5: 0.0},
'standard deviation': 0.0},
'black_prob': {'95% HPD interval': array([ 0.8,  0.8]),
'mc error': 1.1102230246251566e-17,
'mean': 0.80000000000000016,
'n': 10000,
'quantiles': {2.5: 0.80000000000000004,
25: 0.80000000000000004,
50: 0.80000000000000004,
75: 0.80000000000000004,
97.5: 0.80000000000000004},
'standard deviation': 1.1102230246251565e-16},
'urn_0_prob': {'95% HPD interval': array([ 0.2186279 ,  0.99966969]),
'mc error': 0.0062019203981026425,
'mean': 0.6646822946585289,
'n': 10000,
'quantiles': {2.5: 0.15802813648110903,
25: 0.48846888794346444,
50: 0.70758248557324599,
75: 0.87243103830208035,
97.5: 0.98900646106466972},
'standard deviation': 0.24170398816390645}}


Why is the assignment always 0 when the probability that governs the assignment varies pretty widely in [0,1]?