Drawing values from custom distribution

DensityDist should work fine for unobserved variables (and require nothing different). I do recommend installing the new version of PyMC (v4) where we improved several things related to the DensityDist

The API changed just slightly in that logp is a keyword only argument that must follow the distribution inputs.

import pymc as pm

def logp(value, mu, sigma):
  return pm.logp(pm.Normal.dist(mu, sigma), value)

with pm.Model() as m:
  # x = pm.Normal("x", 0, 1, shape=(2,))  # Completely equivalent
  x = pm.DensityDist("x", 0, 1, logp=logp, shape=(2,))
  y = pm.math.exp(x)  # f(x)
  z = pm.Poisson("z", mu=y, observed=[3, 4])

# Probability of variables at initial point
print(m.point_logps())  # {'x': -1.84, 'z': -6.97}