NegativeBinomial takes in 2 parameters - mu and alpha. The final shape of mu is (N,). So I am plugging in an array of N mu’s and would like to return the total log likelihood, i.e., the sum of the log likelihoods for each mu. Is there some way to define a distribution with a logp() that takes in an array (i.e., counts) rather than a scalar/constant?