The Distributions.jl package does not calculate correctly the theoretical mode and/or modes under models like Binomial, Poisson and NegativeBinomial here some examples:

using Distributions
## Example 1
B = Binomial(3, 0.5)
v = collect(0:3);
hcat(v, pdf.(B, v)) # theoretical modes are 1 and 2
pdf(B, 1) == pdf(B, 2) # true, as expected
modes(B) # should return 1 and 2 not just 2
mode(B) # should return 1 not 2
# But surprisingly:
D = DiscreteNonParametric(v, pdf.(B, v))
pdf(D, 1) == pdf(D, 2) # true, as expected
modes(D) # 1 and 2, as expected
mode(D) # 1, as expected
## Example 2
P = Poisson(2.0)
x = collect(0:4);
hcat(x, pdf.(P, x)) # theoretical modes are 1 and 2
pdf(P, 1) == pdf(P, 2) # true, as expected
modes(P) # 1 and 2, as expected
mode(P) # should return 1 not 2
## Example 3
N = NegativeBinomial(3, 0.5)
y = collect(0:3);
hcat(y, pdf.(N, y)) # theoretical modes are 1 and 2
pdf(N, 1) == pdf(N, 2) # should be true
pdf(N, 1) - pdf(N, 2) # very small positive difference
pdf(N, 1) > pdf(N, 2) # so numerically the mode should be 1 but:
mode(N) # oops! 2 instead of 1
modes(N) # should return 1 and 2 not just 2

This is indeed an error. It falls back to a generic method because no method for modes is defined for NegativeBinomial. I can submit an issue if you do not have a github account.