Should we be putting this distribution in the math lib?
real sum_lfac = 0;
for (n in 1:N) sum_lfac += lgamma(y[n] + 1);
makes me realize we should really have a sum-of-lgamma-plus-1 function because it’s the log product of factorials and it’s so common for combinatoric normalization.
The above could be written as
real sum_lfac = sum(lgamma(to_vector(y) + 1));
but that’d actually be less efficient than the way Ben had it because of the extra copy in