Custom priors in PyMC

Question

Say I want to put a custom prior on two variables a and b in PyMC, e.g.:

p(a,b)∝(a+b)^(−5/2)

(for the motivation behind this choice of prior, see this answer)

Can this be done in PyMC? If so how?

As an example, I would like to define such prior on a and b in the model below.

import pymc as pm

# ...
# Code that defines the prior: p(a,b)∝(a+b)^(−5/2)
# ...

theta   = pm.Beta("prior", alpha=a, beta=b)

# Binomials that share a common prior
bins = dict()
for i in xrange(N_cities):
    bins[i] = pm.Binomial('bin_{}'.format(i), p=theta,n=N_trials[i],  value=N_yes[i], observed=True)

mcmc = pm.MCMC([bins, ps])

Update

Following John Salvatier's advice, I try the following (note that I am in PyMC2, although I would be happy to switch to PyMC3), but my questions are:

1. What should I import so that I can properly inherit from Continuous?
2. In PyMC2, do I still need to stick to Theano notation?

3. Finally, how can I later tell my Beta distribution that alpha and beta have a prior from this multivariate distribution?

   import pymc.Multivariate.Continuous

   class CustomPrior(Continuous): """ p(a,b)∝(a+b)^(−5/2)

   :Parameters:
       None
   
   :Support:
       2 positive floats (parameters to a Beta distribution)
   """
   def __init__(self, mu, tau, *args, **kwargs):
       super(CustomPrior, self).__init__(*args, **kwargs)
   
   def logp(self, a,b):
   
   
       return np.log(math.power(a+b),-5./2)
   

Tellez

oh sorry, I thought you were in PyMC3 for some reason. You can do this in pymc2 but its different. I'll update my answer – Nanananak

Answer 1

In PyMC2, the trick is to put the a and b parameters together:

# Code that defines the prior: p(a,b)∝(a+b)^(−5/2)
@pm.stochastic
def ab(power=-2.5, value=[1,1]):
    if np.any(value <= 0):
        return -np.Inf
    return power * np.log(value[0]+value[1])

a = ab[0]
b = ab[1]

This notebook has a full example.

Notarize

Thanks! This is exactly what I was looking for. I have been looking for an answer to this question for a long time. Do you happen to know, by the way, how to do this in PyMC3? – Tellez

You are welcome. I've updated this with a bug fix (return log-probability) and the beginnings of a PyMC3 version. – Notarize

Answer 2

Yup! It's quite possible, and in fact quite straightforward.

If you're in PyMC 2, check out the documentation on the creation of stochastic variables.

@pymc.stochastic(dtype=int)
def switchpoint(value=1900, t_l=1851, t_h=1962):
    """The switchpoint for the rate of disaster occurrence."""
    if value > t_h or value < t_l:
        # Invalid values
        return -np.inf
    else:
        # Uniform log-likelihood
        return -np.log(t_h - t_l + 1)

If you're in PyMC 3, have a look at multivariate.py. Keep in mind the values passed in to init and logp are all theano variables, not numpy arrays. Is that enough to get you started?

For example, this is the Multivariate Normal distribution

class MvNormal(Continuous):
    """
    Multivariate normal

    :Parameters:
        mu : vector of means
        tau : precision matrix

    .. math::
        f(x \mid \pi, T) = \frac{|T|^{1/2}}{(2\pi)^{1/2}} \exp\left\{ -\frac{1}{2} (x-\mu)^{\prime}T(x-\mu) \right\}

    :Support:
        2 array of floats
    """
    def __init__(self, mu, tau, *args, **kwargs):
        super(MvNormal, self).__init__(*args, **kwargs)
        self.mean = self.median = self.mode = self.mu = mu
        self.tau = tau

    def logp(self, value):
        mu = self.mu
        tau = self.tau

        delta = value - mu
        k = tau.shape[0]

        return 1/2. * (-k * log(2*pi) + log(det(tau)) - dot(delta.T, dot(tau, delta)))

Nanananak

Thanks a lot @John. I think I understand the definition above, but I am trying to follow the steps without success. I have updated the OP with more specific questions. – Tellez

Thank you @John. This is very helpful. The only remaining question is how I can use the multivariate prior to get the prior on the individual variables (e.g. alpha and beta in a Binomial distribution, as in the code in the OP). How do I get the prior on the individual variables from the multivariate prior? – Tellez

Oh, like you want a multivariate prior on variables with different names? I don't think there's a way to do that. You can make a multivariate variable and index into it to get the components separately and then name them. Is that what you want? – Nanananak

Thanks @John. In the code example that I have at the top of my OP, I define theta = pm.Beta("prior", alpha=a, beta=b). What I want to do is define my prior on a and b as p(a,b)∝(a+b)^(−5/2). It looks like I can define p(a,b) in PyMC2 with @pymc.stochastic as you mentioned in your answer. What I would later need to pass a and b to the Beta call above. This is the part that I don't know how to do. – Tellez