How do I get parameters from a posterior distribution in PyMC?

Question

I have the following program written in PyMC:

import pymc
from pymc.Matplot import plot as mcplot

def testit( passed, test_p = 0.8, alpha = 5, beta = 2):
    Pi = pymc.Beta( 'Pi', alpha=alpha, beta=beta)
    Tj = pymc.Bernoulli( 'Tj', p=test_p)

    @pymc.deterministic
    def flipper( Pi=Pi, Tj=Tj):
            return Pi if Tj else (1-Pi)
            # Pij = Pi if Tj else (1-Pi)
            # return pymc.Bernoulli( 'Rij', Pij)

    Rij = pymc.Bernoulli( 'Rij', p=flipper, value=passed, observed=True)

    model = pymc.MCMC( [ Pi, Tj, flipper, Rij])
    model.sample(iter=10000, burn=1000, thin=10)

    mcplot(model)

testit( 1.)

It appears to be working properly, but I'd like to extract parameters from the posterior distributions. How can I get the posterior p from Tj and alpha/beta from Pi?

Answer 1

You are very close. If you refactor a little bit so that you have the Pi and Tj objects outside of your function you can get access to the MCMC samples from the (approximate) posterior distribution directly:

import pymc

def testit(passed, test_p = 0.8, alpha = 5, beta = 2):
    Pi = pymc.Beta( 'Pi', alpha=alpha, beta=beta)
    Tj = pymc.Bernoulli( 'Tj', p=test_p)

    @pymc.deterministic
    def flipper( Pi=Pi, Tj=Tj):
            return Pi if Tj else (1-Pi)
            # Pij = Pi if Tj else (1-Pi)
            # return pymc.Bernoulli( 'Rij', Pij)

    Rij = pymc.Bernoulli( 'Rij', p=flipper, value=passed, observed=True)

    return locals()

vars = testit(1.)
model = pymc.MCMC(vars)
model.sample(iter=10000, burn=1000, thin=10)

Then you can examine the marginal posterior distribution of Ti and Pj with the .trace() and .stats() methods:

In [12]: model.Pi.stats()
Out[12]:
{'95% HPD interval': array([ 0.43942434,  0.9910729 ]),
 'mc error': 0.0054870077893956213,
 'mean': 0.7277823553617826,
 'n': 900,
 'quantiles': {2.5: 0.3853555534589701,
  25: 0.62928387568176036,
  50: 0.7453244339604943,
  75: 0.84835518829619661,
  97.5: 0.95826093368693854},
 'standard deviation': 0.15315966296243455}
In [13]: model.Tj.stats()
Out[13]:
{'95% HPD interval': array([ 0.,  1.]),
 'mc error': 0.011249691353790801,
 'mean': 0.89666666666666661,
 'n': 900,
 'quantiles': {2.5: 0.0, 25: 1.0, 50: 1.0, 75: 1.0, 97.5: 1.0},
 'standard deviation': 0.30439375084839554}