Hi all,
With many thanks to maxbiostat, we have a working implementation of the discrete Pareto in Stan with variable (fixed at arbitrary values) y_min. Some code, tests, and examples in both python and R are here:
Please note, as explained in the readme there, that the numerical approximation we are using for the Hurwitz zeta (based on quadrature integration) gives very poor results at alpha << 2. It seems straightforward to replace this with a proper implementation of the Hurwitz zeta using the Euler-Maclaurin summation formula which would relieve this problem, but after a few failed attempts I realized I would need some C++ help to get that working! There is at least one simple pre-existing C implementation (via cephes, as is used in scipy) that should work.