ForwardDiff of logpdf of NegativeBinomial gives wrong results for p==1

[simsurace]

This was discovered and documented in #1568 via tests, but #1579 reverted those tests. I still think this is worth tracking:

using Distributions, ForwardDiff
using Distributions: digamma

mydiffp(r, p, k) = isone(p) && iszero(k) ? r : r/p - k/(1 - p)
mydiffr(r, p, k) = log(p) - inv(k + r) - digamma(r) + digamma(r + k + 1)

f1(_p) = logpdf(NegativeBinomial(1.0, _p), 0)
f2(_r) = logpdf(NegativeBinomial(_r, 1.0), 0)

ForwardDiff.derivative(f1, 1.0) ≈ mydiffp(1.0, 1.0, 0) # false
ForwardDiff.derivative(f2, 1.0) ≈ mydiffr(1.0, 1.0, 0) # false

[simsurace]

Should be solved by JuliaDiff/ForwardDiff.jl#481

[devmotion]

That is not required. I found a fix (and a bug on the way) without it.

[devmotion]

Basically,

• we can eliminate the special case branch alltogether by using xlogy and xlog1py
• we just have to fix the derivatives of these functions, defined in DiffRules and LogExpFunctions

I'll open the three required PRs.

[simsurace]

These tests pass with JuliaDiff/DiffRules.jl#85 and #1583.