using decorators to define models in PyMC

Question

below is one way to define two stochastic Bernoulli random variables, one depending on the other with decorators. the model is meant to be:

p(A) = 0.5
p(B=True|A=True) = 0.75
p(B=True|A=False) = 0.05

using decorators in PyMC it is:

import pymc
from pymc import DiscreteUniform, Exponential, deterministic, Poisson, Uniform
import numpy as np

def make_model():
    @pymc.stochastic(dtype=bool)
    def A(value=False):
        def logp(value):
            return pymc.bernoulli_like(value, 0.5)

        def random():
            return bool(np.random.binomial(1, 10**logp(value)))


    @pymc.stochastic(dtype=bool)
    def B(value=False, A=False):
        def logp(value, A):
            if A:
                logp = pymc.bernoulli_like(value, 0.75)
            else:
                logp = pymc.bernoulli_like(value, 0.05)
            return logp

        def random(A):
            return bool(np.random.binomial(1, 10**logp(True, A)))

    return locals()

test = pymc.Model(make_model())

is this the most correct and compact way to do it? could some code be saved by defining the nodes explicitly as pymc.Bernoulli variables, rather than with the stochastic decorator?

the above code seems redundant having to always define how to sample from a bernoulli RV, but maybe this is unavoidable?

finally is it necessary to invoke sampling functions from numpy like np.random.binomial or are there sampling functions in PyMC?

score 1 · Answer 1 · answered Mar 03 '15 at 04:55

I do not think that this is correct. For example, if I sample from your model with the following code, I find the mean of the samples for A is not 0.5:

In [20]: pymc.MCMC(test).sample(10000)
 [-----------------100%-----------------] 10000 of 10000 complete in 0.6 sec
In [22]: test.A.trace().mean()
Out[22]:
0.19670000000000001

Here is a compact way to achieve your model:

import pymc as pm

def make_model():
    A = pm.Bernoulli('A', .5)
    B_given_A_true = pm.Bernoulli('B_given_A_true', .75)
    B_given_A_false = pm.Bernoulli('B_given_A_false', .05)

    @deterministic
    def B(A=A, B_given_A_true=B_given_A_true,
          B_given_A_false=B_given_A_false):
        if A:
            return B_given_A_true
        else:
            return B_given_A_false

    return locals()

test = pymc.Model(make_model())

why isn't B a stochastic node if its value is determined by a stochastic node A? — lgd, Mar 07 '15 at 16:36
This language is confusing, `B` is deterministic because it is a deterministic function of stochastics. See [here for more detail](http://pymc-devs.github.io/pymc/tutorial.html#two-types-of-variables). — Abraham D Flaxman, Mar 08 '15 at 22:31

using decorators to define models in PyMC

1 Answers1