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])

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)


{'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]?


The reason is that you have observed all 1000 draws from a single urn. This gives overwhelming probability that the assignment=0. But since you have observed only one assignment, you can't determine the value of urn_0_prob. That's like trying to estimate the bias of a coin from observing a single flip.


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.