Defining stochastic and deterministic variables with pymc3

[Zahra Sheikh]

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'

or if I want to define an Exponential prior with the class that I have written:

class LogPrior(object):
    def eval(self, value):
        return 0.
    def __call__(self, value):
        return self.eval(value)
    def sample(self, n=None):
        """ Sample from this prior. The returned array axis=0 is the
            sample axis.

            Parameters
            ----------
            n : int (optional)
                Number of samples to draw
        """
        raise ValueError("Cannot sample from a LogPrior object.")
    def __str__(self):
        return "<LogPrior>"
    def __repr__(self):
        return self.__str__()

class ExponentialPrior(LogPrior):
    def __init__(self, lam, *args, **kwargs):
        super(ExponentialPrior, self).__init__(*args, **kwargs)
        self.lam = lam
        self.mean = 1. / lam
        self.median = self.mean * log(2)
        self.mode = 0
        self.variance = lam ** -2
    def logp(self, value, limits=None):
        if limits:
           lower,upper=limits
           """Log of lognormal prior probability with hard limits."""
           if value >= lower and value <= upper:
          return -log(self.lam)+self.lam*value
           else:
              return -inf
        else:
              """Log of normal prior probability."""
              return -log(self.lam)+self.lam*value
    def cdf(self, value):
        """Cumulative distribution function lognormal function"""
        return (1-exp(-self.lam*value))
    #sampling data with the given distribution   
    def sample(self, n, limits=None):
        res=np.empty(n)
        if limits:
           lower,upper=limits
           j=0
           while (j<n):
               def f(x):
           return self.cdf(x)-np.random.uniform(low=0,high=1,size=1)
           s=scipy.optimize.brenth(f,0,1e20)
           if s >= lower and s <= upper:
          res[j]=s
              j+=1
    else:
       r=np.random.uniform(low=0,high=1,size=n)
       for j in range(n):
               def f(x):
           return self.cdf(x)-r[j]
           s=scipy.optimize.brenth(f,0,1e20)
           res[j]=s
        return res

Using decorator also fails for instance:

@pm.stochastic(dtype=np.float, observed=False, trace=True)
def mass(name='MASS',value=1e14, m_l = 1e13,m_h=1e16):
        MassPrior=ExponentialPrior(1e-15)
        mlogp=MassPrior.logp(value,[m_l,m_h])
        print mlogp

I am wondering how could I define my own priors with pymc3 without using the available pymc distributions?