pl_pymc_model.py

import pymc as mc
from scipy.special import zeta
 
def _model(data, discrete=True, xmin=1.):
    alpha = mc.Exponential('alpha', 1. / 1.5)
 
    @mc.stochastic(observed=True)
    def custom_stochastic(value=data, alpha=alpha, xmin=xmin, discrete=discrete):
        value = value[value >= xmin]
        if discrete == True:
            return np.sum(np.log(value**-alpha / zeta(alpha,xmin)))
        else:
            return np.sum(np.log((alpha-1) * xmin**(alpha-1) * value**-alpha))
 
    return locals()