I guess this is related to the fact that we cannot compute gradient for arbitrary PDFs?
In case someone find this helpful, for the 1D case (in v5), you can use pymc.Interpolated, but I haven’t been able to use the 2D PDFs.
Maybe there’s a way to implement the equivalent multivariate version from scipy.interpolate and implement the required logp, rand etc.