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)))