如何在PyMC中获取后验分布的参数？

Question

我有一个用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.)

看起来运行得不错，但我想从后验分布中提取一些参数。我该怎么从中获取后验的p，以及从Pi中获取alpha和beta呢？

Answer 1

你已经很接近了。如果你稍微调整一下，把Pi和Tj这两个对象放到你的函数外面，这样你就可以直接访问来自（近似）后验分布的MCMC样本了：

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)

然后你可以使用.trace()和.stats()方法来检查Ti和Pj的边际后验分布：

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}