Defining a custom multivariate distribution

Hmmm, this seems to be fine?

with pm.Model(): 
    pm.NegativeBinomial('x', np.ones(10), np.ones(10), shape=10) 

Maybe try specifying the shape in pm.NegativeBinomial.dist(mu=mu, alpha=self.alpha, shape=...)