Search code examples
pythonpymcpymc3

Defining stochastic and deterministic variables with pymc3

I am trying to use write my own stochastic and deterministic variables with pymc3, but old published recipe for pymc2.3 explained how we can parametrize our variables no longer works. For instance I tried to use this direct approach and it failed:

def x_logp(value, x_l, x_h):
    if ((value>x_h) or (value<x_l)):
        return -np.inf
    else:
        return -np.log(x_h-x_l+1)
def x_rand(x_l,x_h):
    return np.round((x_h-x_l)*np.random.random_sample())+x_l

Xpos=pm.stochastic(logp=x_logp,
                   doc="X position of halo center ",
                   observed=False, 
                   trace=True,
                   name='Xpos',
                   random=x_rand,
                   value=25.32,
                   parents={'x_l':0,'x_h'=500},
                   dtype=float64,
                   plot=True,
                   verbose=0)

I got the following error message:

ERROR: AttributeError: 'module' object has no attribute 'Stochastic' [unknown]
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
AttributeError: 'module' object has no attribute 'Stochastic'

I am wondering how I could define my own priors or likelihood in pymc3 without for instance, using decorators and the available pymc distributions?

Solution

  • There are essentially two ways to add custom densities:

    1. Theano expressions (can use gradient based samplers)

      You can use DensityDist for this, e.g.: https://github.com/pymc-devs/pymc/blob/master/pymc/examples/custom_dists.py

    2. Blackbox python functions (only non-gradient samplers like Metropolis, or Slice)

      Theano has a decorator you can use like this:

    
@theano.compile.ops.as_op(itypes=[t.lscalar, t.dscalar, t.dscalar],otypes=[t.dvector])
def rate(switchpoint,early_mean, late_mean):
    ''' Concatenate Poisson means '''
    out = empty(years)
    out[:switchpoint] = early_mean
    out[switchpoint:] = late_mean
    return out

    Taken from this example: https://github.com/pymc-devs/pymc/blob/master/pymc/examples/disaster_model_arbitrary_determinisitc.py

    Deterministics can either be done directly by combining the random variables, or, if you want them to show up in the trace using e.g. pm.Determinstic('sum', alpha + beta).