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:
- What should I import so that I can properly inherit from
Continuous
? - In PyMC2, do I still need to stick to Theano notation?
Finally, how can I later tell my
Beta
distribution thatalpha
andbeta
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)
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.
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
In PyMC2, the trick is to put the
a
andb
parameters together:This notebook has a full example.