BUG (WIP): stats: Handle mu=0 correctly in poisson.pmf. Closes gh-2618.

*Work in progress; do not merge!*

This PR should close gh-2618.  The issue was fixed by changing `k*log(mu)` to `special.xlogy(k, mu)`, and by allowing mu=0 as an argument.